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

using Graphs: SimpleEdge

import Base: *, //, /, ^, +, -, inv, one, zero

struct EquivariantClass
    rule::Expr
    func::Function
end

function Base.show( io::IO, ec::EquivariantClass )::Nothing
        #print(io, "EquivariantClass[ ... ]" )
        return nothing
end

#########################################
### Operations of equivariant classes ###
#########################################

### Products
function *( ec1::EquivariantClass, ec2::EquivariantClass )::EquivariantClass

    rule = quote
        $(ec1.rule)*$(ec2.rule)
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

function *( ec1::EquivariantClass, n::Number )::EquivariantClass

    rule = quote
        $(ec1.rule)*$(n)
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

function *( n::Number, ec1::EquivariantClass )::EquivariantClass

    return ec1*n
end

### Division
function //( ec1::EquivariantClass, ec2::EquivariantClass )::EquivariantClass

    rule = quote
        $(ec1.rule)//$(ec2.rule)
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

function //( ec1::EquivariantClass, n::Number )::EquivariantClass

    rule = quote
        $(ec1.rule)//$(n)
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

### Sums
function +( ec1::EquivariantClass, ec2::EquivariantClass )::EquivariantClass

    rule = quote
        $(ec1.rule)+$(ec2.rule)
    end

    return EquivariantClass( rule, eval(:(( g, c, w, s, m ) -> $rule )))
end

function +( ec1::EquivariantClass, n::Number )::EquivariantClass #it makes sene only if n==0

    rule = quote
        $(ec1.rule)+$(n)
    end

    return EquivariantClass( rule, eval(:(( g, c, w, s, m ) -> $rule )))
end
function +( n::Number, ec1::EquivariantClass )::EquivariantClass #it makes sene only if n==0

    return ec1+n
end
### Minus
function -( ec1::EquivariantClass, ec2::EquivariantClass )::EquivariantClass

    rule = quote
        $(ec1.rule)-$(ec2.rule)
    end

    return EquivariantClass( rule, eval(:(( g, c, w, s, m ) -> $rule )))
end
function -( ec1::EquivariantClass, n::Number )::EquivariantClass #it makes sene only if n==0

    rule = quote
        $(ec1.rule)-$(n)
    end

    return EquivariantClass( rule, eval(:(( g, c, w, s, m ) -> $rule )))
end
function -( n::Number, ec1::EquivariantClass )::EquivariantClass #it makes sene only if n==0

    rule = quote
        $(n)-$(ec1.rule)
    end

    return EquivariantClass( rule, eval(:(( g, c, w, s, m ) -> $rule )))
end
### Exponent
function ^( ec1::EquivariantClass, n::Number )::EquivariantClass

    rule = quote
        $(ec1.rule)^$(n)
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

### Constants
function one(ec1::EquivariantClass)::EquivariantClass

    rule = quote
        1
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function inv(ec1::EquivariantClass)::EquivariantClass

    rule = quote
        $(ec1.rule)^(-1)
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function zero(ec1::EquivariantClass)::EquivariantClass

    rule = quote
        0
    end

    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
#######################################
### List of all equivariant classes ###
#######################################
"""
    Incidency(r)

Equivariant class of the cycle parameterizing curves meeting a linear subspace of codimension `r`.
# Arguments
- `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 = Incidency(3)^2;
julia> AtiyahBottFormula(3,1,0,P);
Result: 1
julia> P = Incidency([2,2,3]);
julia> AtiyahBottFormula(3,1,0,P);
Result: 1
julia> P = Incidency([2,2])*Hypersurface(3);
julia> AtiyahBottFormula(3,3,0,P);
Result: 756
```
!!! warning "Attention!"

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

"""
function Incidency( r )::EquivariantClass

    rule = :(Incidency( g, c, w, s, $r ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function Incidency(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, r::Int64)::fmpq

    local p1::fmpq = fmpq(0); #the final result
    local temp1::fmpq = fmpq(1)
    local temp2::fmpq = fmpq(1)
    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)
            eq!(temp1, scalars[col[src(e)]])
            eq!(temp2, scalars[col[dst(e)]])
            pow_eq!(temp1, t)
            pow_eq!(temp2, r-t)
            mul_eq!(temp1, temp2)
            mul_eq!(temp1, d[e])
            add_eq!(p1, temp1)
            # 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 = fmpq(1)
    local temp1::fmpq = fmpq(1)

    for j in unique(r)
        eq!(temp1, Incidency(g, coloration, weights, scalars, j))
        pow_eq!(temp1, count(x->x==j, r))
        mul_eq!(p1, temp1)
        # p1 *= Hypersurface(g, coloration, weights, scalars, j)^count(x->x==j, b)
    end
    return p1

end

"""
    Hypersurface(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
- `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 = Hypersurface(5);
julia> AtiyahBottFormula(4,1,0,P);
Result: 2875
julia> P = Hypersurface([3,3]);
julia> AtiyahBottFormula(5,2,0,P);
Result: 423549//8
julia> P = Hypersurface(4)*Hypersurface(2);
julia> AtiyahBottFormula(5,3,0,P);
Result: 422690816//27
julia> P = Hypersurface(2)^4;
julia> AtiyahBottFormula(7,4,0,P);
Result: 25705160
```
!!! warning "Attention!"

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

"""
function Hypersurface( b )::EquivariantClass

    rule = :(Hypersurface( g, c, w, s, $b ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

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

    local p1::fmpq = fmpq(1)
    # local q1::fmpq = fmpq(1)
    local temp1::fmpq = fmpq(1)
    local temp2::fmpq = 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])
            eq!(temp1, scalars[col[src(e)]])
            eq!(temp2, scalars[col[dst(e)]])
            mul_eq!(temp1, alph)
            mul_eq!(temp2, b*d[e]-alph)
            add_eq!(temp1, temp2)
            div_eq!(temp1, d[e])
            mul_eq!(p1, temp1)
            # p1 *= (alph*scalars[col[src(e)]]+(b*d[e]-alph)*scalars[col[dst(e)]])//d[e]
        end
    end

    for v in vertices(g)
        eq!(temp1, scalars[col[v]])
        mul_eq!(temp1, b)
        pow_eq!(temp1, 1-length(all_neighbors(g, v)))
        mul_eq!(p1, temp1)
        #q1 *= (b*scalars[col[v]])^(1-length(all_neighbors(g, v)))
    end

    # mul_eq!(p1,q1)

    return p1

end

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

    local p1::fmpq = fmpq(1)
    local temp1::fmpq = fmpq(1)

    for j in unique(b)
        eq!(temp1, Hypersurface(g, coloration, weights, scalars, j))
        pow_eq!(temp1, count(x->x==j, b))
        mul_eq!(p1, temp1)
        # p1 *= Hypersurface(g, coloration, weights, scalars, j)^count(x->x==j, b)
    end
    return p1
end

"""
    Contact()

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.

# 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 = O1_i(1)^2*O1_i(2)^3*Contact();
julia> AtiyahBottFormula(3,1,2,P);
Result: 1
```
"""
function Contact()::EquivariantClass

    rule = :(Contact( g, c, w, s))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function Contact(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq})::fmpq

    local p1::fmpq = fmpq(1)
    # local q1::fmpq = fmpq(1)
    local temp1::fmpq = fmpq(1)
    local temp2::fmpq = 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)
            eq!(temp1, scalars[col[src(e)]])
            eq!(temp2, scalars[col[dst(e)]])
            mul_eq!(temp1, alph)
            mul_eq!(temp2, 2*d[e]-alph)
            add_eq!(temp1, temp2)
            div_eq!(temp1, d[e])
            mul_eq!(p1, temp1)
            # p1 *= (alph*scalars[col[src(e)]]+(2*d[e]-alph)*scalars[col[dst(e)]])//d[e]
        end
    end

    for v in vertices(g)
        eq!(temp1, scalars[col[v]])
        mul_eq!(temp1, 2)
        pow_eq!(temp1, length(all_neighbors(g, v))-1)
        mul_eq!(p1, temp1)
        # q1 *= (2*scalars[col[v]])^(length(all_neighbors(g, v))-1)
    end

    # mul_eq!(p1, q1)
    return p1

end

"""
    O1_i(j)

Equivariant class of the pull-back of ``\\mathcal{O}_{\\mathbb{P}^n}(1)`` with respect to the j-th evaluation map.
# Arguments
- `j::Int64`: the evaluation map.

# Example
The following Gromov-Witten invariants
```math
\\begin{aligned}
\\int_{\\overline{M}_{0,2}(\\mathbb{P}^{1},1)}\\mathrm{ev}_{1}^{*}\\mathcal{O}_{\\mathbb{P}^{1}}(1)\\cdot\\mathrm{ev}_{2}^{*}\\mathcal{O}_{\\mathbb{P}^{1}}(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 = O1_i(1)*O1_i(2);
julia> AtiyahBottFormula(1,1,2,P);
Result: 1
julia> P = O1_i(1)^2*Hypersurface(2);
julia> AtiyahBottFormula(3,1,1,P);
Result: 4
```
!!! warning "Attention!"

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

"""
function O1_i( j::Int64)::EquivariantClass
    rule = :(O1_i(g, c, w, s, m, $j ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

function O1_i(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Vector{Int64}, j::Int64)::fmpq
#function O1_i(col::Vector{UInt8}, scalars::Vector{fmpq}, mark::Vector{Int64}, j::Int64)::fmpq
    return scalars[col[mark[j]]]
end

"""
    O1()

Equivariant class of the pull-back of ``\\mathcal{O}_{\\mathbb{P}^n}(1)`` with respect to the product of all evaluation maps.

This function is equivalent to the product of the function `O1_i(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)^2 &= 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 = O1()^2;
julia> AtiyahBottFormula(2,3,8,P);
Result: 12
julia> P = O1()^2*Hypersurface(3);
julia> AtiyahBottFormula(3,2,1,P);
Result: 81
```

In order to remove `O1_i(j)` for some `j`, it is enough to divide by that function.

# Example
```julia-repl
julia> P = O1()//O1_i(1);
```
Here `P` is the product of all `O1_i(j)` where `j` runs from 2 to `m`.
"""
function O1()::EquivariantClass
    rule = :(O1( g, c, w, s, m ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end

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

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

    for t in 1:length(mark)
        mul_eq!(p1, scalars[col[mark[t]]])
        # p1 *= scalars[col[mark[t]]]
    end

    return p1
end

"""
    R1(k)

The equivariant class of the first derived functor of the pull-back of ``\\mathcal{O}_{\\mathbb{P}^n}(-k)``.
# Arguments
- `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 = R1(1)^2;
julia> AtiyahBottFormula(1,d,0,P);
Result: 1
```
!!! warning "Attention!"

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

"""
function R1( k )::EquivariantClass

    rule = :(R1( g, c, w, s, $k ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function R1(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, k::Int64)::fmpq

    local p1 = fmpq(1)
    # local q1 = fmpq(1)
    local temp1::fmpq = fmpq(1)
    local temp2::fmpq = 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)
            eq!(temp1, scalars[col[src(e)]])
            eq!(temp2, scalars[col[dst(e)]])
            mul_eq!(temp1, alph)
            mul_eq!(temp2, k*d[e]-alph)
            add_eq!(temp1, temp2)
            div_eq!(temp1, d[e])
            neg!(temp1)
            mul_eq!(p1, temp1)
            # p1 *= (-1)*(alph*scalars[col[src(e)]]+(k*d[e]-alph)*scalars[col[dst(e)]])//d[e]
        end
    end

    for v in vertices(g)
        eq!(temp1, scalars[col[v]])
        mul_eq!(temp1, -k)
        pow_eq!(temp1, length(all_neighbors(g, v))-1)
        mul_eq!(p1, temp1)
        # q1 *= (-k*scalars[col[v]])^(length(all_neighbors(g, v))-1)
    end

    # mul_eq!(p1, q1)
    return p1

end

"""
    Psi(a)

Equivariant class of the cycle of ``\\psi``-classes.
# Arguments
- `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 = O1_i(1)^5*O1_i(2)^2*Hypersurface(5)*Psi([1,0]);
julia> AtiyahBottFormula(6,2,2,P);
Result: 495000
julia> P = O1_i(1)^8*O1_i(2)^6*Hypersurface(7)*Psi(2);
julia> AtiyahBottFormula(10,2,2,P);
Result: 71804533752
julia> P = O1()^2*Psi(4);
julia> AtiyahBottFormula(2,2,1,P);
Result: 1//8
julia> P = Incidency(2)^4*O1_i(1)*(O1_i(1) + Psi(1));
julia> AtiyahBottFormula(2,2,1,P); #number of plane conics through four points and tangent to a line
Result: 2
julia> P = O1()^2*(Psi(3)*O1()+Psi(2)*O1()^2);
julia> AtiyahBottFormula(2,2,1,P);
Result: 1//8
julia> P = O1()*(Psi(7)*O1()+Psi(6)*O1()^2);
julia> AtiyahBottFormula(3,2,1,P);
Result: -5//16
```
!!! warning "Psi is singleton!"

    `Psi` cannot be multiplied by itself.
    ```julia-repl
    julia> P = O1()^2*Psi(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 = O1()^2*Psi(3)*Psi(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 = O1()^2*Psi(4);
    julia> AtiyahBottFormula(2,2,1,P);
    Result: 1//8
    ```
"""
function Psi( a )::EquivariantClass

    rule = :(Psi( g, c, w, s, m, $a ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function Psi(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Vector{Int64}, 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::Vector{Int64}, a::Vector{Int64})::fmpq

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

    local q1::fmpq = fmpq(1)
    local temp1::fmpq = fmpq(1)
    local Sum_ai::Int64
    local n::Int64
    local M::Int64
    local d = Dict(edges(g).=> weights) #assign weights to edges
    local inv_marks::Dict{Int64,Vector{Int64}} = invert_marks(mark, nv(g))

    for v in vertices(g)

        a_v = Int64[]
        for i in inv_marks[v]
            (i > length(a) || a[i] == 0) && continue
            push!(a_v, a[i])
        end

        Sum_ai = sum(a_v)
        Sum_ai == 0 && continue #if S contains only zeros, or it is empty, continue

        n = length(all_neighbors(g, v)) + length(inv_marks[v])

        n > 2 && Sum_ai > n - 3 && return fmpq(0)

        #If no previous condition holds, then n>1
        if n == 2 #necessary |S_v| == 1
            M = (-1)^a_v[1]
        else # n>2 and Sum_ai <= n - 3
            M = multinomial(n - 3 - Sum_ai, a_v...)
        end

        local s1 = fmpq(0)

        for w in all_neighbors(g, v)
            e = SimpleEdge(v,w)
            d_e = haskey(d,e) ? d[e] : d[reverse(e)]
            eq!(temp1, scalars[col[w]])
            neg!(temp1)
            add_eq!(temp1, scalars[col[v]])
            div_eq!(temp1, d_e)
            inv!(temp1)
            add_eq!(s1, temp1)
            # s1 += d_e//(scalars[col[v]]-scalars[col[w]])
        end
        pow_eq!(s1, -Sum_ai)
        mul_eq!(q1, s1)
        mul_eq!(q1, M)
        # s1 ^= -Sum_ai
        # q1 *= M*s1

    end

    return q1
end

"""
    Jet(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
- `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(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 = Incidency(2)^4*Jet(1,1);
julia> AtiyahBottFormula(2,2,1,P);
Result: 2
julia> P = Incidency(2)^4*(Jet(1,1)+O1()^2);
julia> AtiyahBottFormula(2,2,1,P);
Result: 3
julia> d=1;k=1; #for other values of d, change this line
julia> P = (O1()^2)//k*Jet(4*d-2,k);
julia> AtiyahBottFormula(3,d,1,P);   #The value of this integral does not depend on k, only on d
```
"""
function Jet( p, q )::EquivariantClass

    rule = :(Jet( g, c, w, s, m, $p, $q ))
    return EquivariantClass( rule, eval( :(( g, c, w, s, m ) -> $rule )))
end
function Jet(g::SimpleGraph, col::Vector{UInt8}, weights::Vector{Int64}, scalars::Vector{fmpq}, mark::Vector{Int64}, p::Int64, q::Int64)::fmpq

    local s1::fmpq = fmpq(0)
    local temp1::fmpq = fmpq(1)

    for h in 0:p
        eq!(temp1, O1_i(g, col, weights, scalars, mark, 1))
        mul_eq!(temp1, q)
        pow_eq!(temp1, p+1-h)
        mul_eq!(temp1, stirlings1(p+1, p+1-h))
        mul_eq!(temp1, Psi(g, col, weights, scalars, mark, [h]))
        add_eq!(s1, temp1)
        # 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::Vector{fmpq}`: the scalars.
- `m::Vector{Int64}`: the marks.

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

    local V::fmpq = fmpq(1)
    #local E::fmpq = fmpq(1)
    local temp1::fmpq = fmpq(1)
    local temp2::fmpq = fmpq(1)
    local q1::fmpq
    local p1::fmpq
    local s1::fmpq
    # 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)
        q1 = fmpq(1)
        for j in 1:max_col
            if j != col[src(e)] && j != col[dst(e)]
                for alph in 0:d[e]
                    eq!(temp1, scalars[col[src(e)]])
                    eq!(temp2, scalars[col[dst(e)]])
                    mul_eq!(temp1, alph)
                    mul_eq!(temp2, d[e]-alph)
                    add_eq!(temp1, temp2)
                    div_eq!(temp1, d[e])
                    sub!(temp1, scalars[j])
                    mul_eq!(q1, temp1)
                    #q1 *= ((alph*scalars[col[src(e)]]+(d[e]-alph)*scalars[col[dst(e)]])//d[e]-scalars[j])
                end
            end
        end

        #eq!(temp1, omega_inv[e])
        eq!(temp1, scalars[col[dst(e)]])
        neg!(temp1)
        add_eq!(temp1, scalars[col[src(e)]])
        inv!(temp1)
        mul_eq!(temp1, d[e])

        pow_eq!(temp1, 2*d[e])
        if isodd(d[e])
            neg!(temp1)
        end
        mul_eq!(q1, factorial(d[e])^2)
        #div_eq!(temp1, factorial(d[e])^2)
        div_eq!(temp1, q1)
        mul_eq!(V, temp1)
        #E *= ((omega_inv[e])^(2*d[e]))*((-1)^d[e])//(factorial(d[e])^2)//q1
    end

    for v in vertices(g)
        nghbrs = all_neighbors(g, v)
        p1 = fmpq(1)
        for j in 1:max_col
            if j != col[v]
                eq!(temp1, scalars[col[v]])
                sub!(temp1, scalars[j])
                mul_eq!(p1, temp1)
                #p1 *= scalars[col[v]]-scalars[j]
            end
        end
        pow_eq!(p1, length(nghbrs)-1)
        mul_eq!(V, p1)
        #p1 ^= length(nghbrs)-1

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

        for w in nghbrs
            eq!(temp1, scalars[col[w]])
            neg!(temp1)
            add_eq!(temp1, scalars[col[v]])
            inv!(temp1)
            e = SimpleEdge(v,w)
            d_e = haskey(d,e) ? d[e] : d[reverse(e)]
            mul_eq!(temp1, d_e)
            # e = SimpleEdge(v,w)
            add_eq!(s1, temp1)
            mul_eq!(V, temp1)
            # s1 += omega_inv[e]
            # r1 *= omega_inv[e]
        end

        pow_eq!(s1, length(nghbrs) + num_marks(mark,v) - 3)
        #s1 ^= length(nghbrs) + num_marks(mark,v) - 3
        #mul_eq!(V, p1)
        mul_eq!(V, s1)
        #V *= p1*s1*r1

    end

    #mul_eq!(V, E)
    return V
end

###############################################################
### List of the class useful for computing the rank (Cycle) ###
###############################################################
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:
using Documenter, IRTest

makedocs(;
    modules=[IRTest],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/IRTest.jl/blob/{commit}{path}#L{line}",
    sitename="IRTest.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
)

deploydocs(;
    repo="github.com/tkf/IRTest.jl",
)

ext: jl


text:
# conversion between Posit8, Posit16 and Posit32
Posit8(x::Posit16) = ccall((:p16_to_p8, SoftPositPath), Posit8, (Posit16,),x)
Posit32(x::Posit16) = ccall((:p16_to_p32, SoftPositPath), Posit32, (Posit16,),x)
Posit8(x::Posit32) = ccall((:p32_to_p8, SoftPositPath), Posit8, (Posit32,),x)
Posit16(x::Posit32) = ccall((:p32_to_p16, SoftPositPath), Posit16, (Posit32,),x)
Posit16(x::Posit8) = ccall((:p8_to_p16, SoftPositPath), Posit16, (Posit8,),x)
Posit32(x::Posit8) = ccall((:p8_to_p32, SoftPositPath), Posit32, (Posit8,),x)

# conversion from Posit8,16,32 to PX1
Posit8_1(x::Posit8) = ccall((:p8_to_pX1, SoftPositPath), Posit8_1, (Posit8,Int64),x,8)
Posit16_1(x::Posit8) = ccall((:p8_to_pX1, SoftPositPath), Posit16_1, (Posit8,Int64),x,16)
Posit24_1(x::Posit8) = ccall((:p8_to_pX1, SoftPositPath), Posit24_1, (Posit8,Int64),x,24)

Posit8_1(x::Posit16) = ccall((:p16_to_pX1, SoftPositPath), Posit8_1, (Posit16,Int64),x,8)
Posit16_1(x::Posit16) = ccall((:p16_to_pX1, SoftPositPath), Posit16_1, (Posit16,Int64),x,16)
Posit24_1(x::Posit16) = ccall((:p16_to_pX1, SoftPositPath), Posit24_1, (Posit16,Int64),x,24)

Posit8_1(x::Posit32) = ccall((:p32_to_pX1, SoftPositPath), Posit8_1, (Posit32,Int64),x,8)
Posit16_1(x::Posit32) = ccall((:p32_to_pX1, SoftPositPath), Posit16_1, (Posit32,Int64),x,16)
Posit24_1(x::Posit32) = ccall((:p32_to_pX1, SoftPositPath), Posit24_1, (Posit32,Int64),x,24)

# conversion from Posit8,16,32 to PX2
Posit8_2(x::Posit8) = ccall((:p8_to_pX2, SoftPositPath), Posit8_2, (Posit8,Int64),x,8)
Posit16_2(x::Posit8) = ccall((:p8_to_pX2, SoftPositPath), Posit16_2, (Posit8,Int64),x,16)
Posit24_2(x::Posit8) = ccall((:p8_to_pX2, SoftPositPath), Posit24_2, (Posit8,Int64),x,24)

Posit8_2(x::Posit16) = ccall((:p16_to_pX2, SoftPositPath), Posit8_2, (Posit16,Int64),x,8)
Posit16_2(x::Posit16) = ccall((:p16_to_pX2, SoftPositPath), Posit16_2, (Posit16,Int64),x,16)
Posit24_2(x::Posit16) = ccall((:p16_to_pX2, SoftPositPath), Posit24_2, (Posit16,Int64),x,24)

Posit8_2(x::Posit32) = ccall((:p32_to_pX2, SoftPositPath), Posit8_2, (Posit32,Int64),x,8)
Posit16_2(x::Posit32) = ccall((:p32_to_pX2, SoftPositPath), Posit16_2, (Posit32,Int64),x,16)
Posit24_2(x::Posit32) = ccall((:p32_to_pX2, SoftPositPath), Posit24_2, (Posit32,Int64),x,24)

# conversion from PX1 to Posit8,16,32
Posit8(x::Posit8_1) = ccall((:pX1_to_p8, SoftPositPath), Posit8, (Posit8_1,Int64),x,8)
Posit8(x::Posit16_1) = ccall((:pX1_to_p8, SoftPositPath), Posit8, (Posit16_1,Int64),x,16)
Posit8(x::Posit24_1) = ccall((:pX1_to_p8, SoftPositPath), Posit8, (Posit24_1,Int64),x,24)

Posit16(x::Posit8_1) = ccall((:pX1_to_p16, SoftPositPath), Posit16, (Posit8_1,Int64),x,8)
Posit16(x::Posit16_1) = ccall((:pX1_to_p16, SoftPositPath), Posit16, (Posit16_1,Int64),x,16)
Posit16(x::Posit24_1) = ccall((:pX1_to_p16, SoftPositPath), Posit16, (Posit24_1,Int64),x,24)

Posit32(x::Posit8_1) = ccall((:pX1_to_p32, SoftPositPath), Posit32, (Posit8_1,Int64),x,8)
Posit32(x::Posit16_1) = ccall((:pX1_to_p32, SoftPositPath), Posit32, (Posit16_1,Int64),x,16)
Posit32(x::Posit24_1) = ccall((:pX1_to_p32, SoftPositPath), Posit32, (Posit24_1,Int64),x,24)

# conversion from PX2 to Posit8,16,32
Posit8(x::Posit8_2) = ccall((:pX2_to_p8, SoftPositPath), Posit8, (Posit8_2,Int64),x,8)
Posit8(x::Posit16_2) = ccall((:pX2_to_p8, SoftPositPath), Posit8, (Posit16_2,Int64),x,16)
Posit8(x::Posit24_2) = ccall((:pX2_to_p8, SoftPositPath), Posit8, (Posit24_2,Int64),x,24)

Posit16(x::Posit8_2) = ccall((:pX2_to_p16, SoftPositPath), Posit16, (Posit8_2,Int64),x,8)
Posit16(x::Posit16_2) = ccall((:pX2_to_p16, SoftPositPath), Posit16, (Posit16_2,Int64),x,16)
Posit16(x::Posit24_2) = ccall((:pX2_to_p16, SoftPositPath), Posit16, (Posit24_2,Int64),x,24)

Posit32(x::Posit8_2) = ccall((:pX2_to_p32, SoftPositPath), Posit32, (Posit8_2,Int64),x,8)
Posit32(x::Posit16_2) = ccall((:pX2_to_p32, SoftPositPath), Posit32, (Posit16_2,Int64),x,16)
Posit32(x::Posit24_2) = ccall((:pX2_to_p32, SoftPositPath), Posit32, (Posit24_2,Int64),x,24)

# conversion between PX1 and PX2
Posit8_1(x::Posit8_2) = ccall((:pX2_to_pX1, SoftPositPath), Posit8_1, (Posit8_2,Int64),x,8)
Posit8_2(x::Posit8_1) = ccall((:pX1_to_pX2, SoftPositPath), Posit8_2, (Posit8_1,Int64),x,8)

Posit16_1(x::Posit16_2) = ccall((:pX2_to_pX1, SoftPositPath), Posit16_1, (Posit16_2,Int64),x,16)
Posit16_2(x::Posit16_1) = ccall((:pX1_to_pX2, SoftPositPath), Posit16_2, (Posit16_1,Int64),x,16)

Posit24_1(x::Posit24_2) = ccall((:pX2_to_pX1, SoftPositPath), Posit24_1, (Posit24_2,Int64),x,24)
Posit24_2(x::Posit24_1) = ccall((:pX1_to_pX2, SoftPositPath), Posit24_2, (Posit24_1,Int64),x,24)


text:
using NonparametricRegression
using LinearAlgebra
using StableRNGs, Distributions
using Test

rng = StableRNG(123)
vn = randn(rng,30)
vu = rand(rng,20)

# write tests here
@testset verbose=true "NonparametricRegression Tests" begin

N = 2500
xs = randn(rng,N)
ys = 2 .+ xs.^2 .+ 0.1 .* randn(rng,N)
xg = collect(-1:0.5:1)

@testset "Exported" begin
    @testset "Kernels" begin
        @testset "NormalKernel" begin
            nk1 = NormalKernel(1)
            @test nk1.(vn,0) ≈ pdf.(Normal(0,1),vn)

            nk01 = NormalKernel(0.1)
            @test nk01.(vn,0) ≈ pdf.(Normal(0,0.1),vn)

            @test bandwidth(nk1) == 1

            @test NonparametricRegression.scaling(nk1) ≈ 1/sqrt(2*pi)
        end

        @testset "UniformKernel" begin
            u1 = UniformKernel(1)
            @test u1.(vu,0) ≈ pdf.(Uniform(-1,1),vu)

            u01 = UniformKernel(0.1)
            @test u01.(vu,0) ≈ pdf.(Uniform(-0.1,0.1),vu)

            @test bandwidth(u1) == 1

            @test NonparametricRegression.scaling(u1) ≈ 1/2
        end

        @testset "EpanechnikovKernel" begin
            ek1 = EpanechnikovKernel(1)
            @test ek1.(vu,0) ≈ 3/4 .* (1 .- vu.^2)

            ek01 = EpanechnikovKernel(0.1)
            @test ek01.(vu,0) ≈ 3/4 .* (1 .- (vu./0.1).^2) ./ 0.1 .* (abs.(vu./0.1) .<= 1)

            @test bandwidth(ek1) == 1

            @test NonparametricRegression.scaling(ek1) ≈ 3/4
        end
    end

    @testset "localconstant" begin
        lc1 = localconstant(xs,ys,xg,1.0)
        @test lc1[3] ≈ 2 rtol=0.5

        loocv_h = leaveoneoutCV(xs,ys; method=:lc)
        # println(loocv_h)
        lc2 = localconstant(xs,ys,xg,loocv_h)
        @test lc2[3] ≈ 2 rtol=1e-2

        aicc_h = optimizeAICc(xs,ys; method=:lc)
        # println(aicc_h)
        lc3 = localconstant(xs,ys,xg,aicc_h)
        @test lc3[3] ≈ 2 rtol=1e-2

        # println([lc2 ;; lc3])

        nplc1 = npregress(xs,ys,xg,1.0; method=:lc)
        @test nplc1 ≈ lc1

        nplc2 = npregress(xs,ys,xg; method=:lc, bandwidthselection=:loocv)
        @test nplc2 ≈ lc2

        nplc3 = npregress(xs,ys,xg; method=:lc, bandwidthselection=:aicc)
        @test nplc3 ≈ lc3

    end

    @testset "locallinear" begin
        ll1 = locallinear(xs,ys,xg,1.0)
        @test ll1[3] ≈ 2 rtol=0.5

        loocv_h = leaveoneoutCV(xs,ys; method=:ll)
        # println(loocv_h)
        ll2 = locallinear(xs,ys,xg,loocv_h)
        @test ll2[3] ≈ 2 rtol=1e-2

        aicc_h = optimizeAICc(xs,ys; method=:ll)
        # println(aicc_h)
        ll3 = locallinear(xs,ys,xg,aicc_h)
        @test ll3[3] ≈ 2 rtol=1e-2

        # println([ll2 ;; ll3])

        npll1 = npregress(xs,ys,xg,1.0; method=:ll)
        @test npll1 ≈ ll1

        npll2 = npregress(xs,ys,xg; method=:ll, bandwidthselection=:loocv)
        @test npll2 ≈ ll2

        npll3 = npregress(xs,ys,xg; method=:ll, bandwidthselection=:aicc)
        @test npll3 ≈ ll3

        ll3ab = llalphabeta(xs,ys,xg,aicc_h)
        @test first(ll3ab) ≈ ll3
    end

    @testset "errors" begin
        @test_throws ErrorException npregress(xs,ys,xg,1; method=:hello)
        @test_throws ErrorException npregress(xs,ys,xg; method=:lc, bandwidthselection=:hello)
        @test_throws ErrorException npregress(xs,ys,xg; method=:hello, bandwidthselection=:loocv)
        @test_throws ErrorException npregress(xs,ys,xg; method=:hello, bandwidthselection=:aicc)
    end

end

@testset "Unexported" begin
    @testset "Utilities" begin
        @test NonparametricRegression.AICc(1,1,100) ≈ 0 + 1 + 2*(1+1)/(100-1-2)
    end

    @testset "KDE" begin
        @test NonparametricRegression.silvermanbw(vn) ≈ 0.4183944776432713

        kd = NonparametricRegression.densityestimator(xs,[0.0,1.0])

        @test kd ≈ pdf.(Normal(0,1),[0.0,1.0]) atol=1e-2
    end
end

end # end all tests


text:
@testset "Ideals in algebras" begin

  @testset "Left / Right" begin
    A = matrix_algebra(GF(3), 2)
    I = A[2]*A

    @test isleft_ideal(I) == false
    @test isright_ideal(I) == true
  end

  @testset "Quotients" begin
    Qx, x = FlintQQ["x"]
    # f = x^2 + 1
    # g = x^3 + 3x^2 + 5x - 5
    f2g3 = x^13 + 9x^12 + 44x^11 + 120x^10 + 205x^9 + 153x^8 + 32x^7 - 168x^6 - 5x^5 - 485x^4 + 500x^3 - 400x^2 + 375x - 125 # = f^2*g^3
    A = AlgAss(f2g3)
    J = radical(A)
    Q, AtoQ = quo(A, J)

    @test dim(Q) == 5
    @test iszero(radical(Q))

    @test iszero(AtoQ(zero(A)))
    @test AtoQ\zero(Q) in J
    @test isone(AtoQ(one(A)))
    @test mod(one(A), J) == mod(AtoQ\one(Q), J)

    for i = 1:5
      c = rand(A, -10:10)
      d = rand(A, -10:10)
      @test mod(AtoQ\(AtoQ(c)), J) == mod(c, J)
      @test mod(AtoQ\(AtoQ(d)), J) == mod(d, J)
      @test AtoQ(c + d) == AtoQ(c) + AtoQ(d)
      @test AtoQ(c*d) == AtoQ(c)*AtoQ(d)

      e = rand(Q, -10:10)
      f = rand(Q, -10:10)
      @test AtoQ(AtoQ\e) == e
      @test AtoQ(AtoQ\f) == f
      @test mod(AtoQ\(e + f), J) == mod((AtoQ\e + AtoQ\f), J)
      @test mod(AtoQ\(e*f), J) == mod(((AtoQ\e)*(AtoQ\f)), J)
    end

    J2 = J^2
    Q, AtoQ = quo(J, J2)

    @test dim(Q) == 5

    @test iszero(AtoQ(zero(A)))
    @test AtoQ\zero(Q) in J2
    @test_throws ErrorException AtoQ(one(A))

    for i = 1:5
      c = rand(J, -10:10)
      d = rand(J, -10:10)
      @test mod(AtoQ\(AtoQ(c)), J2) == mod(c, J2)
      @test mod(AtoQ\(AtoQ(d)), J2) == mod(d, J2)
      @test AtoQ(c + d) == AtoQ(c) + AtoQ(d)

      e = rand(Q, -10:10)
      f = rand(Q, -10:10)
      @test AtoQ(AtoQ\e) == e
      @test AtoQ(AtoQ\f) == f
      @test mod(AtoQ\(e + f), J2) == mod((AtoQ\e + AtoQ\f), J2)
    end
  end

end


text:
# start with julia echoserver.jl
using Sockets
server = listen(8080)               # 1
while true                          # 2
  conn = accept(server)             # 3
  @async begin                      # 4
    try                             # 5
      while true
        line = readline(conn)       # 6
        println(line)
        if chomp(line) == "S"       # 7
          println("Stopping TCP server...")
          close(conn)
          exit(0)
        else
          write(conn, line)         # 8 - here is the echo
        end
      end
    catch ex                        # 9
      print("connection lost with error $ex")
      close(conn)
    end
  end # end coroutine block
end
close(conn)


text:
module SumProductTest

using Test
using ForneyLab
using ForneyLab: generateId, addNode!, associate!, inferUpdateRule!, outboundType, isApplicable, setTargets!, messagePassingSchedule
using ForneyLab: SPClamp

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function MockNode(vars::Vector{Variable}; id=generateId(MockNode))
        n_interfaces = length(vars)
        self = new(id, Array{Interface}(undef, n_interfaces), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:n_interfaces
            self.i[idx] = self.interfaces[idx] = associate!(Interface(self), vars[idx])
        end

        return self
    end
end

@sumProductRule(:node_type     => MockNode,
                :outbound_type => Message{PointMass},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPMockOutPP)

@testset "@SumProductRule" begin
    @test SPMockOutPP <: SumProductRule{MockNode}
    @test isApplicable(SPMockOutPP, [Nothing, Message{PointMass}, Message{PointMass}])
    @test !isApplicable(SPMockOutPP, [Nothing, Message{PointMass}, Message{PointMass}, Message{PointMass}])
end

# Composite definition for inferUpdateRule! testset
@composite TestComposite (b,a) begin
    @RV z ~ GaussianMeanVariance(a, constant(1.0))
    b = constant(2.0) + z
end

@testset "inferUpdateRule!" begin
    FactorGraph()
    nd = MockNode([Variable(), constant(0.0), constant(0.0)])
    inferred_outbound_types = Dict(nd.i[2].partner => Message{PointMass}, nd.i[3].partner => Message{PointMass})

    entry = ScheduleEntry(nd.i[1], SumProductRule{MockNode})
    inferUpdateRule!(entry, entry.message_update_rule, inferred_outbound_types)

    @test entry.message_update_rule == SPMockOutPP

    # Internal msg passing tests
    FactorGraph()
    a = constant(0.0, id=:a)
    b = Variable(id=:b)
    tc = TestComposite(b, a)
    inferred_outbound_types = Dict(tc.i[:a].partner => Message{PointMass})
    entry = ScheduleEntry(tc.i[:b], SumProductRule{Nothing})

    inferUpdateRule!(entry, entry.message_update_rule, inferred_outbound_types)

    @test isdefined(entry, :internal_schedule)
    @test length(entry.internal_schedule) == 4
    @test entry.message_update_rule == entry.internal_schedule[end].message_update_rule
end

@testset "messagePassingSchedule" begin
    fg = FactorGraph()
    x = Variable()
    nd = MockNode([x, constant(0.0), constant(0.0)])

    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([x]))
    schedule = messagePassingSchedule(pf)

    @test length(schedule) == 3
    @test ScheduleEntry(nd.i[2].partner, SPClamp{Univariate}) in schedule
    @test ScheduleEntry(nd.i[3].partner, SPClamp{Univariate}) in schedule
    @test ScheduleEntry(nd.i[1], SPMockOutPP) in schedule
end

@testset "messagePassingAlgorithm" begin
    fg = FactorGraph()
    x = Variable()
    nd = MockNode([x, constant(0.0), constant(0.0)])

    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)

    @test isa(algo, InferenceAlgorithm)
end

end # module


text:
# rrdtool create .rrd

sensors = Dict("Karttula/Olohuone/T/" => "karttula_oh_t.rrd",
               "Karttula/Olohuone/RH/" => "karttula_oh_rh.rrd",
               "Karttula/Olohuone/P/" => "karttula_oh_p.rrd",
               "Karttula/WC/T/" => "karttula_wc_t.rrd",
               "Karttula/WC/RH/" => "karttula_wc_rh.rrd")

st, pr = open(`/usr/bin/mosquitto_sub -h kw.dy.fi -u valpuri -P ammu1001 -v -t "#"`, "r", STDOUT)

while true
    addr, value = split(readline(st), " ")
    if haskey(sensors, addr)
        file = sensors[addr]
        run(`/usr/bin/rrdtool update $file N:$value`)
        println("/usr/bin/rrdtool update $file N:$value")
    end
end

    # rrdtool create karttula.rrd \
    #         --step 30 \
    #         DS:karttula_oh_t:GAUGE:600:U:U \
    #         DS:karttula_oh_rh:GAUGE:600:U:U \
    #         DS:karttula_oh_p:GAUGE:600:U:U \
    #         DS:karttula_wc_t:GAUGE:600:U:U \
    #         DS:karttula_wc_rh:GAUGE:600:U:U \
    #         RRA:AVERAGE:0.5:1:120 \
    #         RRA:AVERAGE:0.5:10:288 \
    #         RRA:AVERAGE:0.5:30:672 \
    #         RRA:AVERAGE:0.5:120:720 \
    #         RRA:AVERAGE:0.5:360:2920 \
    #         RRA:AVERAGE:0.5:720:4380

    # One hour in 30s res, 120 = 60*60/30
    # One day in 5 min res, 5 mins =  5*60/30 = 10, 24*60/5 = 288
    # One week in 15 min res, 15 mins = 15*60/30 = 30, 7*24*60/15 = 672
    # One month in 1 hour res, 1 hours = 60*60/30 = 120, 30*24 = 720
    # One year in 3 hour res, 3 hours = 3*60*60/30 = 360, 365*24/3 = 2920
    # Three years in 6 hour res, 6 hours = 6*60*60/30 = 720, 3*365*24/6 = 4380

max_stars_repo_name: KajWiik/MQTT-monitor

max_issues_repo_name: KajWiik/MQTT-monitor

max_forks_repo_name: KajWiik/MQTT-monitor


text:
module Fail

export FailUnsat
import ..OmegaCore

"Fails on condition"
struct FailUnsatAlg end
const FailUnsat = FailUnsatAlg()

# Set default randalg to fail
OmegaCore.defrandalg(args...) = FailUnsat

"""
`n` samples of` ω ∈ ΩT` from from model.  Throws error if such that `pred`
"""
function randsample(rng,
                    ΩT::Type{OT},
                    pred,
                    n,
                    alg::FailUnsatAlg) where {OT <: Ω}
  ωsamples = ΩT[]
  ωsamples = Vector{ΩT}(undef, n)
  accepted = 0
  for i = 1:n
    ω = ΩT()
    issat = pred(tagrng(ω, rng))
    @inbounds ωsamples[i] = ω
    # lens(Loop, (ω = ω, accepted = accepted, p = float(issat), i = i))
  end
  ωsamples
end

function Base.rand(rng::AbstractRNG,
                   x::RandVar,
                   n::Integer,
                   alg::FailUnsatAlg;
                   ΩT::Type{OT} = defΩ(alg)) where {OT <: Ω}
  pred = Omega.indomain(x)
  ωsamples = rand(rng, ΩT, pred, n, alg)
  map(ω -> apl(x, ω), ωsamples)
end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Autogenerated wrapper script for Gumbo_jll for x86_64-apple-darwin
export libgumbo

JLLWrappers.@generate_wrapper_header("Gumbo")
JLLWrappers.@declare_library_product(libgumbo, "@rpath/libgumbo.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgumbo,
        "lib/libgumbo.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BenchmarkTools, ModelingToolkit, LinearAlgebra, SparseArrays

## ImplySymbolic

@variables x[1:30]

#S
Ind = [collect(1:12); 10; 10; 11; 11; 12; 12]
Jnd = [collect(1:12); 11; 12; 10; 12; 10; 11]
V = [[x[i] for i = 1:14]; x[13]; x[15]; x[14]; x[15]]
S = sparse(Ind, Jnd, V)

#F
Ind = collect(1:9)
Jnd = collect(1:9)
V = fill(1,9)
F = sparse(Ind, Jnd, V, 9, 12)

#A
Ind = collect(1:9)
Jnd = [fill(10, 3); fill(11, 3); fill(12, 3)]
V = [1; x[16]; x[17]; 1; x[18]; x[19]; 1; x[20]; x[21]]
A = sparse(Ind, Jnd, V, 12, 12)

start_val = vcat(
    fill(1, 10),
    [1, 1, 0.5, 0.5, 0.5],
    fill(0.5, 6),
    fill(1, 9))

new_val = fill(1.0, 31)

implysym_test = ImplySymbolic(A, S, F, x, start_val)
@benchmark implysym_test = ImplySymbolic(A, S, F, x, start_val)

implysym_test(new_val)
@benchmark implysym_test(new_val)

## large problem

@variables x[1:200]

Ind = [collect(1:30); 65; 66; collect(31:60); 20; 88; collect(61:90); 21; 50;
    92; 93; 93]
Jnd = [fill(91, 32); fill(92, 32); fill(93, 32); 91; 91; 92]
V = [x[i] for i = 1:99]
A = sparse(Ind, Jnd, V, 93, 93)


text:
"""
    calculateEnergy(r, r_j, r_dist, p_j, m, m_j, mu, mu_sum, ignore_H0)

Calculates the sum of kynetic and potential energy of the N-body hamiltonian system:

```math
\\mathcal{H}_{\\mathrm{N-body}}(\\textbf{q}, \\textbf{p}) = \\sum_{i=1}^{N} \\frac{\\textbf{p'}^{2}_{i}}{2 m_{i}'} - \\sum_{i=1}^{N} \\sum_{j=i+1}^{N} \\frac{Gm_{i}m_{j}}{\\mid \\textbf{q}_{i} - \\textbf{q}_{j} \\mid} \\\\
```

# Args

* `r`: NxD matrix where N is the number of bodies and D is the number of dimensions, each row contains the positions of the bodies in cartesian coordinates.
* `r_j`:  NxD matrix, it contains the positions of the bodies in jacobi coordinates.
* `r_dist`: array of N elements, in which distances from the bodies to the barycenter are saved.
* `p_j`: NxD matrix, it contains the momentum of the bodies in jacobi coordinates.
* `m`: array of N elements, where element number i contains the mass of i-th body.
* `m_j`: array of N elements, each containing the i-th jacobi mass (``m'_{i}``)
```math
m'_{1} = m_{sum_{N}} \\\\
m'_{i} = m_{i} \\frac{m_{sum_{i-1}}}{m_{sum_{i}}} \\\\
\\mathrm{,where } \\hspace{0.5cm} m_{sum_{i}} = \\sum_{j = 1}^{i} m_{j}
```
* `mu`: array of N elements, where element with index i contains the standard gravitational parameter of the i-th body
```math
mu_{i} = G \\cdot m_{j} \\\\
```
* `mu_{sum}`: array of N elements, where the standard gravitational parameter of the bodies from 1 to i is saved in the i-th element
```math
mu_{sum_{i}} = \\sum_{j = 1}^{i} G \\cdot m_{j} \\\\
```
* `ignore_H0`: boolean that determines if Hamiltonian ``\\mathcal{H_{0}}`` will be ignored (when calculating the energy).

# Returns

The energy of the ``\\mathcal{H}_{N-body}`` system.
"""
function calculateEnergy(r, r_j, r_dist, p_j, m, m_j, mu, mu_sum, ignore_H0)
    energy = 0.0

    # H_0
    if !ignore_H0
        energy += (dot(p_j[1,:],p_j[1,:]) / (2.0*m_j[1]))
    end

    # H_Kepler
    for i in 2:size(r)[1]
        energy += (dot(p_j[i,:],p_j[i,:]) / (2.0*m_j[i])) #- (mu_sum[i] * m_j[i] / r_dist[i])
    end
    #println(energy)

    # H_Interaction step 1
    #for i in 2:size(r)[1]
    #    energy += (mu_sum[i] * m_j[i] / r_dist[i])
    #end
    #println(energy)

    # H_Interaction step 2
    for i in 1:size(r)[1]
        for j in (i+1):size(r)[1]
            energy -= (mu[i] * m[j] / norm(r[i,:] - r[j,:]))
        end
    end

    return energy
end

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
module ScientificColourMaps

    using PyPlot, DelimitedFiles

    export acton, bamako, batlow, berlin, bilbao, broc, brocO, buda, cork, corkO, davos, devon, grayC,
            hawaii, imola, lajolla, lapaz, lisbon, nuuk, oleron, oslo, roma, romaO, tofino, tokyo,
            turku, vik, vikO

    include("load_cmaps.jl")

end


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

JLLWrappers.@generate_main_file_header("SFML")
JLLWrappers.@generate_main_file("SFML", UUID("a867eb61-6e20-512e-ac9f-cfafa65876cc"))
end  # module SFML_jll


text:
using FunctionalData, QuickShiftClustering
data = @p map unstack(1:10) (x->10*randn(2,1).+randn(2,1000)) | flatten
@time a = quickshift(data)
Profile.clear_malloc_data()
@time a = quickshift(data)
exit()


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

# This tests if actors generated with SlowActors and Actors
# can communicate.
#
using SlowActors, Test

raise(a, b) = a^b

query(act, a, b) = request!(act, a, b)

A = Actors.spawn(Func(raise))
B = spawn(Func(query, A))
@test request!(A, 2, 2) == 4
@test request!(B, 2, 2) == 4

A = spawn(Func(raise))
B = Actors.spawn(Func(query, A))
@test request!(A, 2, 2) == 4
@test request!(B, 2, 2) == 4


text:
using Distances
using Statistics
using DelimitedFiles
using LinearAlgebra
using QHull
using GeometricalPredicates
using Plots
using Measures
using ProgressBars
using Base

mutable struct Formants
  formants
  medians
  Formants(f, m) = new(f, m)
end

function Base.getindex(f::Formants, inds...)
  return f.formants[inds...]
end

mutable struct VowelSpace
  formants
  f1
  f2
  hull
  hullpts
  density
  ∂z∂f1
  ∂z∂f2
end

function Plots.plot(v::VowelSpace; standardize_axes=true, scale_densities=false, type=:density, kw...)

  f1vals = v.f1
  f2vals = v.f2

  if standardize_axes
    f1vals = f1vals .* v.formants.medians[1] .+ v.formants.medians[1]
    f2vals = f2vals .* v.formants.medians[2] .+ v.formants.medians[2]
  end

  vd = round(vdi(v, unstandardize=standardize_axes), digits=2)
  a = round(area(v, unstandardize=standardize_axes), digits=2)

  density = v.density
  if scale_densities
    density /= maximum(density)
  end

  # density is currently stored with f1 measured as changes in rows
  # while f2 is measured as changes in columns
  # need to transpose the matrix for f1 to be indexed with columns
  # and f2 to be indexed with rows

  if type==:density
    p = Plots.heatmap(f1vals, f2vals, density'; kw...)
  end

  hullpts = Vector()
  for vtx in v.hull.vertices
    push!(hullpts, v.hull.points[vtx, :]')
  end

  push!(hullpts, hullpts[1])

  hullpts = reduce(vcat, hullpts)

  if standardize_axes
    hullpts[:,1] = hullpts[:,1] .* v.formants.medians[1] .+ v.formants.medians[1]
    hullpts[:,2] = hullpts[:,2] .* v.formants.medians[2] .+ v.formants.medians[2]
  end

  if type==:density
    p = Plots.plot!(hullpts[:,1], hullpts[:,2], line=(3, :path), label="Threshold")
  elseif type==:area
    p = Plots.plot(hullpts[:,1], hullpts[:,2], line=(3, :path); kw...)
  end
  return p
end

function Plots.plot!(v::VowelSpace; standardize_axes=true, scale_densities=false, type=:density, kw...)
  p = plot(v, standardize_axes=standardize_axes, scale_densities=scale_densities, type=type; kw...)
  Plots.plot!(p)
end

"""
    vdi(v::VowelSpace; unstandardize=true)

Calculates the vowel dispersion index (vdi) for given `VowelSpace`, `v`. Effectively, it is a modification of the total variation of the vowel space density. Details are given in Kelley & Aalto (2019, Measuring the dispersion of density in head and neck cancer patients' vowel spaces: The vowel dispersion index, *Canadian Acoustics 47*(3), 114-115). Note that the vdi is a dimensionless number.

Args
=====

* `v` The `VowelSpace` for which to calculate the vdi
* `unstandardize` (keyword) Flag to convert the semi-normalized F1 and F2 in `v` to Hz (or whatever units were passed in when `v` was created).
"""
function vdi(v::VowelSpace; unstandardize=true)

  if unstandardize
    f1vals = v.f1[1:2] .* v.formants.medians[1] .+ v.formants.medians[1]
    Δf1 = diff(f1vals)
    f2vals = v.f2[1:2] .* v.formants.medians[2] .+ v.formants.medians[2]
    Δf2 = diff(f2vals)
  else
    Δf1 = diff(v.f1[1:2])
    Δf2 = diff(v.f2[1:2])
  end

  scaled_∂z∂f1 = v.∂z∂f1 ./ Δf1
  scaled_∂z∂f2 = v.∂z∂f2 ./ Δf2

  m = maximum(v.density)
  for (i, d) in enumerate(v.density)
    if d / m < 0.25
      scaled_∂z∂f1[i] = 0
      scaled_∂z∂f2[i] = 0
    end
  end

  return sum(sqrt.(scaled_∂z∂f1.^2  .+ scaled_∂z∂f2.^2))
end

"""
    area(v::VowelSpace; unstandardize=true)

Returns the area of the passed in `VowelSpace`, `v`.

Args
=====

* `v` The `VowelSpace` to calculate the area for
* `unstandardize` (keyword) Flag to convert the semi-normalized F1 and F2 in `v` to Hz (or whatever units were passed in when `v` was created). Will cause the convex hull to be re-calculated to determine the vowel space area in this new space

Returns
========

The area of the convex hull bounding the vowel space. Unit for area when F1 and F2 unstandardized back to Hz is Hz².
"""
function area(v::VowelSpace; unstandardize=true)
  if unstandardize
    hullpts = Vector()
    for vtx in v.hull.vertices
      push!(hullpts, v.hull.points[vtx, :]')
    end

    push!(hullpts, hullpts[1])

    hullpts = reduce(vcat, hullpts)
    hullpts[:,1] = hullpts[:,1] .* v.formants.medians[1] .+ v.formants.medians[1]
    hullpts[:,2] = hullpts[:,2] .* v.formants.medians[2] .+ v.formants.medians[2]
    return(chull(hullpts).area)
  else
    return v.hull.area
  end
end

"""
    VowelSpace(formants::Formants; temporalNorm=true)

Constructor for a `VowelSpace` based on `formants`. Though this function is not an exact implementation, see Story & Bunton (2017, Vowel space density as an indicator of speech performance, *J. Acoust. Soc. Am. 141*(5), EL458-EL464) for more details on calculating the vowel space and the vowel space density.

Args
======

* `formants` A `Formants` object from which to create the `VowelSpace` (and corresponding density)
* `temporalNorm` (keyword) Flag to divide the density counts by the number of discrete time steps in `formants`. Accounts for temporal differences that occur between recordings of different lengths, which will affect the density counts

Returns
=========

A `VowelSpace` object representing the calculated vowel space and density.
"""
function VowelSpace(formants::Formants; temporalNorm=true)

  minF1 = minimum(formants[:,1])
  maxF1 = maximum(formants[:,1])

  minF2 = minimum(formants[:,2])
  maxF2 = maximum(formants[:,2])

  dur = size(formants.formants, 1)

  densityGrid = createDensityGrid(formants.formants)

  local hullGrid

  if temporalNorm
    normGrid = densityGrid ./ size(formants.formants, 1)
    hullGrid = map(x -> x < 0.25 ? zero(x) : x, normGrid ./ maximum(normGrid))
  else
    normGrid = densityGrid ./ maximum(densityGrid)
    hullGrid = map(x -> x < 0.25 ? zero(x) : x, densityGrid)
  end

  f1vals = collect(range(minF1, maxF1, length=size(densityGrid, 1)))
  f2vals = collect(range(minF2, maxF2, length=size(densityGrid, 2)))

  densityPoints = Vector()
  for i in CartesianIndices(hullGrid)
    if hullGrid[i] > 0
      f1 = f1vals[i[1]]
      f2 = f2vals[i[2]]
      push!(densityPoints, [f1 f2])
    end
  end

  densityPoints = reduce(vcat, densityPoints)
  hull = chull(densityPoints)
  h1pts = Array{CartesianIndex{2},1}()

  h1poly = Polygon2D([Point(Float64.(hull.points[i,:])...) for i in hull.vertices]...)

  den_max = maximum(densityGrid)
  for i in CartesianIndices(densityGrid)
    densityGrid[i] / den_max < 0.25 && continue
    currX = f1vals[i[1]]
    currY = f2vals[i[2]]
    currPt = Point2D(currX, currY)
    if inpolygon(h1poly, currPt)
      push!(h1pts, i)
    end
  end

  # moving from row 1 to row 2, etc., changes f1
  # so, diff over dim1
  ∂z∂f1 = diff(densityGrid, dims=1)
  ∂z∂f1 = vcat(∂z∂f1, zeros(typeof(densityGrid[1]), size(densityGrid, 2))')

  # moving from column 1 to column 2, etc., changes f2
  # so, diff over dim2
  ∂z∂f2 = diff(densityGrid, dims=2)
  ∂z∂f2 = hcat(∂z∂f2, zeros(typeof(densityGrid[2]), size(densityGrid, 1)))

  return VowelSpace(formants, f1vals, f2vals, hull, h1pts, densityGrid, ∂z∂f1, ∂z∂f2)
end

function threeMdnNorm(x)
  filtered = [median([x[1], x[1], x[2]])]
  for i=1:(length(x)-2)
    push!(filtered, median(x[i:i+2]))
  end
  push!(filtered, median([x[end-1], x[end], x[end]]))
  return filtered
end

function fiveSpanMovingAvg(x)
  filtered = [x[1]]
  push!(filtered, mean(x[1:3]))
  for i=3:length(x)-2
    push!(filtered, mean(x[i-2:i+2]))
  end
  push!(filtered, mean(x[end-2:end]))
  push!(filtered, x[end])
end

function createDensityGrid(formants, gridInterval=0.01, densityRadius=0.05)

  pointsTree = TextVPTree([formants[i,:] for i=1:size(formants, 1)], euclidean)

  minF1 = minimum(formants[:,1])
  maxF1 = maximum(formants[:,1])

  minF2 = minimum(formants[:,2])
  maxF2 = maximum(formants[:,2])

  nXvals = round(Int64, (maxF1 - minF1) / gridInterval)
  nYvals = round(Int64, (maxF2 - minF2) / gridInterval)

  densityGrid = zeros(nXvals, nYvals)

  for i in 1:size(densityGrid, 1)
    xVal = minF1 + (i-1)*gridInterval
    for j=1:size(densityGrid, 2)
      yVal = minF2 + (j-1)*gridInterval
      densityGrid[i,j] = length(radiusSearch(pointsTree, [xVal, yVal], densityRadius))
    end
  end
  return densityGrid
end

"""
    Formants(fname; mdnFilt=true, avgFilt=true, outlierFilt=false)

Constructs a `Formants` object from tabular formant data contained in the file `fname` refers to. See Story & Bunton (2017, Vowel space density as an indicator of speech performance, *J. Acoust. Soc. Am. 141*(5), EL458-EL464) for more details on processing the formant data.

Args
=====

* `fname` A `String` containing the name of the file of formant data to read in. The file is expected to be tab-delimited, where the first column has F1 measurements, and the second column has F2 measurements.
* `mdnFilt` (keyword) Flag to perform a moving median filter on the formant data (current window size is set to 3)
* `avgFilt` (keyword) Flag to perform a moving mean filter on the formant data (current window size is set to five)
* `outlierFilt` (keyword) Flag to remove outlying formant observations (current threshold is set to remove rows where F1 >= 800 or F2 >= 2300); not specified in Story & Bunton (2017)

Returns
========

A `Formants` object containing the array of formant values, as well as the values needed to undo the median transformation that is applied to the formant data.
"""
function Formants(fname; mdnFilt=true, avgFilt=true, outlierFilt=false)
  open(fname, "r") do f
    formants = DelimitedFiles.readdlm(f, '\t')
  end

  if mdnFilt
    formants = mapslices(threeMdnNorm, formants, dims=1)
  end

  if avgFilt
    formants = mapslices(fiveSpanMovingAvg, formants, dims=1)
  end

  if outlierFilt
    formants = formants[findall(x -> x < 800, formants[:,1]),:]
    formants = formants[findall(x -> x < 2300, formants[:,2]),:]
  end

  md1 = median(formants[:,1])
  md2 = median(formants[:,2])
  formants = mapslices(x -> (x .- median(x)) ./ median(x), formants, dims=1)
  return Formants(formants, [md1, md2])
end


text:
using Random, Printf, JLD
const seed = 420  # Random seed to use for all pseudorandom number generators.

mutable struct JLDOutputWriter{N} <: OutputWriter
                 dir :: String
              prefix :: String
           fieldsets :: NTuple{N, Symbol}
    output_frequency :: Int
             padding :: Int
end

ext(::JLDOutputWriter) = ".jld"

function JLDOutputWriter(; dir=".", prefix="", fieldsets=(:velocities, :tracers, :G), frequency=1, padding=9)
    return JLDOutputWriter(dir, prefix, fieldsets, frequency, padding)
end

filename(iter, fw::JLDOutputWriter) = joinpath(fw.dir, fw.prefix * lpad(iter, fw.padding, "0") * ext(fw))

function Oceananigans.write_output(model, fw::JLDOutputWriter)
    filepath = filename(model.clock.iteration, fw)
    write_jld_output(Tuple(getproperty(model, s) for s in fw.fieldsets), filepath)
end

function Oceananigans.read_output(name, iter, fw::JLDOutputWriter)
    filepath = filename(iter, fw)
    data = load(filepath, name)
    return data
end

function write_jld_output(sets, filepath)
    allvars = Dict{String, Any}()
    for s in sets
        svars = Dict((String(fld), Array(getproperty(s, fld).data)) for fld in propertynames(s))
        merge!(allvars, svars)
    end

    save(filepath, allvars)

    return nothing
end

function run_thermal_bubble_regression_tests(arch)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=arch, ν=4e-2, κ=4e-2)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    nc_writer = NetCDFOutputWriter(dir=".",
                                   prefix="thermal_bubble_regression_",
                                   frequency=10, padding=2)

    # Uncomment to include a NetCDF output writer that produces the regression.
    # push!(model.output_writers, nc_writer)

    time_step!(model, 10, Δt)

    u = read_output(nc_writer, "u", 10)
    v = read_output(nc_writer, "v", 10)
    w = read_output(nc_writer, "w", 10)
    T = read_output(nc_writer, "T", 10)
    S = read_output(nc_writer, "S", 10)

    field_names = ["u", "v", "w", "T", "S"]
    fields = [model.velocities.u.data, model.velocities.v.data, model.velocities.w.data, model.tracers.T.data, model.tracers.S.data]
    fields_gm = [u, v, w, T, S]
    for (field_name, φ, φ_gm) in zip(field_names, fields, fields_gm)
        φ_min = minimum(Array(φ) - φ_gm)
        φ_max = maximum(Array(φ) - φ_gm)
        φ_mean = mean(Array(φ) - φ_gm)
        φ_abs_mean = mean(abs.(Array(φ) - φ_gm))
        φ_std = std(Array(φ) - φ_gm)
        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, φ_min, φ_max, φ_mean, φ_abs_mean, φ_std))
    end

    # Now test that the model state matches the regression output.
    @test all(Array(model.velocities.u.data) .≈ u)
    @test all(Array(model.velocities.v.data) .≈ v)
    @test all(Array(model.velocities.w.data) .≈ w)
    @test all(Array(model.tracers.T.data)    .≈ T)
    @test all(Array(model.tracers.S.data)    .≈ S)
end

function run_deep_convection_regression_tests()
    Nx, Ny, Nz = 32, 32, 16
    Lx, Ly, Lz = 2000, 2000, 1000
    Δt = 20

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2)

    function cooling_disk(grid, u, v, w, T, S, i, j, k)
        if k == 1
            x = i*grid.Δx
            y = j*grid.Δy
            r² = (x - grid.Lx/2)^2 + (y - grid.Ly/2)^2
            if r² < 600^2
                return -4.5e-6
            else
                return 0
            end
        else
            return 0
        end
    end

    model.forcing = Forcing(nothing, nothing, nothing, cooling_disk, nothing)

    rng = MersenneTwister(seed)
    model.tracers.T.data[:, :, 1] .+= 0.01*rand(rng, Nx, Ny)

    nc_writer = NetCDFOutputWriter(dir=".",
                                   prefix="deep_convection_regression_",
                                   frequency=10, padding=2)

    # Uncomment to include a NetCDF output writer that produces the regression.
    # push!(model.output_writers, nc_writer)

    time_step!(model, 10, Δt)

    u = read_output(nc_writer, "u", 10)
    v = read_output(nc_writer, "v", 10)
    w = read_output(nc_writer, "w", 10)
    T = read_output(nc_writer, "T", 10)
    S = read_output(nc_writer, "S", 10)

    field_names = ["u", "v", "w", "T", "S"]
    fields = [model.velocities.u.data, model.velocities.v.data, model.velocities.w.data, model.tracers.T.data, model.tracers.S.data]
    fields_gm = [u, v, w, T, S]
    for (field_name, φ, φ_gm) in zip(field_names, fields, fields_gm)
        φ_min = minimum(Array(φ) - φ_gm)
        φ_max = maximum(Array(φ) - φ_gm)
        φ_mean = mean(Array(φ) - φ_gm)
        φ_abs_mean = mean(abs.(Array(φ) - φ_gm))
        φ_std = std(Array(φ) - φ_gm)
        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, φ_min, φ_max, φ_mean, φ_abs_mean, φ_std))
    end

    # Now test that the model state matches the regression output.
    @test_skip all(model.velocities.u.data .≈ u)
    @test_skip all(model.velocities.v.data .≈ v)
    @test_skip all(model.velocities.w.data .≈ w)
    @test_skip all(model.tracers.T.data    .≈ T)
    @test_skip all(model.tracers.S.data    .≈ S)
end

function run_rayleigh_benard_regression_test(arch)

    #
    # Parameters
    #
          α = 2                 # aspect ratio
          n = 1                 # resolution multiple
         Ra = 1e6               # Rayleigh number
    Nx = Ny = 8n * α            # horizontal resolution
    Lx = Ly = 1.0 * α           # horizontal extent
         Nz = 16n               # vertical resolution
         Lz = 1.0               # vertical extent
         Pr = 0.7               # Prandtl number
          a = 1e-1              # noise amplitude for initial condition
         Δb = 1.0               # buoyancy differential

    # Rayleigh and Prandtl determine transport coefficients
    ν = sqrt(Δb * Pr * Lz^3 / Ra)
    κ = ν / Pr

    #
    # Model setup
    #

    model = Model(
         arch = arch,
            N = (Nx, Ny, Nz),
            L = (Lx, Ly, Lz),
            ν = ν,
            κ = κ,
          eos = LinearEquationOfState(βT=1., βS=0.),
    constants = PlanetaryConstants(g=1., f=0.)
    )

    # Constant buoyancy boundary conditions on "temperature"
    model.boundary_conditions.T.z.top = BoundaryCondition(Value, 0.0)
    model.boundary_conditions.T.z.bottom = BoundaryCondition(Value, Δb)

    # Force salinity as a passive tracer (βS=0)
    S★(x, z) = exp(4z) * sin(2π/Lx * x)
    FS(grid, u, v, w, T, S, i, j, k) = 1/10 * (S★(grid.xC[i], grid.zC[k]) - S[i, j, k])
    model.forcing = Forcing(FS=FS)

    ArrayType = typeof(model.velocities.u.data)
    Δt = 0.01 * min(model.grid.Δx, model.grid.Δy, model.grid.Δz)^2 / ν

    spinup_steps = 1000
      test_steps = 100

    prefix = "data_rayleigh_benard_regression_"
    outputwriter = JLDOutputWriter(dir=".", prefix=prefix, frequency=test_steps)

    #
    # Initial condition and spinup steps for creating regression test data
    #

    #=
    ξ(z) = a * rand() * z * (Lz + z) # noise, damped at the walls
    b₀(x, y, z) = (ξ(z) - z) / Lz

    x, y, z = model.grid.xC, model.grid.yC, model.grid.zC
    x, y, z = reshape(x, Nx, 1, 1), reshape(y, 1, Ny, 1), reshape(z, 1, 1, Nz)

    model.tracers.T.data .= ArrayType(b₀.(x, y, z))

    println("Spinning up... ")

    @time begin
        time_step!(model, spinup_steps-test_steps, Δt)
        push!(model.output_writers, outputwriter)
    end

    time_step!(model, 2test_steps, Δt)
    =#

    #
    # Regression test
    #

    # Load initial state
    u₀ = read_output("u",  spinup_steps, outputwriter)
    v₀ = read_output("v",  spinup_steps, outputwriter)
    w₀ = read_output("w",  spinup_steps, outputwriter)
    T₀ = read_output("T",  spinup_steps, outputwriter)
    S₀ = read_output("S",  spinup_steps, outputwriter)

    Gu = read_output("Gu", spinup_steps, outputwriter)
    Gv = read_output("Gv", spinup_steps, outputwriter)
    Gw = read_output("Gw", spinup_steps, outputwriter)
    GT = read_output("GT", spinup_steps, outputwriter)
    GS = read_output("GS", spinup_steps, outputwriter)

    model.velocities.u.data .= ArrayType(u₀)
    model.velocities.v.data .= ArrayType(v₀)
    model.velocities.w.data .= ArrayType(w₀)
       model.tracers.T.data .= ArrayType(T₀)
       model.tracers.S.data .= ArrayType(S₀)

    model.G.Gu.data .= ArrayType(Gu)
    model.G.Gv.data .= ArrayType(Gv)
    model.G.Gw.data .= ArrayType(Gw)
    model.G.GT.data .= ArrayType(GT)
    model.G.GS.data .= ArrayType(GS)

    model.clock.iteration = spinup_steps
    model.clock.time = spinup_steps * Δt

    # Step the model forward and perform the regression test
    time_step!(model, test_steps, Δt)

    u₁ = read_output("u", spinup_steps + test_steps, outputwriter)
    v₁ = read_output("v", spinup_steps + test_steps, outputwriter)
    w₁ = read_output("w", spinup_steps + test_steps, outputwriter)
    T₁ = read_output("T", spinup_steps + test_steps, outputwriter)
    S₁ = read_output("S", spinup_steps + test_steps, outputwriter)

    field_names = ["u", "v", "w", "T", "S"]
    fields = [model.velocities.u.data, model.velocities.v.data, model.velocities.w.data, model.tracers.T.data, model.tracers.S.data]
    fields_gm = [u₁, v₁, w₁, T₁, S₁]
    for (field_name, φ, φ_gm) in zip(field_names, fields, fields_gm)
        φ_min = minimum(Array(φ) - φ_gm)
        φ_max = maximum(Array(φ) - φ_gm)
        φ_mean = mean(Array(φ) - φ_gm)
        φ_abs_mean = mean(abs.(Array(φ) - φ_gm))
        φ_std = std(Array(φ) - φ_gm)
        @info(@sprintf("Δ%s: min=%.6g, max=%.6g, mean=%.6g, absmean=%.6g, std=%.6g\n",
                       field_name, φ_min, φ_max, φ_mean, φ_abs_mean, φ_std))
    end

    @test all(Array(model.velocities.u.data) .≈ u₁)
    @test all(Array(model.velocities.v.data) .≈ v₁)
    @test all(Array(model.velocities.w.data) .≈ w₁)
    @test all(Array(model.tracers.T.data)    .≈ T₁)
    @test all(Array(model.tracers.S.data)    .≈ S₁)

    return nothing
end


text:
# using DSGE, ModelConstructors, Dates, CSV, FileIO, Random, JLD2, Test, Nullables
# Load data to use for tests and prep model objet
path = dirname(@__FILE__)
data = JLD2.jldopen("$path/../reference/load_data_out.jld2", "r") do file
    read(file, "data")
end
fred = CSV.read("$path/../reference/fred_160812.csv")
custom_settings = Dict{Symbol, Setting}(
    :data_vintage             => Setting(:data_vintage, "160812"),
    :cond_vintage             => Setting(:cond_vintage, "160812"),
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :date_conditional_end     => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :n_anticipated_shocks     => Setting(:n_anticipated_shocks, 6))

m = Model990(custom_settings = custom_settings, testing = true)
m <= Setting(:rate_expectations_source, :ois)
hpfilter_out = JLD2.jldopen("$path/../reference/transformations_out.jld2", "r") do file
    read(file, "hpfilter_out")
end
fred[!,:filtered_population] = hpfilter_out
@testset "Check quarter and percapita helper functions" begin
    @test annualtoquarter(4.) == 1.
    @test quartertoannual(1.) == 4.
    @test quartertoannualpercent(1.) == 400.
    @test @test_matrix_approx_eq (nominal_to_real(:GDP, fred, deflator_mnemonic = :GDPDEF)) (fred[!,:GDP] ./ fred[!,:GDPDEF])
    @test (percapita(:GDP, fred, :CNP16OV)) ≈ (fred[!,:GDP] ./ fred[!,:CNP16OV])
    @test @test_matrix_approx_eq (percapita(m, :GDP, fred))  (fred[!,:GDP] ./ fred[!,:filtered_population])
    m <= Setting(:hpfilter_population, false)
    @test @test_matrix_approx_eq (percapita(m, :GDP, fred))  (fred[!,:GDP] ./ fred[!,:CNP16OV])
    m <= Setting(:population_mnemonic, Nullable())
    @test_throws ErrorException percapita(m, :GDP, fred)
    @test DSGE.lag(fred[!,:GDP], 10)[11:end] ==  fred[!,:GDP][1:end-10]
    hpfilter_test, ~ = hpfilter(fred[!,:CNP16OV], 1600.)
    @test @test_matrix_approx_eq hpfilter_test hpfilter_out
    @test @test_matrix_approx_eq difflog([1., 3., 1.]) vec([missing, log(3), -log(3)])
    @test @test_matrix_approx_eq oneqtrpctchange([1., 3., 1.])  100 * vec([missing, log(3), -log(3)])
end

@testset "Check loggrowth and loglevel transforms" begin
    @test loggrowthtopct([1., 3.]) == 100 .* (exp.([.01, .03]) .- 1.)
    @test loggrowthtopct_percapita([1., 3.], [1., 3.]) == 100. .* ((exp.([.01, .03]) .*
                                                                   exp.([1., 3.]).^4) .- 1.)
    @test loggrowthtopct_percapita([1. 3.; 1. 3.], [1., 3.]) == (100. .* ((exp.([.01, .03]) .*
                                                                           exp.([1., 3.]).^4) .- 1.))' .* ones(2)
    @test_throws AssertionError loggrowthtopct_percapita([1., 3.], [1., 3., 4.])
    @test loggrowthtopct_annualized([1., 3.]) == 100. .* (exp.([.01, .03]).^4 .- 1.)
    @test loggrowthtopct_annualized_percapita([2., 3.], [1.1, 1.2]) == 100. * (exp.([2., 3.] ./ 100. .+  [1.1, 1.2]).^4 .- 1.)
    @test loggrowthtopct_annualized_percapita([2. 3.; 2. 3.], [1.1, 1.2]) == ones(2) * (100. * (exp.([2., 3.] ./ 100. .+
                                                                                [1.1, 1.2]).^4 .- 1.))'
    @test_throws AssertionError loggrowthtopct_annualized_percapita([1., 3.], [1., 3., 4.])

    @test isnan(logleveltopct_annualized_approx([.1, .2])[1])
    @test logleveltopct_annualized_approx([.1, .2], 0.) == [.4, .4]
    @test sum(isnan.(logleveltopct_annualized_approx([.1 .2; 1. 2.]))) == 2
    @test logleveltopct_annualized_approx([.1 .2; 1. 2.], 0.) == [.4 .4; 4 4]
end

@testset "Check retrieval functions for transforms" begin
    @test DSGE.get_nopop_transform(loggrowthtopct_annualized_percapita) == loggrowthtopct_annualized
    @test DSGE.get_nopop_transform(logleveltopct_annualized_percapita) == logleveltopct_annualized
    @test DSGE.get_nopop_transform(sum) == sum
    @test DSGE.get_irf_transform(loggrowthtopct_annualized_percapita) == loggrowthtopct_annualized
    @test DSGE.get_irf_transform(logleveltopct_annualized_percapita) == logleveltopct_annualized
    @test DSGE.get_irf_transform(sum) == sum
    @test DSGE.get_transform4q(loggrowthtopct_annualized_percapita) == DSGE.loggrowthtopct_4q_percapita
    @test DSGE.get_transform4q(loggrowthtopct_annualized) == DSGE.loggrowthtopct_4q
    @test DSGE.get_transform4q(logleveltopct_annualized_percapita) == DSGE.logleveltopct_4q_percapita
    @test DSGE.get_transform4q(logleveltopct_annualized) == DSGE.logleveltopct_4q
    @test DSGE.get_transform4q(identity) == identity
    @test_throws ErrorException DSGE.get_transform4q(sum)

    @test DSGE.get_scenario_transform(loggrowthtopct_annualized_percapita) == quartertoannual
    @test DSGE.get_scenario_transform(loggrowthtopct_annualized) == quartertoannual
    @test DSGE.get_scenario_transform(quartertoannual) == quartertoannual
    @test DSGE.get_scenario_transform(logleveltopct_annualized_percapita) == logleveltopct_annualized_approx
    @test DSGE.get_scenario_transform(logleveltopct_annualized) == logleveltopct_annualized_approx
    @test DSGE.get_scenario_transform(identity) == identity
    @test_throws ErrorException DSGE.get_scenario_transform(sum)
end

@testset "Check 4q transforms" begin
    @test DSGE.prepend_data([1., 2.], [1., 2.]) == [1., 2., 1., 2.]
    @test DSGE.prepend_data([1. 2.; 1. 2.], [1.; 2. ]) == [1. 2. 1. 2. ; 1. 2. 1. 2.]
    @test sum(isnan.(DSGE.loggrowthtopct_4q([1., 2.]))) == 2
    @test DSGE.loggrowthtopct_4q([1., 2.], [1., 1., 1.]) == 100. .* (exp.([4., 5.] ./ 100) .- 1.)
    @test sum(isnan.(DSGE.loggrowthtopct_4q([1. 2.; 1. 2.]))) == 4
    @test DSGE.loggrowthtopct_4q([1. 2.; 1. 2.], [1., 1., 1.]) == 100. .* (exp.([4. 5.; 4. 5.] ./ 100) .- 1.)
    @test sum(isnan.(DSGE.loggrowthtopct_4q_percapita([1.], [1., 1., 1., 1.]))) == 1
    @test DSGE.loggrowthtopct_4q_percapita([1.], [1., 1., 1., 1.], [1., 1., 1.]) == 100. .* (exp.([4.] ./ 100. .+ [4.]) .- 1.)
    @test sum(isnan.(DSGE.loggrowthtopct_4q_percapita([1. 1.; 1. 1.], [1., 1., 1., 1., 1.]))) == 4
    @test DSGE.loggrowthtopct_4q_percapita([1. 1.; 1. 1.], [1., 1., 1., 1., 1.], [1., 1., 1.]) == ones(2) * 100. .* (exp.([4.; 4.] ./
                                                                                100. .+ [4.; 4.]) .- 1.)'
    @test sum(isnan.(DSGE.logleveltopct_4q([1., 2., 3., 4.]))) == 4
    @test DSGE.logleveltopct_4q([1., 2., 3., 4. ], [1., 1., 1., 1.]) == 100. .* (exp.([0., 1., 2., 3.] ./ 100) .- 1.)
    @test sum(isnan.(DSGE.logleveltopct_4q([1. 2. 3. 4. ; 1. 2. 3. 4. ]))) == 8
    @test DSGE.logleveltopct_4q([1. 2. 3. 4. ; 1. 2. 3. 4.], [1., 1., 1., 1.]) == 100. .* (exp.([0. 1. 2. 3.; 0. 1. 2. 3.] ./ 100) .- 1.)
    @test sum(isnan.(DSGE.logleveltopct_4q_percapita([1., 2., 3., 4.], ones(7)))) == 4
    @test DSGE.logleveltopct_4q_percapita([1., 2., 3., 4.], ones(7), zeros(4)) == 100. .* (exp.([1., 2.,
                                                                                3., 4.] ./ 100. .+
                                                                                [4., 4., 4., 4.]) .- 1.)
    @test sum(isnan.(DSGE.logleveltopct_4q_percapita([1. 2. 3. 4.; 1. 2. 3. 4.], ones(7)))) == 8
    @test DSGE.logleveltopct_4q_percapita([1. 2. 3. 4.; 1. 2. 3. 4.], ones(7), zeros(4)) == 100. .*
              (exp.([1. 2. 3. 4.; 1. 2. 3. 4.] ./ 100. + 4 * ones(2,4)) .- 1.)
    @test_throws AssertionError DSGE.logleveltopct_4q_percapita([1., 2., 3., 4.], [1., 2., 3.])

    @test sum(isnan.(loggrowthtopct_4q_approx([.1, .2]))) == 2
    @test @test_matrix_approx_eq loggrowthtopct_4q_approx([.1, .2], zeros(3)) [.1, .3]
    @test @test_matrix_approx_eq loggrowthtopct_4q_approx([.1 .2; .1 .2], zeros(3)) [.1 .3; .1 .3]
    @test_throws AssertionError loggrowthtopct_4q_approx([.1, .2], [.1, .2])

    @test sum(isnan.(DSGE.logleveltopct_4q_approx([.1, .2, .3, .4]))) == 4
    @test @test_matrix_approx_eq DSGE.logleveltopct_4q_approx([.1, .2, .3, .4], [.1, .1, .1, .1]) [0., .1, .2, .3]
    @test_throws AssertionError DSGE.logleveltopct_4q_approx([.1, .2], [.1, .2])
end

nothing


text:
using Oceananigans: fields

"""
    RungeKutta3TimeStepper{FT, TG} <: AbstractTimeStepper

Holds parameters and tendency fields for a low storage, third-order Runge-Kutta-Wray
time-stepping scheme described by Le and Moin (1991).
"""
struct RungeKutta3TimeStepper{FT, TG} <: AbstractTimeStepper
    γ¹ :: FT
    γ² :: FT
    γ³ :: FT
    ζ² :: FT
    ζ³ :: FT
    Gⁿ :: TG
    G⁻ :: TG
end

"""
    RungeKutta3TimeStepper(arch, grid, tracers,
                           Gⁿ = TendencyFields(arch, grid, tracers),
                           G⁻ = TendencyFields(arch, grid, tracers))

Return an `RungeKutta3TimeStepper` object with tendency fields on `arch` and
`grid`. The tendency fields can be specified via optional kwargs.
"""
function RungeKutta3TimeStepper(arch, grid, tracers;
                                Gⁿ = TendencyFields(arch, grid, tracers),
                                G⁻ = TendencyFields(arch, grid, tracers))

    γ¹ = 8 // 15
    γ² = 5 // 12
    γ³ = 3 // 4

    ζ² = -17 // 60
    ζ³ = -5 // 12

    return RungeKutta3TimeStepper{eltype(grid), typeof(Gⁿ)}(γ¹, γ², γ³, ζ², ζ³, Gⁿ, G⁻)
end

#####
##### Time steppping
#####

"""
    time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 3rd-order Runge-Kutta method.
The 3rd-order Runge-Kutta method takes three intermediate substep stages to
achieve a single timestep. A pressure correction step is applied at each intermediate
stage.
"""
function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt)
    Δt == 0 && @warn "Δt == 0 may cause model blowup!"

    # Be paranoid and update state at iteration 0, in case run! is not used:
    model.clock.iteration == 0 && update_state!(model)

    γ¹ = model.timestepper.γ¹
    γ² = model.timestepper.γ²
    γ³ = model.timestepper.γ³

    ζ² = model.timestepper.ζ²
    ζ³ = model.timestepper.ζ³

    first_stage_Δt  = γ¹ * Δt
    second_stage_Δt = (γ² + ζ²) * Δt
    third_stage_Δt  = (γ³ + ζ³) * Δt

    #
    # First stage
    #

    calculate_tendencies!(model)

    correct_immersed_tendencies!(model, Δt, γ¹, 0)

    rk3_substep!(model, Δt, γ¹, nothing)

    calculate_pressure_correction!(model, first_stage_Δt)
    pressure_correct_velocities!(model, first_stage_Δt)

    tick!(model.clock, first_stage_Δt; stage=true)
    store_tendencies!(model)
    update_state!(model)
    update_particle_properties!(model, first_stage_Δt)

    #
    # Second stage
    #

    calculate_tendencies!(model)

    correct_immersed_tendencies!(model, Δt, γ², ζ²)

    rk3_substep!(model, Δt, γ², ζ²)

    calculate_pressure_correction!(model, second_stage_Δt)
    pressure_correct_velocities!(model, second_stage_Δt)

    tick!(model.clock, second_stage_Δt; stage=true)
    store_tendencies!(model)
    update_state!(model)
    update_particle_properties!(model, second_stage_Δt)

    #
    # Third stage
    #

    calculate_tendencies!(model)

    correct_immersed_tendencies!(model, Δt, γ³, ζ³)

    rk3_substep!(model, Δt, γ³, ζ³)

    calculate_pressure_correction!(model, third_stage_Δt)
    pressure_correct_velocities!(model, third_stage_Δt)

    tick!(model.clock, third_stage_Δt)
    update_state!(model)
    update_particle_properties!(model, third_stage_Δt)

    return nothing
end

#####
##### Time stepping in each substep
#####

function rk3_substep!(model, Δt, γⁿ, ζⁿ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    substep_field_kernel! = rk3_substep_field!(device(model.architecture), workgroup, worksize)

    model_fields = fields(model)

    events = []

    for (i, field) in enumerate(model_fields)

        field_event = substep_field_kernel!(field, Δt, γⁿ, ζⁿ,
                                            model.timestepper.Gⁿ[i],
                                            model.timestepper.G⁻[i],
                                            dependencies=barrier)

        push!(events, field_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end

"""
Time step fields via the 3rd-order Runge-Kutta method

    `U^{m+1} = U^m + Δt (γⁿ G^{m} + ζⁿ G^{m-1})`,

where `m` denotes the substage.
"""

"""
Time step velocity fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_field!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U[i, j, k] += Δt * (γⁿ * Gⁿ[i, j, k] + ζⁿ * G⁻[i, j, k])
    end
end

"""
Time step velocity fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_field!(U, Δt, γ¹, ::Nothing, G¹, G⁰)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U[i, j, k] += Δt * γ¹ * G¹[i, j, k]
    end
end


text:
# upper
X = [1 0.1 0.4 0.2; 0.1 1 0.3 0.5; 0.4 0.3 1 0.6; 0.2 0.5 0.6 1]
@test ChromNetPaper.upper(X) == [0.1, 0.4, 0.3, 0.2, 0.5, 0.6]

# network_enrichment
T = round(Bool, [1 1 1 0; 1 1 0 1; 1 0 1 0; 0 1 0 1])
@test network_enrichment(X, T) == (2/3)/0.5
M = round(Bool, [1 1 1 0; 1 1 1 0; 1 1 1 0; 0 0 0 0])
@test network_enrichment(X, T, M) == (1/2)/(2/3)

# id2uniprot
@test id2uniprot("ENCSR177HDZ") == "P01100"
@test id2uniprot("ENCSR664POU") == "Q04206"
@test id2uniprot("ENCSR459FTB") == "P17480"

# id2celltype
@test id2celltype("ENCSR177HDZ") == "HepG2"

# id2target
@test id2target("ENCSR177HDZ") == "FOS"

# id2treatments
@test id2treatments("ENCSR177HDZ") == "None"

# id2truth
@test id2truth("ENCSR177HDZ", "ENCSR664POU")
@test !id2truth("ENCSR177HDZ", "ENCSR459FTB")

# truth_matrix
ans = [
    true  true  false;
    true  true  false;
    false false true;
]
@test all(truth_matrix(["ENCSR177HDZ", "ENCSR664POU", "ENCSR459FTB"]) .== ans)

# ishistone
@test ishistone("ENCSR449AYM")
@test !ishistone("ENCSR177HDZ")

# mask_matrix
ans = [
    false false false;
    false false true;
    false true  false;
]
@test all(mask_matrix("within_all", ["ENCSR177HDZ", "ENCSR664POU", "ENCSR459FTB"]) .== ans)

# enrichment_rank
mask = ChromNetPaper.upper(M)
scores = ChromNetPaper.upper(X)[mask]
truth = ChromNetPaper.upper(T)[mask]
x,y = enrichment_rank(truth, scores)
@test maximum(abs(x .- [1,2,3])) < 1e-8
@test maximum(abs(y .- [(1/1)/(2/3), (1/2)/(2/3), (2/3)/(2/3)])) < 1e-8

# network_enrichment_density
x2,y2 = network_enrichment_density(X, T, M)
@test maximum(abs(y2 .- y)) < 1e-8
@test abs(x2[end] - 1.0) < 1e-8

# network_enrichment_rank
x2,y2 = network_enrichment_rank(X, T, M)
@test maximum(abs(y2 .- y)) < 1e-8
@test abs(x2[end] - length(x2)) < 1e-8


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.freqspaceM)),NamedTuple{(:dim,), Tuple{Int64}},typeof(freqspaceM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.ndgridM)),NamedTuple{(:dim,), Tuple{Int64}},typeof(ndgridM),StepRange{Int64, Int64}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,), Tuple{Matrix{Int64}}},typeof(randM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortperm)),NamedTuple{(:dims,), Tuple{Int64}},typeof(sortperm),Matrix{Int64}})
    Base.precompile(Tuple{typeof(blkdiagM),Matrix{Float64},Matrix{Float64},Vararg{Matrix{Float64}, N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64, 3},Vararg{Array{Float64, 3}, N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Matrix{Float64},Vararg{Matrix{Float64}, N} where N})
    Base.precompile(Tuple{typeof(diagM),Vector{Int64}})
    Base.precompile(Tuple{typeof(eyeM),Int64})
    Base.precompile(Tuple{typeof(flipM),Matrix{Int64}})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64})
    Base.precompile(Tuple{typeof(meshgridM),StepRange{Int64, Int64}})
    Base.precompile(Tuple{typeof(meshgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(ndgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(onesM),Type{Int32},Vector{Int64}})
    Base.precompile(Tuple{typeof(randM),Symbol,Int64})
    Base.precompile(Tuple{typeof(randM),Type,Int64,Int64})
    Base.precompile(Tuple{typeof(randM),Type{Int32},Vector{Int64}})
    Base.precompile(Tuple{typeof(repelemM),Matrix{Int64},Int64})
    Base.precompile(Tuple{typeof(repmatM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(repmatM),Vector{Int64},Vector{Int64}})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Float64, 3},Int64,Int64})
    Base.precompile(Tuple{typeof(sortM),Array{Integer, 3},Int64})
    Base.precompile(Tuple{typeof(sortM),Matrix{Int64}})
    Base.precompile(Tuple{typeof(sortM),Vector{Int64}})
    Base.precompile(Tuple{typeof(squeezeM),Array{Float64, 3}})
    Base.precompile(Tuple{typeof(trueM),Vector{Int64}})
    let fbody = try __lookup_kwbody__(which(sortperm, (Matrix{Int64},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Int64,Bool,typeof(sortperm),Matrix{Int64},))
        end
    end
end

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang