text:
import .Tracker
using .Tracker: Tracker,
                TrackedReal,
                TrackedVector,
                TrackedMatrix,
                TrackedArray,
                TrackedVecOrMat,
                @grad,
                track,
                data,
                param
using LinearAlgebra

# Different from Tracker.istracked
_istracked(x::AbstractArray{<:TrackedReal}) = true
_istracked(::TrackedArray) = false

Tracker.dual(x::Bool, p) = x
Base.prevfloat(r::TrackedReal) = track(prevfloat, r)
@grad function prevfloat(r::Real)
    prevfloat(data(r)), Δ -> Δ
end
Base.nextfloat(r::TrackedReal) = track(nextfloat, r)
@grad function nextfloat(r::Real)
    nextfloat(data(r)), Δ -> Δ
end
for i = 0:2, c = Tracker.combinations([:AbstractArray, :TrackedArray, :TrackedReal, :Number], i), f = [:hcat, :vcat]
    if :TrackedReal in c
        cnames = map(_ -> gensym(), c)
        @eval Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::Union{TrackedArray,TrackedReal}, xs::Union{AbstractArray,Number}...) =
            track($f, $(cnames...), x, xs...)
    end
end
@grad function vcat(x::Real)
    vcat(data(x)), (Δ) -> (Δ[1],)
end
@grad function vcat(x1::Real, x2::Real)
    vcat(data(x1), data(x2)), (Δ) -> (Δ[1], Δ[2])
end
@grad function vcat(x1::AbstractVector, x2::Real)
    vcat(data(x1), data(x2)), (Δ) -> (Δ[1:length(x1)], Δ[length(x1)+1])
end

function Base.copy(
    A::TrackedArray{T, 2, <:Adjoint{T, <:AbstractTriangular{T, <:AbstractMatrix{T}}}},
) where {T <: Real}
    return track(copy, A)
end
@grad function Base.copy(
    A::TrackedArray{T, 2, <:Adjoint{T, <:AbstractTriangular{T, <:AbstractMatrix{T}}}},
) where {T <: Real}
    return copy(data(A)), ∇ -> (copy(∇),)
end

Base.:*(A::TrackedMatrix, B::AbstractTriangular) = track(*, A, B)
Base.:*(A::AbstractTriangular{T}, B::TrackedVector) where {T} = track(*, A, B)
Base.:*(A::AbstractTriangular{T}, B::TrackedMatrix) where {T} = track(*, A, B)
Base.:*(A::Adjoint{T, <:AbstractTriangular{T}}, B::TrackedMatrix) where {T} = track(*, A, B)
Base.:*(A::Adjoint{T, <:AbstractTriangular{T}}, B::TrackedVector) where {T} = track(*, A, B)

_eps(::Type{<:TrackedReal{T}}) where {T} = _eps(T)

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:TrackerAD}, Inverse{<:ADBijector{<:TrackerAD}}},
    x::Real
)
    return data(Tracker.gradient(b, x)[1])
end
function jacobian(
    b::Union{<:ADBijector{<:TrackerAD}, Inverse{<:ADBijector{<:TrackerAD}}},
    x::AbstractVector{<:Real}
)
    # We extract `data` so that we don't return a `Tracked` type
    return data(Tracker.jacobian(b, x))
end

# implementations for Shift bijector
function _logabsdetjac_shift(a::TrackedReal, x::Real, ::Val{0})
    return param(_logabsdetjac_shift(data(a), data(x), Val(0)))
end
function _logabsdetjac_shift(a::TrackedReal, x::AbstractVector{<:Real}, ::Val{0})
    return param(_logabsdetjac_shift(data(a), data(x), Val(0)))
end
function _logabsdetjac_shift(
    a::Union{TrackedReal, TrackedVector{<:Real}},
    x::AbstractVector{<:Real},
    ::Val{1}
)
    return param(_logabsdetjac_shift(data(a), data(x), Val(1)))
end
function _logabsdetjac_shift(
    a::Union{TrackedReal, TrackedVector{<:Real}},
    x::AbstractMatrix{<:Real},
    ::Val{1}
)
    return param(_logabsdetjac_shift(data(a), data(x), Val(1)))
end

# Log bijector

@grad function logabsdetjac(b::Log{1}, x::AbstractVector)
    return -sum(log, data(x)), Δ -> (nothing, -Δ ./ data(x))
end
@grad function logabsdetjac(b::Log{1}, x::AbstractMatrix)
    return -vec(sum(log, data(x); dims = 1)), Δ -> (nothing, .- Δ' ./ data(x))
end

# implementations for Scale bijector
# Adjoints for 0-dim and 1-dim `Scale` using `Real`
function _logabsdetjac_scale(a::TrackedReal, x::Real, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::Real, ::Val{0})
    return _logabsdetjac_scale(data(a), data(x), Val(0)), Δ -> (inv(data(a)) .* Δ, nothing, nothing)
end
# Need to treat `AbstractVector` and `AbstractMatrix` separately due to ambiguity errors
function _logabsdetjac_scale(a::TrackedReal, x::AbstractVector, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0})
    da = data(a)
    J = fill(inv.(da), length(x))
    return _logabsdetjac_scale(da, data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedReal, x::AbstractMatrix, ::Val{0})
    return track(_logabsdetjac_scale, a, data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{0})
    da = data(a)
    J = fill(size(x, 1) / da, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
# adjoints for 1-dim and 2-dim `Scale` using `AbstractVector`
function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    return track(_logabsdetjac_scale, a, data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    da = data(a)
    J = inv.(da)
    return _logabsdetjac_scale(da, data(x), Val(1)), Δ -> (J .* Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    return track(_logabsdetjac_scale, a, data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    da = data(a)
    Jᵀ = repeat(inv.(da), 1, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end
# TODO: implement analytical gradient for scaling a vector using a matrix
# function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     track(_logabsdetjac_scale, a, data(x), Val{1})
# end
# @grad function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     throw
# end

# implementations for Stacked bijector
function logabsdetjac(b::Stacked, x::TrackedMatrix{<:Real})
    return mapvcat(eachcol(x)) do c
        logabsdetjac(b, c)
    end
end
# TODO: implement custom adjoint since we can exploit block-diagonal nature of `Stacked`
function (sb::Stacked)(x::TrackedMatrix{<:Real})
    return eachcolmaphcat(sb, x)
end

# Simplex adjoints

function _simplex_bijector(X::TrackedVecOrMat, b::SimplexBijector)
    return track(_simplex_bijector, X, b)
end
function _simplex_inv_bijector(Y::TrackedVecOrMat, b::SimplexBijector)
    return track(_simplex_inv_bijector, Y, b)
end

@grad function _simplex_bijector(X::AbstractVector, b::SimplexBijector)
    Xd = data(X)
    return _simplex_bijector(Xd, b), Δ -> (simplex_link_jacobian(Xd)' * Δ, nothing)
end
@grad function _simplex_inv_bijector(Y::AbstractVector, b::SimplexBijector)
    Yd = data(Y)
    return _simplex_inv_bijector(Yd, b), Δ -> (simplex_invlink_jacobian(Yd)' * Δ, nothing)
end

@grad function _simplex_bijector(X::AbstractMatrix, b::SimplexBijector)
    Xd = data(X)
    return _simplex_bijector(Xd, b), Δ -> begin
        maphcat(eachcol(Xd), eachcol(Δ)) do c1, c2
            simplex_link_jacobian(c1)' * c2
        end, nothing
    end
end
@grad function _simplex_inv_bijector(Y::AbstractMatrix, b::SimplexBijector)
    Yd = data(Y)
    return _simplex_inv_bijector(Yd, b), Δ -> begin
        maphcat(eachcol(Yd), eachcol(Δ)) do c1, c2
            simplex_invlink_jacobian(c1)' * c2
        end, nothing
    end
end

replace_diag(::typeof(log), X::TrackedMatrix) = track(replace_diag, log, X)
@grad function replace_diag(::typeof(log), X)
    Xd = data(X)
    f(i, j) = i == j ? log(Xd[i, j]) : Xd[i, j]
    out = f.(1:size(Xd, 1), (1:size(Xd, 2))')
    out, ∇ -> begin
        g(i, j) = i == j ? ∇[i, j]/Xd[i, j] : ∇[i, j]
        return (nothing, g.(1:size(Xd, 1), (1:size(Xd, 2))'))
    end
end

replace_diag(::typeof(exp), X::TrackedMatrix) = track(replace_diag, exp, X)
@grad function replace_diag(::typeof(exp), X)
    Xd = data(X)
    f(i, j) = ifelse(i == j, exp(Xd[i, j]), Xd[i, j])
    out = f.(1:size(Xd, 1), (1:size(Xd, 2))')
    out, ∇ -> begin
        g(i, j) = ifelse(i == j, ∇[i, j]*exp(Xd[i, j]), ∇[i, j])
        return (nothing, g.(1:size(Xd, 1), (1:size(Xd, 2))'))
    end
end

logabsdetjac(b::SimplexBijector, x::TrackedVecOrMat) = track(logabsdetjac, b, x)
@grad function logabsdetjac(b::SimplexBijector, x::AbstractVector)
    xd = data(x)
    return logabsdetjac(b, xd), Δ -> begin
        (nothing, simplex_logabsdetjac_gradient(xd) * Δ)
    end
end
@grad function logabsdetjac(b::SimplexBijector, x::AbstractMatrix)
    xd = data(x)
    return logabsdetjac(b, xd), Δ -> begin
        (nothing, maphcat(eachcol(xd), Δ) do c, g
            simplex_logabsdetjac_gradient(c) * g
        end)
    end
end

for header in [
    (:(u::TrackedArray), :w),
    (:u, :(w::TrackedArray)),
    (:(u::TrackedArray), :(w::TrackedArray)),
]
    @eval begin
        function get_u_hat($(header...))
            if u isa TrackedArray
                T = typeof(u)
            else
                T = typeof(w)
            end
            x = w' * u
            return (u .+ (planar_flow_m(x) - x) .* w ./ sum(abs2, w))::T
        end
    end
end

for header in [
    (:(z::TrackedArray), :w, :b),
    (:z, :(w::TrackedArray), :b),
    (:z, :w, :(b::TrackedReal)),
    (:(z::TrackedArray), :(w::TrackedArray), :b),
    (:(z::TrackedArray), :w, :(b::TrackedReal)),
    (:z, :(w::TrackedArray), :(b::TrackedReal)),
    (:(z::TrackedArray), :(w::TrackedArray), :(b::TrackedReal)),
]
    @eval begin
        function ψ($(header...))
            if z isa AbstractMatrix
                if z isa TrackedMatrix
                    T = typeof(z)
                elseif w isa TrackedVector
                    T = matrixof(typeof(w))
                else
                    T = matrixof(typeof(b))
                end
            else
                if z isa TrackedVector
                    T = typeof(z)
                elseif w isa TrackedVector
                    T = typeof(w)
                else
                    T = vectorof(typeof(b))
                end
            end
            return ((1 .- tanh.(w' * z .+ b).^2) .* w)::T    # for planar flow from eq(11)
        end
    end
end

for header in [
    (:(u::TrackedArray), :w, :b, :(z::AbstractVecOrMat)),
    (:u, :(w::TrackedArray), :b, :(z::AbstractVecOrMat)),
    (:u, :w, :(b::TrackedReal), :(z::AbstractVecOrMat)),
    (:u, :w, :b, :(z::TrackedVecOrMat)),
    (:(u::TrackedArray), :(w::TrackedArray), :b, :(z::AbstractVecOrMat)),
    (:(u::TrackedArray), :w, :(b::TrackedReal), :(z::AbstractVecOrMat)),
    (:(u::TrackedArray), :w, :b, :(z::TrackedVecOrMat)),
    (:u, :(w::TrackedArray), :(b::TrackedReal), :(z::AbstractVecOrMat)),
    (:u, :(w::TrackedArray), :b, :(z::TrackedVecOrMat)),
    (:u, :w, :(b::TrackedArray), :(z::TrackedVecOrMat)),
    (:(u::TrackedArray), :(w::TrackedArray), :(b::TrackedReal), :(z::AbstractVecOrMat)),
    (:(u::TrackedArray), :(w::TrackedArray), :b, :(z::TrackedVecOrMat)),
    (:(u::TrackedArray), :w, :(b::TrackedReal), :(z::TrackedVecOrMat)),
    (:u, :(w::TrackedArray), :(b::TrackedReal), :(z::TrackedVecOrMat)),
    (:(u::TrackedArray), :(w::TrackedArray), :(b::TrackedReal), :(z::TrackedVecOrMat)),
]
    @eval begin
        function _planar_transform($(header...))
            u_hat = get_u_hat(u, w)
            if z isa AbstractVector
                temp = w' * z + b + zero(eltype(u_hat))
                if z isa TrackedVector
                    T = typeof(z)
                elseif u_hat isa TrackedVector
                    T = typeof(u_hat)
                else
                    T = vectorof(typeof(temp))
                end
            else
                temp = w' * z .+ (b + zero(eltype(u_hat)))
                if z isa TrackedMatrix
                    T = typeof(z)
                elseif u_hat isa TrackedVector
                    T = matrixof(typeof(u_hat))
                else
                    T = matrixof(typeof(temp'))
                end
            end
            transformed::T = z .+ u_hat .* tanh.(temp) # from eq(10)
            return (transformed = transformed, u_hat = u_hat)
        end
    end
end

for header in [
    (:(α_::TrackedReal), :β, :z_0, :(z::AbstractVector)),
    (:α_, :(β::TrackedReal), :z_0, :(z::AbstractVector)),
    (:α_, :β, :(z_0::TrackedVector), :(z::AbstractVector)),
    (:α_, :β, :z_0, :(z::TrackedVector)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::AbstractVector)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::AbstractVector)),
    (:(α_::TrackedReal), :β, :z_0, :(z::TrackedVecOrMat)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::AbstractVector)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::TrackedVector)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::TrackedVector)),
    (:α_, :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedVector)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedVector)),
]
    @eval begin
        function _radial_transform($(header...))
            α = softplus(α_)            # from A.2
            β_hat = -α + softplus(β)    # from A.2
            if β_hat isa TrackedReal
                TV = vectorof(typeof(β_hat))
                T = vectorof(typeof(β_hat))
            elseif z_0 isa TrackedVector
                TV = typeof(z_0)
                T = typeof(z_0)
            else
                T = TV = typeof(z)
            end
            Tr = promote_type(eltype(z), eltype(z_0))
            r::Tr = norm((z .- z_0)::TV)
            transformed::T = z .+ β_hat ./ (α .+ r') .* (z .- z_0)   # from eq(14)
            return (transformed = transformed, α = α, β_hat = β_hat, r = r)
        end
    end
end

for header in [
    (:(α_::TrackedReal), :β, :z_0, :(z::AbstractMatrix)),
    (:α_, :(β::TrackedReal), :z_0, :(z::AbstractMatrix)),
    (:α_, :β, :(z_0::TrackedVector), :(z::AbstractMatrix)),
    (:α_, :β, :z_0, :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::AbstractMatrix)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::AbstractMatrix)),
    (:(α_::TrackedReal), :β, :z_0, :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::AbstractMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :z_0, :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :β, :(z_0::TrackedVector), :(z::TrackedMatrix)),
    (:α_, :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedMatrix)),
    (:(α_::TrackedReal), :(β::TrackedReal), :(z_0::TrackedVector), :(z::TrackedMatrix)),
]
    @eval begin
        function _radial_transform($(header...))
            α = softplus(α_)            # from A.2
            β_hat = -α + softplus(β)    # from A.2
            if β_hat isa TrackedReal
                TV = vectorof(typeof(β_hat))
                T = matrixof(TV)
            elseif z_0 isa TrackedVector
                TV = typeof(z_0)
                T = matrixof(TV)
            else
                T = typeof(z)
                TV = vectorof(T)
            end
            r::TV = mapvcat(norm, eachcol(z .- z_0))
            transformed::T = z .+ β_hat ./ (α .+ r') .* (z .- z_0)   # from eq(14)
            return (transformed = transformed, α = α, β_hat = β_hat, r = r)
        end
    end
end

function matrixof(::Type{TrackedArray{T, 1, Vector{T}}}) where {T <: Real}
    return TrackedArray{T, 2, Matrix{T}}
end
function matrixof(::Type{TrackedReal{T}}) where {T <: Real}
    return TrackedArray{T, 2, Matrix{T}}
end
function vectorof(::Type{TrackedArray{T, 2, Matrix{T}}}) where {T <: Real}
    return TrackedArray{T, 1, Vector{T}}
end
function vectorof(::Type{TrackedReal{T}}) where {T <: Real}
    return TrackedArray{T, 1, Vector{T}}
end

(b::Exp)(x::TrackedVector) = exp.(x)::vectorof(float(eltype(x)))
(b::Log)(x::TrackedVector) = log.(x)::vectorof(float(eltype(x)))
logabsdetjac(b::Log{0}, x::TrackedVector) = .-log.(x)::vectorof(float(eltype(x)))
logabsdetjac(b::Log{1}, x::TrackedMatrix) = - vec(sum(log.(x); dims = 1))

_logabsdetjac_shift(a::Real, x::TrackedVector{T}, ::Val{0}) where {T<:Real} = zeros(T, length(x))
_logabsdetjac_shift(a::T1, x::TrackedVector{T2}, ::Val{1}) where {T1<:Union{Real, TrackedVector}, T2<:Real} = zero(T2)
_logabsdetjac_shift(a::T1, x::TrackedMatrix{T2}, ::Val{1}) where {T1<:Union{Real, TrackedVector}, T2<:Real} = zeros(T2, size(x, 2))


text:
using LaserTag
using POMDPGifs
using QMDP
using Random
using ParticleFilters

rng = MersenneTwister(7)

m = gen_lasertag(rng=rng, robot_position_known=true)
policy = solve(QMDPSolver(verbose=true), m)
filter = SIRParticleFilter(m, 10000, rng=rng)

@show makegif(m, policy, filter, filename="out.gif", rng=rng)

ext: jl


text:
using Test

@testset "Mutts.jl" begin
    include("Mutts.jl")
end


text:
module TestCommon
# tests that can be run by any physics module

using Base.Test

using ODLCommonTools
using PDESolver
using Utils
using OptimizationInterface
using PdePumiInterface

# fd_tests.jl
export testMeshDerivatives, testResidualDerivatives, testDJDx

include("fd_tests.jl")

end  # end module


text:
using YaoBase
export Concentrator, concentrate

"""
    Concentrator{N, T, BT <: AbstractBlock} <: AbstractContainer{BT, N, T}

concentrates serveral lines together in the circuit, and expose
it to other blocks.
"""
struct Concentrator{N,BT<:AbstractBlock,C} <: AbstractContainer{BT,N}
    content::BT
    locs::NTuple{C,Int}
end

function Concentrator{N}(block::BT, locs::NTuple{C,Int}) where {N,M,C,BT<:AbstractBlock{M}}
    if !(length(locs) == M && N >= M)
        throw(LocationConflictError("length of locs must be equal to the size of block, and smaller than size of itself."))
    end
    return Concentrator{N,BT,C}(block, locs)
end

"""
    concentrate(n, block, locs)

Create a [`Concentrator`](@ref) block with total number of current active qubits `n`,
which concentrates given wire location together to `length(locs)` active qubits,
and relax the concentration afterwards.

# Example

Concentrator is equivalent to [`put`](@ref) a block on given position mathematically, but more efficient
and convenient for large blocks.

```jldoctest; setup=:(using YaoBlocks; using YaoArrayRegister)
julia> r = rand_state(3)
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 3/3

julia> apply!(copy(r), concentrate(X, 1)) ≈ apply!(copy(r), put(1=>X))
true
```

It works for in-contigious locs as well

```jldoctest; setup=:(using YaoBlocks; using YaoArrayRegister)
julia> r = rand_state(4)
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 4/4

julia> cc = concentrate(4, kron(X, Y), (1, 3))
nqubits: 4
Concentrator: (1, 3)
└─ kron
   ├─ 1=>X gate
   └─ 2=>Y gate

julia> pp = chain(4, put(1=>X), put(3=>Y))
nqubits: 4
chain
├─ put on (1)
│  └─ X gate
└─ put on (3)
   └─ Y gate

julia> apply!(copy(r), cc) ≈ apply!(copy(r), pp)
true
```
"""
function concentrate(n::Int, block::AbstractBlock, locs)
    return Concentrator{n}(block, Tuple(locs))
end

# support lazy qubits
concentrate(n::Int, block::Function, locs) = concentrate(n, parse_block(length(locs), block), locs)

"""
    concentrate(block, locs) -> f(n)

Lazy curried version of [`concentrate`](@ref).
"""
concentrate(block::AbstractBlock, locs) = @λ(n -> concentrate(n, block, locs))
concentrate(block::Function, locs) = @λ(n -> concentrate(n, block, locs))

occupied_locs(c::Concentrator) = map(i -> c.locs[i], c.content |> occupied_locs)
chsubblocks(pb::Concentrator{N}, blk::AbstractBlock) where {N} = Concentrator{N}(blk, pb.locs)
PreserveTrait(::Concentrator) = PreserveAll()

function apply!(r::AbstractRegister, c::Concentrator)
    _check_size(r, c)
    focus!(r, c.locs)
    apply!(r, c.content)
    relax!(r, c.locs, to_nactive = nqubits(c))
    return r
end

function mat(::Type{T}, c::Concentrator{N,<:AbstractBlock}) where {N,T}
    mat(T, PutBlock{N}(c.content, c.locs))
end

Base.adjoint(blk::Concentrator{N}) where {N} = Concentrator{N}(adjoint(blk.content), blk.locs)

function Base.:(==)(a::Concentrator{N,BT}, b::Concentrator{N,BT}) where {N,BT}
    return a.content == b.content && a.locs == b.locs
end

YaoBase.nqubits(::Concentrator{N}) where {N} = N
YaoBase.nactive(c::Concentrator) = length(c.locs)

function YaoBase.iscommute(x::Concentrator{N}, y::Concentrator{N}) where {N}
    isempty(setdiff(occupied_locs(x), occupied_locs(y))) && return true
    if x.locs == y.locs
        return iscommute(x.content, y.content)
    else
        return iscommute_fallback(x, y)
    end
end


text:
function reduce_index(a, j)
    return a < j ? a : a - 1
end

function recover_index(a, j)
    return a < j ? a : a + 1
end

function reduce_omega(omega::Array{<:Real}, i::Int64, j::Int64)::Array{Float64}
    omega_tilde = copy(omega)
    omega_tilde[i] += omega_tilde[j]
    omega_tilde = omega_tilde[[1:(j - 1); (j + 1):end]]
    return copy(omega_tilde)
end

function reduce_lambda(lambda::Array{<:Real}, j::Int64)::Array{Float64}
    new_lambda = copy(lambda)[[1:(j - 1); (j + 1):end]]
    return copy(new_lambda)
end

function recover_lambda(lambda_tilde::Array{<:Real}, i::Int64, j::Int64, epsilon::Real)::Array{Float64}
    lambda = [lambda_tilde[1:(j - 1)]; lambda_tilde[reduce_index(i, j)] + epsilon ; lambda_tilde[j:end]]
    return copy(lambda)
end

function reduce_mu(mu::Array{<:Real}, omega::Array{<:Real}, i::Int64, j::Int64, epsilon::Real)::Array{Float64}
    mu_tilde = mu[[1:(j - 1); (j + 1):end]]
    mu_tilde[reduce_index(i, j)] = (omega[i] * mu[i] + omega[j] * mu[j]) / (omega[i] + omega[j]) - (epsilon * omega[i]) / (omega[i] + omega[j])
    return copy(mu_tilde)
end

function reduce_laplacian(Laplacian::Array{<:Real}, i::Int64, j::Int64)::Array{Float64}
    Laplacian_tilde = copy(Laplacian)
    Laplacian_tilde[i,:] += Laplacian_tilde[j,:]
    Laplacian_tilde[:,i] += Laplacian_tilde[:,j]
    Laplacian_tilde = Laplacian_tilde[[1:(j - 1); (j + 1):end], [1:(j - 1); (j + 1):end]]
    # Laplacian_tilde[reduce_index(i, j), reduce_index(i, j)] -= sum(Laplacian_tilde[reduce_index(i, j), :])
    return copy(Laplacian_tilde)
end

function reduce_laplacian_j(Laplacian::Array{<:Real}, i::Int64, j::Int64)::Array{Float64}
    Laplacian_tilde_j = Laplacian[:, j]
    Laplacian_tilde_j[j] += Laplacian_tilde_j[i]
    Laplacian_tilde_j = Laplacian_tilde_j[[(1:i - 1); (i + 1):end]]
    return copy(Laplacian_tilde_j)
end

max_stars_repo_name: tomaskocak/BAI

max_issues_repo_name: tomaskocak/BAI

max_forks_repo_name: tomaskocak/BAI


text:
const source_dir_name = "source";
const desc_file_name = "description.myresults.txt";
const log_file_name = "log.myresults.txt";
const data_file_name = "data.myresults.jld2";

struct ResDirTree
    dirname::AbstractString;
    basename::AbstractString;
    ResDirTree(dir_name, base_name) =
        new(dir_name, basename(base_name));
    ResDirTree(res_dir) =
        new(dirname(res_dir), basename(res_dir));
end

Base.basename(tree::ResDirTree) = tree.basename;
Base.dirname(tree::ResDirTree) = tree.dirname;
resdir(tree::ResDirTree) = joinpath(dirname(tree), basename(tree));
Base.isdir(tree::ResDirTree) = isdir(resdir(tree))
sourcedir(tree::ResDirTree) = joinpath(resdir(tree), source_dir_name);
descfile(tree::ResDirTree) = joinpath(resdir(tree), desc_file_name);
logfile(tree::ResDirTree) = joinpath(resdir(tree), log_file_name);
datafile(tree::ResDirTree) = joinpath(resdir(tree), data_file_name);

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
#===== stuff =====#

t2 = time()

@testset "parsing all the things" begin include("parsing.jl") end

@testset "tests from Einsum.jl" begin include("einsum.jl") end

@info @sprintf("Basic tests took %.1f seconds", time()-t2)

@testset "internal pieces" begin include("utils.jl") end

@testset "matrix multiplication" begin
    # size 200 is big enough to test block_halves even with MINIBLOCK = 64^3
    @testset "size $N, elements $T" for N in [2, 20, 200], T in [1:99, Float32, Float64, ComplexF64]
        for f in [identity, adjoint]
            A = f(rand(T, N,N));
            B = f(rand(T, N,N));
            @test A * B ≈ @tullio C[i,k] := A[i,j] * B[j,k]
        end
        if N < 200
            X = rand(T, N,N+1);
            Y = rand(T, N+1,N+2);
            Z = rand(T, N+2,N+1);
            @test X * Y * Z ≈ @tullio C[a,d] := X[a,b] * Y[b,c] * Z[c,d]
        end
    end
    @testset "@allocated" begin
        m!(C,A,B) = @tullio C[i,k] = A[i,j] * B[j,k] threads=false
        C1, A1, B1 = rand(4,4), rand(4,4), rand(4,4)
        @allocated m!(C1, A1, B1)
        @test 0 == @allocated m!(C1, A1, B1)
    end
end

#===== Tracker =====#

t3 = time()
using Tracker

GRAD = :Tracker
_gradient(x...) = Tracker.gradient(x...)

@tullio grad=Base
@testset "gradients: Tracker + DiffRules" begin include("gradients.jl") end

@tullio grad=Dual
@testset "gradients: Tracker + ForwardDiff" begin include("gradients.jl") end

@info @sprintf("Tracker tests took %.1f seconds", time()-t3)


text:
#=
A lot of visualization forms in GLVisualize are realised in the form of instanced
particles. This is because they can be handled very efficiently by OpenGL.
There are quite a few different ways to feed instances with different attributes.
The main constructor for particles is a tuple of (Primitive, Position), whereas
position can come in all forms and shapes. You can leave away the primitive.
In that case, GLVisualize will fill in some default that is anticipated to make
the most sense for the datatype.
=#

#3D primitives
const Primitives3D = Union{AbstractGeometry{3}, AbstractMesh}
#2D primitives AKA sprites, since they are shapes mapped onto a 2D rectangle
const Sprites = Union{AbstractGeometry{2}, Shape, Char, Type}
const AllPrimitives = Union{AbstractGeometry, Shape, Char, AbstractMesh}

using AbstractPlotting: RectanglePacker

# There is currently no way to get the two following two signatures
# under one function, which is why we delegate to meshparticle
function _default(
        p::Tuple{TOrSignal{Pr}, VectorTypes{P}}, s::Style, data::Dict
    ) where {Pr <: Primitives3D, P <: Point}
    return meshparticle(p, s, data)
end

function to_meshcolor(color::TOrSignal{Vector{T}}) where T <: Colorant
    TextureBuffer(color)
end

function to_meshcolor(color::TOrSignal{Matrix{T}}) where T <: Colorant
    Texture(color)
end
function to_meshcolor(color)
    color
end

function to_mesh(mesh::TOrSignal{<: GeometryPrimitive})
    return NativeMesh(const_lift(GeometryBasics.normal_mesh, mesh))
end

function to_mesh(mesh::TOrSignal{<: GeometryBasics.Mesh})
    return NativeMesh(mesh)
end

using AbstractPlotting
using AbstractPlotting: get_texture_atlas

vec2quaternion(rotation::StaticVector{4}) = rotation

function vec2quaternion(r::StaticVector{2})
    vec2quaternion(Vec3f0(r[1], r[2], 0))
end
function vec2quaternion(rotation::StaticVector{3})
    AbstractPlotting.rotation_between(Vec3f0(0, 0, 1), Vec3f0(rotation))
end

vec2quaternion(rotation::Vec4f0) = rotation
vec2quaternion(rotation::VectorTypes) = const_lift(x-> vec2quaternion.(x), rotation)
vec2quaternion(rotation::Node) = lift(vec2quaternion, rotation)
vec2quaternion(rotation::AbstractPlotting.Quaternion)= Vec4f0(rotation.data)
GLAbstraction.gl_convert(rotation::AbstractPlotting.Quaternion)= Vec4f0(rotation.data)

"""
This is the main function to assemble particles with a GLNormalMesh as a primitive
"""
function meshparticle(p, s, data)
    rot = get!(data, :rotation, Vec4f0(0, 0, 0, 1))
    rot = vec2quaternion(rot)
    delete!(data, :rotation)
    @gen_defaults! data begin
        primitive = p[1] => to_mesh
        position = p[2] => TextureBuffer
        position_x = nothing => TextureBuffer
        position_y = nothing => TextureBuffer
        position_z = nothing => TextureBuffer

        scale = Vec3f0(1) => TextureBuffer
        scale_x = nothing => TextureBuffer
        scale_y = nothing => TextureBuffer
        scale_z = nothing => TextureBuffer

        rotation = rot => TextureBuffer
        texturecoordinates = nothing
        shading = true
    end

    @gen_defaults! data begin
        color_map = nothing => Texture
        color_norm = nothing
        intensity = nothing
        image = nothing
        color = if color_map == nothing
            default(RGBA{Float32}, s)
        else
            nothing
        end => to_meshcolor
        vertex_color = Vec4f0(1)
        fetch_pixel = false
        uv_scale = Vec2f0(1)

        instances = const_lift(length, position)
        shading = true
        shader = GLVisualizeShader(
            "util.vert", "particles.vert", "fragment_output.frag", "standard.frag",
            view = Dict(
                "position_calc" => position_calc(position, position_x, position_y, position_z, TextureBuffer),
                "light_calc" => light_calc(shading)
            )
        )
    end
    if AbstractPlotting.to_value(intensity) != nothing
        if AbstractPlotting.to_value(position) != nothing
            data[:intensity] = intensity_convert_tex(intensity, position)
            data[:len] = const_lift(length, position)
        else
            data[:intensity] = intensity_convert_tex(intensity, position_x)
            data[:len] = const_lift(length, position_x)
        end
    end
    data
end

to_pointsize(x::Number) = Float32(x)
to_pointsize(x) = Float32(x[1])

struct PointSizeRender
    size::Observable
end
(x::PointSizeRender)() = glPointSize(to_pointsize(x.size[]))
"""
This is the most primitive particle system, which uses simple points as primitives.
This is supposed to be the fastest way of displaying particles!
"""
function _default(position::VectorTypes{T}, s::style"speed", data::Dict) where T <: Point
    @gen_defaults! data begin
        vertex       = position => GLBuffer
        color_map    = nothing  => Texture
        color        = (color_map == nothing ? default(RGBA{Float32}, s) : nothing) => GLBuffer
        color_norm   = nothing
        scale        = 2f0
        shader       = GLVisualizeShader("fragment_output.frag", "dots.vert", "dots.frag")
        gl_primitive = GL_POINTS
    end
    data[:prerender] = PointSizeRender(data[:scale])
    data
end

"""
returns the Shape for the distancefield algorithm
"""
primitive_shape(::Char) = DISTANCEFIELD
primitive_shape(x::X) where {X} = primitive_shape(X)
primitive_shape(::Type{T}) where {T <: Circle} = CIRCLE
primitive_shape(::Type{T}) where {T <: Rect2D} = RECTANGLE
primitive_shape(x::Shape) = x

"""
Extracts the scale from a primitive.
"""
primitive_scale(prim::GeometryPrimitive) = Vec2f0(widths(prim))
primitive_scale(::Union{Shape, Char}) = Vec2f0(40)
primitive_scale(c) = Vec2f0(0.1)

"""
Extracts the offset from a primitive.
"""
primitive_offset(x, scale::Nothing) = Vec2f0(0) # default offset
primitive_offset(x, scale) = const_lift(/, scale, -2f0)  # default offset

"""
Extracts the uv offset and width from a primitive.
"""
primitive_uv_offset_width(c::Char) = glyph_uv_width!(c)
primitive_uv_offset_width(x) = Vec4f0(0,0,1,1)

"""
Gets the texture atlas if primitive is a char.
"""
primitive_distancefield(x) = nothing
primitive_distancefield(::Char) = get_texture!(get_texture_atlas())
primitive_distancefield(::Node{Char}) = get_texture!(get_texture_atlas())

function _default(
        p::Tuple{TOrSignal{Matrix{C}}, VectorTypes{P}}, s::Style, data::Dict
    ) where {C <: Colorant, P <: Point}
    data[:image] = p[1] # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = lift(x-> Vec2f0(size(x)), p[1])
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end
function _default(
        p::Tuple{TOrSignal{Matrix{C}}, VectorTypes{P}}, s::Style, data::Dict
    ) where {C <: AbstractFloat, P <: Point}
    data[:distancefield] = p[1] # we don't want this to be overwritten by user
    @gen_defaults! data begin
        scale = lift(x-> Vec2f0(size(x)), p[1])
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end

function _default(
        p::Tuple{VectorTypes{Matrix{C}}, VectorTypes{P}}, s::Style, data::Dict
    ) where {C <: Colorant, P <: Point}
    images = to_value(p[1])
    isempty(images) && error("Can not display empty vector of images as primitive")
    sizes = map(size, images)
    if !all(x-> x == sizes[1], sizes) # if differently sized
        # create texture atlas
        maxdims = sum(map(Vec{2, Int}, sizes))
        rectangles = map(x->Rect2D(0, 0, x...), sizes)
        rpack = RectanglePacker(Rect2D(0, 0, maxdims...))
        uv_coordinates = [push!(rpack, rect).area for rect in rectangles]
        max_xy = maximum(maximum.(uv_coordinates))
        texture_atlas = Texture(C, (max_xy...,))
        for (area, img) in zip(uv_coordinates, images)
            texture_atlas[area] = img #transfer to texture atlas
        end
        data[:uv_offset_width] = map(uv_coordinates) do uv
            m = max_xy .- 1
            mini = reverse((minimum(uv)) ./ m)
            maxi = reverse((maximum(uv) .- 1) ./ m)
            return Vec4f0(mini..., maxi...)
        end
        images = texture_atlas
    end
    data[:image] = images # we don't want this to be overwritten by user
    @gen_defaults! data begin
        shape = RECTANGLE
        offset = Vec2f0(0)
    end
    sprites(p, s, data)
end

# There is currently no way to get the two following two signatures
# under one function, which is why we delegate to sprites
_default(p::Tuple{TOrSignal{Pr}, VectorTypes{P}}, s::Style, data::Dict) where {Pr <: Sprites, P<:Point} =
    sprites(p,s,data)

function _default(
            p::Tuple{TOrSignal{Pr}, G}, s::Style, data::Dict
        ) where {Pr <: Sprites, G <: Tuple}
    @gen_defaults! data begin
        shape      = const_lift(primitive_shape, p[1])
        position   = nothing => GLBuffer
        position_x = p[2][1] => GLBuffer
        position_y = p[2][2] => GLBuffer
        position_z = length(p[2]) > 2 ? p[2][3] : 0f0 => GLBuffer
    end
    sprites(p, s, data)
end

"""
Main assemble functions for sprite particles.
Sprites are anything like distance fields, images and simple geometries
"""
function sprites(p, s, data)
    rot = get!(data, :rotation, Vec4f0(0, 0, 0, 1))
    rot = vec2quaternion(rot)
    delete!(data, :rotation)

    @gen_defaults! data begin
        shape       = const_lift(x-> Int32(primitive_shape(x)), p[1])
        position    = p[2]    => GLBuffer
        position_x  = nothing => GLBuffer
        position_y  = nothing => GLBuffer
        position_z  = nothing => GLBuffer

        scale       = const_lift(primitive_scale, p[1]) => GLBuffer
        scale_x     = nothing => GLBuffer
        scale_y     = nothing => GLBuffer
        scale_z     = nothing => GLBuffer

        rotation    = rot => GLBuffer
        image       = nothing => Texture
    end
    # TODO don't make this dependant on some shady type dispatch
    if isa(to_value(p[1]), Char) && !isa(to_value(scale), Union{StaticVector, AbstractVector{<: StaticVector}}) # correct dimensions
        data[:scale] = const_lift(correct_scale, p[1], scale)
    end

    @gen_defaults! data begin
        offset          = primitive_offset(p[1], scale) => GLBuffer
        intensity       = nothing => GLBuffer
        color_map       = nothing => Texture
        color_norm      = nothing
        color           = (color_map == nothing ? default(RGBA, s) : nothing) => GLBuffer

        glow_color      = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_color    = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_width    = 0f0
        glow_width      = 0f0
        uv_offset_width = const_lift(primitive_uv_offset_width, p[1]) => GLBuffer

        distancefield   = primitive_distancefield(p[1]) => Texture
        indices         = const_lift(length, p[2]) => to_index_buffer
        # rotation and billboard don't go along
        billboard        = rotation == Vec4f0(0,0,0,1) => "if `billboard` == true, particles will always face camera"
        preferred_camera = :orthographic_pixel
        fxaa             = false
        shader           = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "sprites.geom",
            "sprites.vert", "distance_shape.frag",
            view = Dict("position_calc"=>position_calc(position, position_x, position_y, position_z, GLBuffer))
        )
        gl_primitive = GL_POINTS
    end
    # Exception for intensity, to make it possible to handle intensity with a
    # different length compared to position. Intensities will be interpolated in that case
    if position != nothing
        data[:intensity] = intensity_convert(intensity, position)
        data[:len] = const_lift(length, position)
    else
        data[:intensity] = intensity_convert(intensity, position_x)
        data[:len] = const_lift(length, position_x)
    end
    return data
end

"""
Transforms text into a particle system of sprites, by inferring the
texture coordinates in the texture atlas, widths and positions of the characters.
"""
function _default(main::Tuple{TOrSignal{S}, P}, s::Style, data::Dict) where {S <: AbstractString, P}
    data[:position] = main[2]
    _default(main[1], s, data)
end

function _default(main::TOrSignal{S}, s::Style, data::Dict) where S <: AbstractString
    @gen_defaults! data begin
        relative_scale  = 4 #
        start_position  = Point2f0(0)
        atlas           = get_texture_atlas()
        distancefield   = get_texture!(atlas)
        stroke_width    = 0f0
        glow_width      = 0f0
        font            = to_font("default")
        scale_primitive = true
        position        = const_lift(calc_position, main, start_position, relative_scale, font, atlas)
        offset          = const_lift(calc_offset, main, relative_scale, font, atlas)
        prerender       = () -> begin
            glDisable(GL_DEPTH_TEST)
            glDepthMask(GL_TRUE)
            glDisable(GL_CULL_FACE)
            enabletransparency()
        end
        uv_offset_width = const_lift(main) do str
            Vec4f0[glyph_uv_width!(atlas, c, font) for c = str]
        end
        scale = const_lift(main, relative_scale) do str, s
            Vec2f0[glyph_scale!(atlas, c, font, s) for c = str]
        end
    end
    delete!(data, :font)
    _default((DISTANCEFIELD, position), s, data)
end


text:
"""
Calculates the reference variable for ``μ`` when the elements in vector ``x`` comes from a normal distribution.
See also: [`sample_mean`](@ref), [`sample_variance`](@ref).
"""
function normal_reference_mu_unknown_variance(μ::AbstractFloat, x::Vector)
    n = length(x)
    x̅ = sample_mean(x)
    s = √(sample_variance(x))
    T = (x̅ - μ) / (s / √(n))
    return T
end

"""
Calculates the reference variable for ``μ`` when the elements in vector ``x`` comes from a normal distribution.
This time the variance ``σ²`` has to be known for the distribution.
See also: [`sample_mean`](@ref).
"""
function normal_reference_mu_known_variance(μ::AbstractFloat, σ²::AbstractFloat, x::Vector)
    n = length(x)
    x̅ = sample_mean(x)
    σ = √(σ²)
    T = (x̅ - μ) / (σ / √(n))
    return T
end

"""
Reference variable for ``σ²`` when the elements in vector ``x`` comes from a normal distribution.
See also: [`sample_variance`](@ref).
"""
function normal_reference_sigma(σ²::AbstractFloat, x::Vector)
    n = length(x)
    s² = sample_variance(x)
    T = (n - 1) * s² / σ²
    return T
end

"""
Calculates the reference variable for tests of the difference between two means.
This variable is used when the variance is known and different for the two samples.
See also: [`sample_mean`](@ref).
"""
function normal_reference_two_sample_mu_known_variance(μᵪ::AbstractFloat, μᵧ::AbstractFloat, σ²ᵪ::AbstractFloat, σ²ᵧ::AbstractFloat, x::Vector, y::Vector)
    nᵪ = length(x)
    nᵧ = length(y)
    x̅ = sample_mean(x)
    y̅ = sample_mean(x)
    T = ((x̅ - y̅) - (μᵪ - μᵧ)) / √(σ²ᵪ / nᵪ + σ²ᵧ / nᵧ)
    return T
end

"""
Calculates the reference variable for tests of the difference between two means.
This variable is used when the variance is unknown but equal for the two samples.
See also: [`sample_mean`](@ref), [`two_sample_variance`](@ref).
"""
function normal_reference_two_sample_mu_unknown_equal_variance(μᵪ::AbstractFloat, μᵧ::AbstractFloat, x::Vector, y::Vector)
    nᵪ = length(x)
    nᵧ = length(y)
    x̅ = sample_mean(x)
    y̅ = sample_mean(x)
    s = √(two_sample_variance(x, y))
    T = ((x̅ - y̅) - (μᵪ - μᵧ)) / (s * √(1 / nᵪ + 1 / nᵧ))
    return T
end

"""
Calculates the reference variable for tests of the difference between two means.
This variable is used when the variance is unknown but not necessarily equal for the two samples.
See also: [`sample_mean`](@ref), [`sample_variance`](@ref), [`two_sample_degrees_of_freedom_unknown_variance`](@ref).
"""
function normal_reference_two_sample_mu_unknown_variance(μᵪ::AbstractFloat, μᵧ::AbstractFloat, x::Vector, y::Vector)
    nᵪ = length(x)
    nᵧ = length(y)
    x̅ = sample_mean(x)
    y̅ = sample_mean(x)
    s²ᵪ = sample_variance(x)
    s²ᵧ = sample_variance(y)
    T = ((x̅ - y̅) - (μᵪ - μᵧ)) / √(s²ᵪ / nᵪ + s²ᵧ / nᵧ)
    return T
end

"""
Calculates the degrees of freedom of the t-distribution that the reference variable above is from.
See also: [`sample_variance`](@ref), [`normal_reference_two_sample_mu_unknown_variance`](@ref).
"""
function two_sample_degrees_of_freedom_unknown_variance(x::Vector, y::Vector)
    nᵪ = length(x)
    nᵧ = length(y)
    s²ᵪ = sample_variance(x)
    s²ᵧ = sample_variance(y)
    first = 1 / (nᵪ - 1)
    second = (nᵧ * s²ᵪ)^2 / (nᵧ * s²ᵪ + nᵪ * s²ᵧ)^2
    third = 1 / (nᵧ - 1)
    fourth = (nᵪ * s²ᵧ)^2 / (nᵧ * s²ᵪ + nᵪ * s²ᵧ)^2
    f = 1 / (first * second + third * fourth)
    return f
end

"""
Calculates the interval for the estimated mean difference of two paired samples.
See also: [`sample_mean`](@ref), [`sample_variance`](@ref).
"""
function normal_interval_two_sample_paired(x::Vector, y::Vector, α::AbstractFloat)
    n = length(x)
    z = y .- x
    z̅ = sample_mean(z)
    dist = tDistribution(n-1)
    t = quantile_finder(dist, α)
    s = √(sample_variance(z))
    I = [z̅ - t * s / √(n), z̅ + t * s / √(n)]
    return I
end

function normal_interval_mu_unknown_variance(μ::AbstractFloat, x::Vector, α::AbstractFloat)

end

function normal_interval_mu_known_variance(μ::AbstractFloat, σ²::AbstractFloat, x::Vector, α::AbstractFloat)

end

function normal_interval_sigma(x::Vector, α::AbstractFloat)

end

function normal_interval_two_sample_mu_known_variance(μᵪ::AbstractFloat, μᵧ::AbstractFloat, σ²ᵪ::AbstractFloat, σ²ᵧ::AbstractFloat, x::Vector, y::Vector, α::AbstractFloat)

end

function normal_interval_two_sample_mu_unknown_equal_variance(μᵪ::AbstractFloat, μᵧ::AbstractFloat, x::Vector, y::Vector, α::AbstractFloat)

end

function exp_reference_mu(μ::AbstractFloat, x::Vector)
    n = length(n)
    x̅ = sample_mean(x)
    T = (x̅ - μ) / (μ / √(n))
    return T
end

function exp_interval_mu(μ::AbstractFloat, x::Vector, α::AbstractFloat)

end

function binomial_reference_p()

end

function binomial_interval_p()

end

function hypergeometric_reference_p()

end

function hypergeometric_interval_p()

end

function binomial_interval_comparison_of_proportions()

end

function poisson_reference_mu()

end

function poisson_interval_mu()

end

function test()

end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
type SoftMaxLayer <: Layer
    inputActivation::Vol
    outputActivation::Vol

    function SoftMaxLayer(numOutputs::Int)
        dimensions = (1, 1, numOutputs)
        new(Vol(dimensions, 0.0), Vol(dimensions, 0.0))
    end
end

function forward(self::SoftMaxLayer, inputActivation::Vol)
    @assert inputSize(self) == outputSize(self)
    self.inputActivation = inputActivation
    exponentials = exp(self.inputActivation.w .- maximum(self.inputActivation.w))
    exponentials /= sum(exponentials)

    self.outputActivation.w = exponentials
    @assert abs(sum(self.outputActivation.w) - 1.0) < 1E-5
    return self.outputActivation
end

function backward(self::SoftMaxLayer, class::Int)
    @assert inputSize(self) == outputSize(self)
    fill!(self.inputActivation.dw, 0.0)
    for i in 1:length(self.outputActivation)
        indicator = ((i == class) ? 1.0 : 0.0)
        mul = -(indicator - self.outputActivation.w[1, 1, i])
        self.inputActivation.dw[1, 1, i] = mul
    end
    @assert self.inputActivation.dw |> sum |> abs < 1E-5
    return -log(self.outputActivation.w[1, 1, class])
end

gradientStates(self::SoftMaxLayer) = []


text:
"""
    ScaleTransform(l::Real)

Multiply every element of the input by `l`
```
    l = 2.0
    tr = ScaleTransform(l)
```
"""
struct ScaleTransform{T<:Real} <: Transform
    s::Vector{T}
end

function ScaleTransform(s::T=1.0) where {T<:Real}
    ScaleTransform{T}([s])
end

set!(t::ScaleTransform,ρ::Real) = t.s .= [ρ]

(t::ScaleTransform)(x) = first(t.s) .* x

_map(t::ScaleTransform, x::AbstractVector{<:Real}) = first(t.s) .* x
_map(t::ScaleTransform, x::ColVecs) = ColVecs(first(t.s) .* x.X)
_map(t::ScaleTransform, x::RowVecs) = RowVecs(first(t.s) .* x.X)

Base.isequal(t::ScaleTransform,t2::ScaleTransform) = isequal(first(t.s),first(t2.s))

Base.show(io::IO,t::ScaleTransform) = print(io,"Scale Transform (s = ", first(t.s), ")")


text:
# https://github.com/SGL-UT/GPSTk/tree/master/core/lib/GNSSCore/ObsID.hpp
# https://github.com/SGL-UT/GPSTk/tree/master/core/lib/GNSSCore/ObsIDInitializer.cpp

@enum ObservationType begin
    otUnknown
    otAny       # < Used to match any observation type
    otRange     # < pseudorange, in meters
    otPhase     # < accumulated phase, in cycles
    otDoppler   # < Doppler, in Hz
    otSNR       # < Signal strength, in dB-Hz
    otChannel   # < Channel number
    otDemodStat # < Demodulator status
    otIono      # < Ionospheric delay (see RINEX3 section 5.12)
    otSSI       # < Signal Strength Indicator (RINEX)
    otLLI       # < Loss of Lock Indicator (RINEX)
    otTrackLen  # < Number of continuous epochs of 'good' tracking
    otNavMsg    # < Navigation Message data
    otRngStdDev # < pseudorange standard deviation, in meters
    otPhsStdDev # < phase standard deviation, in cycles
    otFreqIndx  # < GLONASS frequency offset index [-6..7]
    otUndefined # < Undefined
    otLast      # < Used to verify that all items are described at compile time
end

# The frequency band this obs was collected from.
@enum CarrierBand begin
    cbUnknown
    cbAny  # < Used to match any carrier band
    cbZero # < Used with the channel observation type (see RINEx3 section 5.13)
    cbL1   # < GPS L1, Galileo E2-L1-E1, SBAS L1, QZSS L1
    cbL2   # < GPS L2, QZSS L2
    cbL5   # < GPS L5, Galileo E5a, SBAS L5, QZSS L5, INRSS L5
    cbG1   # < Glonass G1
    cbG2   # < Glonass G2
    cbG3   # < Glonass G3
    cbE5b  # < Galileo E5b, BeiDou L7
    cbE5ab # < Galileo E5a+b
    cbE6   # < Galileo E6, QZSS LEX
    cbB1   # < BeiDou L1
    cbB2   # < BeiDou L7
    cbB3   # < BeiDou L6
    cbI9   # < IRNSS S-band (RINEX '9')
    cbL1L2 # < Combined L1L2 (like an ionosphere free obs)
    cbUndefined
    cbLast # < Used to verify that all items are described at compile time
end

#= /** The code used to collect the observation. Each of these
    * should uniquely identify a code that was correlated
    * against to track the signal. While the notation generally
    * follows section 5.1 of RINEX 3, due to ambiguities in that
    * specification some extensions are made. Note that as
    * concrete specifications for the codes are released, this
    * list may need to be adjusted. Specifically, this lists
    * assumes that the same I & Q codes will be used on all
    * three of the Galileo carriers. If that is not true, more
    * identifiers need to be allocated */ =#
@enum TrackingCode begin
    tcUnknown
    tcAny     # < Used to match any tracking code
    tcCA      # < Legacy GPS civil code
    tcP       # < Legacy GPS precise code
    tcY       # < Encrypted legacy GPS precise code
    tcW       # < Encrypted legacy GPS precise code, codeless Z tracking
    tcN       # < Encrypted legacy GPS precise code, squaring codeless tracking
    tcD       # < Encrypted legacy GPS precise code, other codeless tracking
    tcM       # < Modernized GPS military unique code
    tcC2M     # < Modernized GPS L2 civil M code
    tcC2L     # < Modernized GPS L2 civil L code
    tcC2LM    # < Modernized GPS L2 civil M+L combined tracking (such as Trimble NetRS, Septrentrio, and ITT)
    tcI5      # < Modernized GPS L5 civil in-phase
    tcQ5      # < Modernized GPS L5 civil quadrature
    tcIQ5     # < Modernized GPS L5 civil I+Q combined tracking
    tcG1P     # < Modernized GPS L1C civil code tracking (pilot)
    tcG1D     # < Modernized GPS L1C civil code tracking (data)
    tcG1X     # < Modernized GPS L1C civil code tracking (pilot + data)

    tcGCA     # < Legacy Glonass civil signal
    tcGP      # < Legacy Glonass precise signal
    tcIR3     # < Glonass L3 I code
    tcQR3     # < Glonass L3 Q code
    tcIQR3    # < Glonass L3 I+Q combined tracking

    tcA       # < Galileo L1 PRS code
    tcB       # < Galileo OS/CS/SoL code
    tcC       # < Galileo Dataless code
    tcBC      # < Galileo B+C combined tracking
    tcABC     # < Galileo A+B+C combined tracking
    tcIE5     # < Galileo E5 I code
    tcQE5     # < Galileo E5 Q code
    tcIQE5    # < Galileo E5 I+Q combined tracking
    tcIE5a    # < Galileo E5a I code
    tcQE5a    # < Galileo E5a Q code
    tcIQE5a   # < Galileo E5a I+Q combined tracking
    tcIE5b    # < Galileo E5b I code
    tcQE5b    # < Galileo E5b Q code
    tcIQE5b   # < Galileo E5b I+Q combined tracking

    tcSCA     # < SBAS civil code
    tcSI5     # < SBAS L5 I code
    tcSQ5     # < SBAS L5 Q code
    tcSIQ5    # < SBAS L5 I+Q code

    tcJCA     # < QZSS civil code
    tcJD1     # < QZSS L1C(D)
    tcJP1     # < QZSS L1C(P)
    tcJX1     # < QZSS L1C(D+P)
    tcJZ1     # < QZSS L1-SAIF
    tcJM2     # < QZSS L2C(M)
    tcJL2     # < QZSS L2C(L)
    tcJX2     # < QZSS L2C(M+L)
    tcJI5     # < QZSS L5 in-phase
    tcJQ5     # < QZSS L5 quadrature
    tcJIQ5    # < QZSS L5 I+Q combined tracking
    tcJI6     # < QZSS LEX(6) in-phase
    tcJQ6     # < QZSS LEX(6) quadrature
    tcJIQ6    # < QZSS LEX(6) I+Q combined tracking

    tcCI1     # < BeiDou B1 I code
    tcCQ1     # < BeiDou B1 Q code
    tcCIQ1    # < BeiDou B1 I code
    tcCI7     # < BeiDou B2 I+Q code
    tcCQ7     # < BeiDou B2 Q code
    tcCIQ7    # < BeiDou B2 I+Q code
    tcCI6     # < BeiDou B3 I code
    tcCQ6     # < BeiDou B3 Q code
    tcCIQ6    # < BeiDou B3 I+Q code

        #  Nomenclature follows RiNEX 3.03 Table 10
    tcIA5     # < IRNSS L5 SPS
    tcIB5     # < IRNSS L5 RS(D)
    tcIC5     # < IRNSS L5 RS(P)
    tcIX5     # < IRNSS L5 B+C
    tcIA9     # < IRNSS S-band SPS
    tcIB9     # < IRNSS S=band RS(D)
    tcIC9     # < INRSS S-band RS(P)
    tcIX9     # < IRNSS S-band B+C

    tcUndefined
    tcLast    # < Used to verify that all items are described at compile time
end

# The following definitions really should only describe the items that are
# in the Rinex 3 specification. If an application needs additional ObsID
# types to be able to be translated to/from Rinex3, the additional types
# must be added by the application.

# ObsID::char2[a-z]{2}\[('.')\] = ObsID::(.+);
const char2ot = Dict{Char,ObservationType}(
        ' ' => otUnknown,
        '*' => otAny,
        'C' => otRange,
        'L' => otPhase,
        'D' => otDoppler,
        'S' => otSNR,
        '-' => otUndefined,
)

const char2cb = Dict{Char,CarrierBand}(
        ' ' => cbUnknown,
        '*' => cbAny,
        '1' => cbL1,
        '2' => cbL2,
        '3' => cbG3,
        '5' => cbL5,
        '6' => cbE6,
        '7' => cbE5b,
        '8' => cbE5ab,
        '9' => cbI9,
        '-' => cbUndefined,
)

const char2tc = Dict{Char,TrackingCode}(
        ' ' => tcUnknown,
        '*' => tcAny,
        'C' => tcCA,
        'P' => tcP,
        'W' => tcW,
        'Y' => tcY,
        'M' => tcM,
        'N' => tcN,
        'D' => tcD,
        'S' => tcC2M,
        'L' => tcC2L,
        'X' => tcC2LM,
        'I' => tcI5,
        'Q' => tcQ5,
        'A' => tcA,
        'B' => tcB,
        'Z' => tcABC,
        '-' => tcUndefined,
)

# ObsID::otDesc\[ObsID::(ot[A-Za-z]+)\]( +)= +("[A-Za-z]+");( +)//(.+)
# $1$2=>$3,$4#$5
const otDesc = Dict{ObservationType,String}(
    otUnknown => "UnknownType",   # Rinex (sp)
    otAny       => "AnyType",       # Rinex *
    otRange     => "pseudorange",   # Rinex C
    otPhase     => "phase",         # Rinex L
    otDoppler   => "doppler",       # Rinex D
    otSNR       => "snr",           # Rinex S
    otChannel   => "channel",       # Rinex
    otDemodStat => "demodStatus",   # test
    otIono      => "iono",          # Rinex
    otSSI       => "ssi",           # Rinex
    otLLI       => "lli",           # Rinex
    otTrackLen  => "tlen",          # Rinex
    otNavMsg    => "navmsg",        # Rinex
    otRngStdDev => "rngSigma",      # test
    otPhsStdDev => "phsSigma",      # test
    otFreqIndx  => "freqIndx",      # test
    otUndefined => "undefined",     # Rinex -
)
# ObsID::cbDesc\[ObsID::(cb[A-Za-z0-9]+)\]( +)= +("[A-Za-z0-9\+]+");( +)//(.+)
# $1$2=>$3,$4#$5
const cbDesc = Dict{CarrierBand,String}(
    cbUnknown   => "UnknownBand",   # Rinex (sp)
    cbAny       => "AnyBand",       # Rinex *
    cbZero      => "0",              # Rinex
    cbL1        => "L1",            # Rinex 1
    cbL2        => "L2",            # Rinex 2
    cbL5        => "L5",            # Rinex 5
    cbG1        => "G1",            # Rinex 1
    cbG2        => "G2",            # Rinex 2
    cbG3        => "G3",            # Rinex 3
    cbE5b       => "E5b",           # Rinex 7
    cbE5ab      => "E5a+b",         # Rinex 8
    cbE6        => "E6",            # Rinex 6
    cbB1        => "B1",            # Rinex 1  2 in RINEX 3.0[013]
    cbB2        => "B2",            # Rinex 7
    cbB3        => "B3",            # Rinex 6
    cbI9        => "I9",            # Rinex 9
    cbL1L2      => "comboL1L2",     # Rinex 3
    cbUndefined => "undefined",     # Rinex -
)
# ObsID::tcDesc\[ObsID::(tc[A-Za-z0-9]+)\]( +)= +(".+?");( +)//(.+)
# $1$2=>$3,$4#$5
const tcDesc = Dict{TrackingCode,String}(
    tcUnknown   => "UnknownCode",   # Rinex (sp)
    tcAny       => "AnyCode",       # Rinex *
    tcCA        => "GPSC/A",        # Rinex C    // GPScivil
    tcP         => "GPSP",          # Rinex P    // GPSprecise
    tcY         => "GPSY",          # Rinex Y    // GPSprecise_encrypted
    tcW         => "GPScodelessZ",  # Rinex W    // GPSprecise_encrypted_codeless_Z
    tcN         => "GPSsquare",     # Rinex N    // GPSprecise_encrypted_codeless_squaring
    tcD         => "GPScodeless",   # Rinex D    // GPSprecise_encrypted_codeless_other
    tcM         => "GPSM",          # Rinex M    // GPSmilitary
    tcC2M       => "GPSC2M",        # Rinex S    // GPScivil_M
    tcC2L       => "GPSC2L",        # Rinex L    // GPScivil_L
    tcC2LM      => "GPSC2L+M",      # Rinex X    // GPScivil_L+M
    tcI5        => "GPSI5",         # Rinex I    // GPScivil_I
    tcQ5        => "GPSQ5",         # Rinex Q    // GPScivil_Q
    tcIQ5       => "GPSI+Q5",       # Rinex X    // GPScivil_I+Q
    tcG1P       => "GPSC1P",        # Rinex L    // GPScivil_L1P
    tcG1D       => "GPSC1D",        # Rinex S    // GPScivil_L1D
    tcG1X       => "GPSC1(D+P)",    # Rinex X    // GPScivil_L1D+P

    tcGCA       => "GLOC/A",        # Rinex C    // GLOcivil
    tcGP        => "GLOP",          # Rinex P    // GLOprecise
    tcIR3       => "GLOIR5",        # Rinex I    // GLO L3 I code
    tcQR3       => "GLOQR5",        # Rinex Q    // GLO L3 Q code
    tcIQR3      => "GLOI+QR5",      # Rinex X    // GLO L3 I+Q code

    tcA         => "GALA",          # Rinex A    // GAL
    tcB         => "GALB",          # Rinex B    // GAL
    tcC         => "GALC",          # Rinex C    // GAL
    tcBC        => "GALB+C",        # Rinex X    // GAL
    tcABC       => "GALA+B+C",      # Rinex Z    // GAL
    tcIE5       => "GALI5",         # Rinex I    // GAL
    tcQE5       => "GALQ5",         # Rinex Q    // GAL
    tcIQE5      => "GALI+Q5",       # Rinex X    // GAL
    tcIE5a      => "GALI5a",        # Rinex I    // GAL
    tcQE5a      => "GALQ5a",        # Rinex Q    // GAL
    tcIQE5a     => "GALI+Q5a",      # Rinex X    // GAL
    tcIE5b      => "GALI5b",        # Rinex I    // GAL
    tcQE5b      => "GALQ5b",        # Rinex Q    // GAL
    tcIQE5b     => "GALI+Q5b",      # Rinex X    // GAL

    tcSCA       => "SBASC/A",       # Rinex C    // SBAS civil code
    tcSI5       => "SBASI5",        # Rinex I    // SBAS L5 I code
    tcSQ5       => "SBASQ5",        # Rinex Q    // SBAS L5 Q code
    tcSIQ5      => "SBASI+Q5",      # Rinex X    // SBAS L5 I+Q code

    tcJCA       => "QZSSC/A",       # Rinex C    // QZSS L1 civil code
    tcJD1       => "QZSSL1C(D)",    # Rinex S    // QZSS L1C(D)
    tcJP1       => "QZSSL1C(P)",    # Rinex L    // QZSS L1C(P)
    tcJX1       => "QZSSL1C(D+P)",  # Rinex X    // QZSS L1C(D+P)
    tcJZ1       => "QZSSL1-SAIF",   # Rinex Z    // QZSS L1-SAIF
    tcJM2       => "QZSSL2C(M)",    # Rinex S    // QZSS L2 M code
    tcJL2       => "QZSSL2C(L)",    # Rinex L    // QZSS L2 L code
    tcJX2       => "QZSSL2C(M+L)",  # Rinex X    // QZSS L2 M+L code
    tcJI5       => "QZSSL5I",       # Rinex I    // QZSS L5 I code
    tcJQ5       => "QZSSL5Q",       # Rinex Q    // QZSS L5 Q code
    tcJIQ5      => "QZSSL5I+Q",     # Rinex X    // QZSS L5 I+Q code
    tcJI6       => "QZSSL6I",       # Rinex S    // QZSS LEX(6) I code
    tcJQ6       => "QZSSL6Q",       # Rinex L    // QZSS LEX(6) Q code
    tcJIQ6      => "QZSSL6I+Q",     # Rinex X    // QZSS LEX(6) I+Q code

    tcCI1       => "BDSIB1",        # Rinex I    // BeiDou L1 I code
    tcCQ1       => "BDSQB1",        # Rinex Q    // BeiDou L1 Q code
    tcCIQ1      => "BDSI+QB1",      # Rinex X    // BeiDou L1 I+Q code
    tcCI7       => "BDSIB2",        # Rinex I    // BeiDou B2 I code
    tcCQ7       => "BDSQB2",        # Rinex Q    // BeiDou B2 Q code
    tcCIQ7      => "BDSI+QB2",      # Rinex X    // BeiDou B2 I+Q code
    tcCI6       => "BDSIB3",        # Rinex I    // BeiDou B3 I code
    tcCQ6       => "BDSQB3",        # Rinex Q    // BeiDou B3 Q code
    tcCIQ6      => "BDSI+QB3",      # Rinex X    // BeiDou B3 I+Q code

    tcIA5       => "IRNSSL5A",      # Rinex A    // IRNSS L5 SPS
    tcIB5       => "IRNSSL5B",      # Rinex B    // IRNSS L5 RS(D)
    tcIC5       => "IRNSSL5C",      # Rinex C    // IRNSS L5 RS(P)
    tcIX5       => "IRNSSL5B+C",    # Rinex X    // IRNSS L5 B+C
    tcIA9       => "IRNSSL9A",      # Rinex A    // IRNSS S-band SPS
    tcIB9       => "IRNSSL9B",      # Rinex B    // IRNSS S-band RS(D)
    tcIC9       => "IRNSSL9C",      # Rinex C    // IRNSS S-band RS(P)
    tcIX9       => "IRNSSL9B+C",    # Rinex X    // IRNSS S-band B+C
    tcUndefined => "undefined",     # test
)

export ObsID
struct ObsID
    str::String
    ot::ObservationType
    cb::CarrierBand
    tc::TrackingCode
    ss::SatelliteSystemType
    function ObsID(str::Union{String,SubString}, sys::SatelliteSystemType=sstUnknown)
        length(str) == 3  || throw(error("ObsID |$str| error"))
        cot, ccb, ctc = str
        haskey(char2cb, ccb) || throw(error("no such CarrierBand: $ccb"))
        haskey(char2ot, cot) || throw(error("no such ObservationType: $cot"))
        haskey(char2tc, ctc) || throw(error("no such TrackingCode: $ctc"))
        type = char2ot[cot]
        band = char2cb[ccb]
        code = char2tc[ctc]
        # / This next block takes care of fixing up the codes that are reused
        # / between the various signals
        if sys == sstGPS
            if band == cbL5
                if code == tcC2LM code = tcIQ5 end
            elseif band == cbL1
                if (code == tcC2LM)     code = tcG1X
                elseif (code == tcC2M)  code = tcG1D
                elseif (code == tcC2L)  code = tcG1P
                end
            end
        elseif sys == sstGalileo
            if band == cbL1 || band == cbE6
                if code == tcCA       code = tcC
                elseif code == tcC2LM code = tcBC
                end
            elseif band == cbL5
                if code == tcI5       code = tcIE5a
                elseif code == tcQ5   code = tcQE5a
                elseif code == tcC2LM code = tcIQE5a
                end
            elseif band == cbE5b
                if code == tcI5       code = tcIE5b
                elseif code == tcQ5   code = tcQE5b
                elseif code == tcC2LM code = tcIQE5b
                end
            elseif band == cbE5ab
                if code == tcI5       code = tcIE5
                elseif code == tcQ5   code = tcQE5
                elseif code == tcC2LM code = tcIQE5
                end
            end
        elseif sys == sstGlonass
            if code == tcCA      code = tcGCA
            elseif code == tcP   code = tcGP
            elseif code == tcI5  code = tcIR3
            elseif code == tcQ5  code = tcQR3
            elseif code == tcC2LM || code == tcG1X
                code = tcIQR3
            end
            if band == cbL1 band = cbG1
            elseif band == cbL2 band = cbG2
            end

        elseif sys == sstGeosync
            if code == tcCA     code = tcSCA     # 'C'
            elseif code == tcI5 code = tcSI5     # 'I'
            elseif code == tcQ5 code = tcSQ5     # 'Q'
            elseif code == tcC2LM || code == tcG1X
                code = tcSIQ5  # 'X'
            end
        elseif sys == sstQZSS
            if band == cbL1
                if code == tcCA
                    code = tcJCA    # 'C'
                elseif code == tcC2M  || code == tcG1D
                    code = tcJD1    # 'S'
                elseif code == tcC2L  || code == tcG1P
                    code = tcJP1    # 'L'
                elseif code == tcC2LM || code == tcG1X
                    code = tcJX1    # 'X'
                elseif code == tcABC
                    code = tcJZ1    # 'Z'
                end
            elseif band == cbL2
                if code == tcC2M      || code == tcG1D
                    code = tcJM2  # 'S'
                elseif code == tcC2L  || code == tcG1P
                    code = tcJL2  # 'L'
                elseif code == tcC2LM || code == tcG1X
                    code = tcJX2  # 'X'
                end
            elseif band == cbL5
                if code == tcI5       code = tcJI5    # 'I'
                elseif code == tcQ5   code = tcJQ5    # 'Q'
                elseif code == tcC2LM code = tcJIQ5   # 'X'
                end
            elseif band == cbE6
                if code == tcC2M  || code == tcG1D
                    code = tcJI6     # 'S'
                elseif code == tcC2L  || code == tcG1P
                    code = tcJQ6     # 'L'
                elseif code == tcC2LM || code == tcG1X
                    code = tcJIQ6    # 'X'
                end
            end
        elseif sys == sstBeiDou
            if band == cbL1     band = cbB1         # RINEX 3.02
            elseif band == cbL2 band = cbB1         # RINEX 3.0[013]
            elseif band == cbE6 band = cbB3
            elseif band == cbB1
                if code == tcI5     code = tcCI1    # 'I'
                elseif code == tcQ5 code = tcCQ1    # 'Q'
                elseif code == tcC2LM || code == tcG1X
                    code = tcCIQ1                   # 'X'
                end
            elseif band == cbB3
                if code == tcI5      code = tcCI7    # 'I'
                elseif code == tcQ5  code = tcCQ7    # 'Q'
                elseif code == tcC2LM || code == tcG1X
                    code = tcCIQ7                    # 'X'
                end
            elseif band == cbE5b
                if code == tcI5     code = tcCI6     # 'I'
                elseif code == tcQ5 code = tcCQ6     # 'Q'
                elseif code == tcC2LM || code == tcG1X
                    code = tcCIQ6                    # 'X'
                end
            end
        elseif sys == sstIRNSS # IRNSS
            if band == cbL5
                if code == tcCA     code = tcIA5   # 'A'
                elseif code == tcA  code = tcIB5   # 'B'
                elseif code == tcB  code = tcIC5   # 'B'
                elseif code == tcC2LM || code == tcG1X
                    code = tcIX5                   # 'X'
                end
            end
        end
        new(str,  type, band, code, sys)
    end
end


text:
using Primes

function carmichael(pmax::Integer)
    if pmax ≤ 0 throw(DomainError("pmax must be strictly positive")) end
    car = Vector{typeof(pmax)}(0)
    for p in primes(pmax)
        for h₃ in 2:(p-1)
            m = (p - 1) * (h₃ + p)
            pmh = mod(-p ^ 2, h₃)
            for Δ in 1:(h₃+p-1)
                if m % Δ != 0 || Δ % h₃ != pmh continue end
                q = m ÷ Δ + 1
                if !isprime(q) continue end
                r = (p * q - 1) ÷ h₃ + 1
                if !isprime(r) || mod(q * r, p - 1) == 1 continue end
                append!(car, [p, q, r])
            end
        end
    end
    return reshape(car, 3, length(car) ÷ 3)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
for n in [:DateTime, :Date, :Time, :TimeDelta, :TZInfo, :TimeZone]
    p = lowercase(string(n))
    t = Symbol(:Py, n, :_Type)
    r = Symbol(t, :__ref)
    c = Symbol(:Py, n, :_Check)
    @eval $r = Ref(PyNULL)
    @eval $t(doimport::Bool = true) = begin
        ptr = $r[]
        isnull(ptr) || return ptr
        m = Py_DateTimeModule(doimport)
        isnull(m) && return PyNULL
        o = PyObject_GetAttrString(m, $p)
        isnull(o) && return PyNULL
        $r[] = o
    end
    @eval $c(o) = begin
        t = $t(false)
        isnull(t) && return (PyErr_IsSet() ? Cint(-1) : Cint(0))
        PyObject_IsInstance(o, t)
    end
end

PyDateTime_FromParts(
    year::Integer = 1,
    month::Integer = 1,
    day::Integer = 1,
    hour::Integer = 0,
    minute::Integer = 0,
    second::Integer = 0,
    microsecond::Integer = 0;
    tzinfo = nothing,
) = begin
    t = PyDateTime_Type()
    isnull(t) && return PyNULL
    PyObject_CallNice(
        t,
        year,
        month,
        day,
        hour,
        minute,
        second,
        microsecond,
        tzinfo,
    )
end

PyDateTime_From(x::DateTime) = PyDateTime_FromParts(
    year(x),
    month(x),
    day(x),
    hour(x),
    minute(x),
    second(x),
    millisecond(x) * 1000,
)
PyDateTime_From(x::Date) = PyDateTime_FromParts(year(x), month(x), day(x))

PyDateTime_IsAware(o::PyPtr) = begin
    tzinfo = PyObject_GetAttrString(o, "tzinfo")
    isnull(tzinfo) && return Cint(-1)
    PyNone_Check(tzinfo) && (Py_DecRef(tzinfo); return Cint(0))
    utcoffset = PyObject_GetAttrString(tzinfo, "utcoffset")
    Py_DecRef(tzinfo)
    isnull(utcoffset) && return Cint(-1)
    off = PyObject_CallNice(utcoffset, PyObjectRef(o))
    Py_DecRef(utcoffset)
    isnull(off) && return Cint(-1)
    r = PyNone_Check(off)
    Py_DecRef(off)
    r ? Cint(0) : Cint(1)
end

PyDateTime_TryConvertRule_datetime(o::PyPtr, ::Type{S}) where {S<:DateTime} = begin
    # TODO: worry about fold?
    # can only convert non-aware times
    aware = PyDateTime_IsAware(o)
    aware == -1 && return -1
    aware != 0 && return 0
    # year
    p = PyObject_GetAttrString(o, "year")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    year = takeresult(Int)
    # month
    p = PyObject_GetAttrString(o, "month")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    month = takeresult(Int)
    # day
    p = PyObject_GetAttrString(o, "day")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    day = takeresult(Int)
    # hour
    p = PyObject_GetAttrString(o, "hour")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    hour = takeresult(Int)
    # minute
    p = PyObject_GetAttrString(o, "minute")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    min = takeresult(Int)
    # second
    p = PyObject_GetAttrString(o, "second")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    sec = takeresult(Int)
    # microsecond
    p = PyObject_GetAttrString(o, "microsecond")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    micro = takeresult(Int)
    # millisecond resolution
    mod(micro, 1000) == 0 || return 0
    # success
    putresult(DateTime(year, month, day, hour, min, sec, div(micro, 1000)))
end

PyDate_FromParts(year::Integer = 1, month::Integer = 1, day::Integer = 1) = begin
    t = PyDate_Type()
    isnull(t) && return PyNULL
    PyObject_CallNice(t, year, month, day)
end

PyDate_From(x::Date) = PyDate_FromParts(year(x), month(x), day(x))

PyDate_TryConvertRule_date(o::PyPtr, ::Type{S}) where {S<:Date} = begin
    # year
    p = PyObject_GetAttrString(o, "year")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    year = takeresult(Int)
    # month
    p = PyObject_GetAttrString(o, "month")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    month = takeresult(Int)
    # day
    p = PyObject_GetAttrString(o, "day")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    day = takeresult(Int)
    # success
    putresult(Date(year, month, day))
end

PyTime_FromParts(
    hour::Integer = 0,
    minute::Integer = 0,
    second::Integer = 0,
    microsecond::Integer = 0;
    tzinfo = nothing,
) = begin
    t = PyTime_Type()
    isnull(t) && return PyNULL
    PyObject_CallNice(t, hour, minute, second, microsecond, tzinfo)
end

PyTime_From(x::Time) =
    if iszero(nanosecond(x))
        PyTime_FromParts(
            hour(x),
            minute(x),
            second(x),
            millisecond(x) * 1000 + microsecond(x),
        )
    else
        PyErr_SetString(
            PyExc_ValueError(),
            "cannot create 'datetime.time' with resolution less than microseconds",
        )
        PyNULL
    end

PyTime_IsAware(o::PyPtr) = begin
    tzinfo = PyObject_GetAttrString(o, "tzinfo")
    isnull(tzinfo) && return Cint(-1)
    PyNone_Check(tzinfo) && (Py_DecRef(tzinfo); return Cint(0))
    utcoffset = PyObject_GetAttrString(tzinfo, "utcoffset")
    Py_DecRef(tzinfo)
    isnull(utcoffset) && return Cint(-1)
    off = PyObject_CallNice(utcoffset, PyObjectRef(Py_None()))
    Py_DecRef(utcoffset)
    isnull(off) && return Cint(-1)
    r = PyNone_Check(off)
    Py_DecRef(off)
    r ? Cint(0) : Cint(1)
end

PyTime_TryConvertRule_time(o::PyPtr, ::Type{S}) where {S<:Time} = begin
    # TODO: worry about fold?
    # can only convert non-aware times
    aware = PyTime_IsAware(o)
    aware == -1 && return -1
    aware != 0 && return 0
    # hour
    p = PyObject_GetAttrString(o, "hour")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    hour = takeresult(Int)
    # minute
    p = PyObject_GetAttrString(o, "minute")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    min = takeresult(Int)
    # second
    p = PyObject_GetAttrString(o, "second")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    sec = takeresult(Int)
    # microsecond
    p = PyObject_GetAttrString(o, "microsecond")
    isnull(p) && return -1
    r = PyObject_TryConvert(p, Int)
    Py_DecRef(p)
    r < 1 && return r
    micro = takeresult(Int)
    # success
    putresult(Time(hour, min, sec, div(micro, 1000), mod(micro, 1000)))
end

PyTimeDelta_FromParts(;
    days = 0,
    seconds = 0,
    microseconds = 0,
    milliseconds = 0,
    minutes = 0,
    hours = 0,
    weeks = 0,
) = begin
    t = PyTimeDelta_Type()
    isnull(t) && return PyNULL
    PyObject_CallNice(t, days, seconds, microseconds, milliseconds, minutes, hours, weeks)
end

PyTimeDelta_From(x::Second) = PyTimeDelta_FromParts(seconds = Dates.value(x))
PyTimeDelta_From(x::Millisecond) = PyTimeDelta_FromParts(milliseconds = Dates.value(x))
PyTimeDelta_From(x::Microsecond) = PyTimeDelta_FromParts(microseconds = Dates.value(x))
PyTimeDelta_From(x::Nanosecond) =
    if mod(Dates.value(x), 1000) == 0
        PyTimeDelta_FromParts(microseconds = div(Dates.value(x), 1000))
    else
        PyErr_SetString(
            PyExc_ValueError(),
            "cannot create 'datetime.timedelta' with resolution less than microseconds",
        )
        PyNULL
    end
PyTimeDelta_From(x::Union{Minute,Hour,Day,Week,Month,Year}) = begin
    PyErr_SetString(
        PyExc_ValueError(),
        "cannot create 'datetime.timedelta' from a Julia '$(typeof(x))' because the latter does not specify a fixed period of time",
    )
    PyNULL
end
PyTimeDelta_From(x::Period) = begin
    PyErr_SetString(
        PyExc_ValueError(),
        "cannot create 'datetime.timedelta' from a Julia '$(typeof(x))'",
    )
    PyNULL
end
PyTimeDelta_From(x::Dates.CompoundPeriod) = begin
    r = PyTimeDelta_FromParts()
    isnull(r) && return PyNULL
    for p in x.periods
        po = PyTimeDelta_From(p)
        isnull(po) && (Py_DecRef(r); return PyNULL)
        r2 = PyNumber_Add(r, po)
        Py_DecRef(po)
        Py_DecRef(r)
        isnull(r2) && return PyNULL
        r = r2
    end
    r
end

PyTimeDelta_AsParts(o::PyPtr) = begin
    # days
    p = PyObject_GetAttrString(o, "days")
    isnull(p) && return (-1, -1, -1)
    r = PyObject_Convert(p, Int)
    Py_DecRef(p)
    r == -1 && return (-1, -1, -1)
    days = takeresult(Int)
    # seconds
    p = PyObject_GetAttrString(o, "seconds")
    isnull(p) && return (-1, -1, -1)
    r = PyObject_Convert(p, Int)
    Py_DecRef(p)
    r == -1 && return (-1, -1, -1)
    seconds = takeresult(Int)
    # microseconds
    p = PyObject_GetAttrString(o, "microseconds")
    isnull(p) && return (-1, -1, -1)
    r = PyObject_Convert(p, Int)
    Py_DecRef(p)
    r == -1 && return (-1, -1, -1)
    microseconds = takeresult(Int)
    return (days, seconds, microseconds)
end

PyTimeDelta_TryConvertRule_compoundperiod(
    o::PyPtr,
    ::Type{S},
) where {S<:Dates.CompoundPeriod} = begin
    # unpack
    days, seconds, microseconds = PyTimeDelta_AsParts(o)
    microseconds == -1 && return -1
    # check overflow
    (
        cld(typemin(Int), 60 * 60 * 24) ≤
        days ≤
        fld(typemax(Int) - (60 * 60 * 24 - 1), 60 * 60 * 24)
    ) || return 0
    # make result
    putresult(
        S(
            Period[
                Second(seconds + 60 * 60 * 24 * days),
                Millisecond(div(microseconds, 1000)),
                Microsecond(mod(microseconds, 1000)),
            ],
        ),
    )
end

PyTimeDelta_TryConvertRule_period(o::PyPtr, ::Type{S}) where {S<:Period} = begin
    # unpack
    days, seconds, microseconds = PyTimeDelta_AsParts(o)
    microseconds == -1 && return -1
    # checks for overflow
    overflowcheck(x, B) = cld(typemin(Int), B) ≤ x ≤ fld(typemax(Int) - (B - 1), B)
    # make result
    if Microsecond <: S && overflowcheck(days, 1000_000 * 60 * 60 * 24)
        putresult(
            Microsecond(microseconds + 1000_000 * seconds + 1000_000 * 60 * 60 * 24 * days),
        )
    elseif Nanosecond <: S && overflowcheck(days, 1000_000_000 * 60 * 60 * 24)
        putresult(
            Nanosecond(
                1000 * microseconds +
                1000_000_000 * seconds +
                1000_000_000 * 60 * 60 * 24 * days,
            ),
        )
    elseif Millisecond <: S &&
           mod(microseconds, 1000) == 0 &&
           overflowcheck(days, 1000 * 60 * 60 * 24)
        putresult(
            Millisecond(
                div(microseconds, 1000) + 1000 * seconds + 1000 * 60 * 60 * 24 * days,
            ),
        )
    elseif Second <: S && microseconds == 0 && overflowcheck(days, 60 * 60 * 24)
        putresult(Second(seconds + 60 * 60 * 24 * days))
    else
        0
    end
end


text:
import JLD.readas
import JLD.writeas

#=
Assoc Serialized for saving
Note that saving would convert row and col types to number.
=#

type AssocSerial
    rowstr::AbstractString
    colstr::AbstractString
    valstr::AbstractString

    rownum::AbstractString
    colnum::AbstractString
    valnum::AbstractString

    A::AbstractSparseMatrix
end

function JLD.writeas(data::Assoc)
    #Get parts from the assoc
    row = Row(data)
    col = Col(data)
    val = Val(data)
    A   = Adj(data)

    #Split the mapping array into string and number arraries for storage
    rowindex = searchsortedfirst(row,AbstractString,lt=isa)
    colindex = searchsortedfirst(col,AbstractString,lt=isa)
    valindex = searchsortedfirst(val,AbstractString,lt=isa)
    #test if there are string elements in the array
    if (rowindex == 1)
        rowstr = ""
    else
        rowstr = join(row[1:(rowindex-1)],",");
    end

    if (colindex == 1)
        colstr = ""
    else
        colstr = join(col[1:(colindex-1)],",");
    end

    if (valindex == 1)
        valstr = ""
    else
        valstr = join(val[1:(valindex-1)],",");
    end

    #test if there are number elements in the array.  Note that by serialization all numbers will be converted to float.

    if (rowindex == size(row,1)+1)
        rownum = ""
    else
        rownum = join(row[rowindex:end],",");
    end

    if (colindex == size(col,1)+1)
        colnum = ""
    else
        colnum = join(col[colindex:end],",");
    end

    if (valindex == size(val,1)+1)
        valnum = ""
    else
        valnum = join(val[valindex:end],",");
    end

    #Reconstitute converted parts into Serialized

    return AssocSerial(rowstr,colstr,valstr,rownum,colnum,valnum,A)
end

function JLD.readas(serData::AssocSerial)
    row = []
    if serData.rowstr != ""
       row = vcat(row, split(serData.rowstr,",")) ;
    end
    if serData.rownum != ""
        row = vcat(row, map(float,split(serData.rownum,",")));
    end
    row = Array{Union{AbstractString,Number}}(row)

    col = []
    if serData.colstr != ""
       col = vcat(col, split(serData.colstr,",")) ;
    end
    if serData.colnum != ""
        col = vcat(col, map(float,split(serData.colnum,",")));
    end
    col = Array{Union{AbstractString,Number}}(col)

    val = []
    if serData.valstr != ""
       val = vcat(val, split(serData.valstr,",")) ;
    end
    if serData.valnum != ""
        val = vcat(val, map(float,split(serData.valnum,",")));
    end
    val = Array{Union{AbstractString,Number}}(val)

    return Assoc(row,col,val,serData.A)
end


text:
using Replay

instructions = """
using TestImages
using ImageInTerminal
testimage("mandril_color")
"""

replay(instructions, julia_project=@__DIR__)


text:
# =========================================================================================
"""
    readmcsf{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64;
        # kwargs
        domain::Symbol=:none
    )

Reads a volume model from the given GAMer-generated mcsf file.

!!! note
    This file type does not support charge models! Hence, the charge list of the returning
    `Model` object is empty and has to be set separately.

# Arguments
 * `domain` Element domain (solute `:Ω`, solvent `:Σ`, or `:none`)

# Return type
`Model{T, Tetrahedron{T}}`

# Aliases

    readmcsf{T}(
        fname ::String,
              ::Type{T}=Float64;
        # kwargs
        domain::Symbol=:none
    )

Reads the model using a file name rather than a `IOStream` object.

    readmcsf{T}(
        fnameΩ::String,
        fnameΣ::String,
              ::Type{T}=Float64
    )

Reads the contents of two separate files (given by name), sets the element domains to `:Ω`
or `:Σ`, respectively, and returns a single (merged) `Model` object.
"""
function readmcsf(
        stream::IOStream,
              ::Type{T}=Float64;
        domain::Symbol=:none
    ) where T <: AbstractFloat
    nodes = readmcsf_nodes(stream, T)
    Model(nodes, readmcsf_elements(stream, nodes, domain=domain))
end

function readmcsf(
        fname ::String,
              ::Type{T}=Float64;
        domain::Symbol=:none
    ) where T
    open(fh -> readmcsf(fh, T, domain=domain), fname)
end

function readmcsf(
        fnameΩ::String,
        fnameΣ::String,
              ::Type{T}=Float64
    ) where T
    meshunion(readmcsf(fnameΩ, T, domain=:Ω), readmcsf(fnameΣ, T, domain=:Σ))
end

# =========================================================================================
"""
    readmcsf_nodes{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64
    )

Reads all nodes from the given GAMer-generated mcsf file.

# Return type
`Vector{Vector{T}}`
"""
function readmcsf_nodes(
        stream::IOStream,
              ::Type{T}=Float64
    ) where T <: AbstractFloat
    nodes = Vector{T}[]
    seek(stream, "vert=[")
    for line in eachline(stream)
        startswith(line, "%") && continue        # skip comments
        startswith(line, "];") && break          # all nodes read
        push!(nodes, [parse(T, e) for e in split(line)[3:5]])
    end
    nodes
end

# =========================================================================================
"""
    readmcsf_elements{T <: AbstractFloat}(
        stream::IOStream,
        nodes ::Vector{Vector{T}};
        # kwargs
        domain::Symbol=:none
    )

Reads all elements from the given GAMer-generated mcsf file.

# Return type
`Vector{Tetrahedron{T}}`
"""
function readmcsf_elements(
        stream::IOStream,
        nodes ::Vector{Vector{T}};
        domain::Symbol=:none
    ) where T <: AbstractFloat
    elements = Tetrahedron{T}[]
    seek(stream, "simp=[")
    for line in eachline(stream)
        startswith(line, "%") && continue        # skip comments
        startswith(line, "];") && break          # all simplices read
        push!(
            elements,
            Tetrahedron([nodes[parse(Int, e)+1] for e in split(line)[8:end]]..., domain)
        )
    end
    elements
end


text:
export hasweights, weights
export setscope!, setobs!, addobs!, scope, obs
export nobs, nscope
export isnormalized
export addscope!
export removescope!
export hasscope, hasobs
export logweights
export updatescope!
export classes, children, parents, length, add!, remove!, logpdf!, logpdf
export rand

function isnormalized(node::Node)
    if !hasweights(node)
        return mapreduce(child -> isnormalized(child), &, children(node))
    else
        return sum(weights(node)) ≈ 1.0
    end
end
isnormalized(node::Leaf) = true

hasweights(node::SumNode) = true
hasweights(node::Node) = false

weights(node::SumNode) = exp.(node.logweights)
logweights(node::SumNode) = node.logweights

getindex(node::Node, i...) = getindex(node.children, i...)

function setscope!(node::SPNNode, scope::Vector{Int})
    if length(scope) > 0

        if maximum(scope) > length(node.scopeVec)
            @warn "New scope is larger than original scope, resize node scope..."
            resize!(node.scopeVec, maximum(scope))
        end

        fill!(node.scopeVec, false)
        node.scopeVec[scope] .= true
    else
        fill!(node.scopeVec, false)
    end
end

function setscope!(node::SPNNode, scope::Int)

    if scope > length(node.scopeVec)
        @warn "New scope is larger than original scope, resize node scope..."
        resize!(node.scopeVec, scope)
    end

    fill!(node.scopeVec, false)
    node.scopeVec[scope] = true
end

function addscope!(node::Node, scope::Int)
    @assert scope <= length(node.scopeVec)
    node.scopeVec[scope] = true
end

function removescope!(node::Node, scope::Int)
    @assert scope <= length(node.scopeVec)
    node.scopeVec[scope] = false
end

@inline scope(node::Node) = findall(node.scopeVec)
@inline scope(node::Leaf) = node.scope
@inline nscope(node::Node) = sum(node.scopeVec)
@inline nscope(node::Leaf) = length(node.scope)
@inline hasscope(node::Node) = sum(node.scopeVec) > 0
@inline hasscope(node::Leaf) = true
@inline hassubscope(n1::ProductNode, n2::SPNNode) = any(n1.scopeVec .& n2.scopeVec)

function addobs!(node::SPNNode, obs::Int)
    @assert obs <= length(node.obsVec)
    node.obsVec[obs] = true
end

function setobs!(node::Node, obs::AbstractVector{Int})
    if length(obs) > 0
        if maximum(obs) > length(node.obsVec)
            @warn "New obs is larger than original obs field, resize node obs..."
            resize!(node.obsVec, maximum(obs))
        end

        fill!(node.obsVec, false)
        node.obsVec[obs] .= true
    else
        fill!(node.obsVec, false)
    end
end

function setobs!(node::Node, obs::Int)
    if obs > length(node.obsVec)
        @warn "New obs is larger than original obs field, resize node obs..."
        resize!(node.obsVec, obs)
    end

    fill!(node.obsVec, false)
    node.obsVec[obs] = true
end

@inline nobs(node::Node) = sum(node.obsVec)
@inline obs(node::Node) = findall(node.obsVec)
@inline hasobs(node::Node) = any(node.obsVec)

"""
    updatescope!(spn)
Update the scope of all nodes in the SPN.
"""
updatescope!(spn::SumProductNetwork) = updatescope!(spn.root)

function updatescope!(node::SumNode)
    for child in children(node)
        updatescope!(child)
    end
    setscope!(node, scope(first(children(node))))
    return node
end

function updatescope!(node::ProductNode)
    for child in children(node)
        updatescope!(child)
    end
    setscope!(node, mapreduce(c -> scope(c), vcat, children(node)))
    return node
end

updatescope!(node::Leaf) = node

"""

children(node) -> children::SPNNode[]

Returns the children of an internal node.

##### Parameters:
* `node::Node`: Internal SPN node to be evaluated.
"""
function children(node::Node)
    node.children
end

"""

parents(node) -> parents::SPNNode[]

Returns the parents of an SPN node.

##### Parameters:
* `node::SPNNode`: SPN node to be evaluated.
"""
function parents(node::SPNNode)
    node.parents
end

"""
Add a node to a finite sum node with given weight in place.
add!(node::FiniteSumNode, child::SPNNode, weight<:Real)
"""
function add!(parent::SumNode, child::SPNNode, logweight::T) where T <: Real
    if !(child in parent.children)
        push!(parent.children, child)
        push!(parent.logweights, logweight)
        push!(child.parents, parent)
    end
end

function add!(parent::ProductNode, child::SPNNode)
    if !(child in parent.children)
        push!(parent.children, child)
        push!(child.parents, parent)
    end
end

"""
Remove a node from the children list of a sum node in place.
remove!(node::FiniteSumNode, index::Int)
"""
function remove!(parent::SumNode, index::Int)
    pid = findfirst(map(p -> p == parent, parent.children[index].parents))

    @assert pid > 0 "Could not find parent ($(node.id)) in list of parents ($(parent.children[index].parents))!"
    deleteat!(parent.children[index].parents, pid)
    deleteat!(parent.children, index)
    deleteat!(parent.logweights, index)
end

function remove!(parent::ProductNode, index::Int)
    pid = findfirst(parent .== parent.children[index].parents)
    @assert pid > 0 "Could not find parent ($(node.id)) in list of parents ($(parent.children[index].parents))!"
    deleteat!(parent.children[index].parents, pid)
    deleteat!(parent.children, index)
end

function length(node::Node)
    Base.length(node.children)
end

function logpdf!(spn::SumProductNetwork, X::AbstractMatrix{<:Real}, llhvals::AxisArray)

    fill!(llhvals, -Inf)

    # Call inplace functions.
    for layer in spn.layers
        Threads.@threads for node in layer
            logpdf!(node, X, llhvals)
        end
    end

    spn.info[:llh] = mean(llhvals[:,spn.root.id])

    return llhvals[:,spn.root.id]
end

function logpdf(spn::SumProductNetwork, X::AbstractMatrix{<:Real})
    idx = Axis{:id}(collect(keys(spn)))
    llhvals = AxisArray(Matrix{Float64}(undef, size(X, 1), length(idx)), 1:size(X, 1), idx)
    return logpdf!(spn, X, llhvals)
end

function logpdf!(spn::SumProductNetwork, x::AbstractVector{<:Real}, llhvals::AxisArray)

    fill!(llhvals, -Inf)

    # Call inplace functions.
    for layer in spn.layers
        Threads.@threads for node in layer
            logpdf!(node, x, llhvals)
        end
    end

    spn.info[:llh] = llhvals[spn.root.id]

    return llhvals[spn.root.id]
end

function logpdf(spn::SumProductNetwork, x::AbstractVector{<:Real})
    idx = Axis{:id}(collect(keys(spn)))
    llhvals = AxisArray(Vector{Float64}(undef, length(idx)), idx)
    return logpdf!(spn, x, llhvals)
end

"""
    logpdf(n::SumNode, x)

Compute the logpdf of a sum node.
"""
logpdf(n::SumNode, x::AbstractArray{<:Real}; lw::AbstractVector{<:Real}=logweights(n)) = n(x, lw=lw)

function logpdf!(n::SumNode, x::AbstractVector{<:Real}, llhvals::AxisArray{U,1}) where {U<:Real}
    @inbounds y = map(c -> llhvals[c.id], children(n))
    llhvals[n.id] = map(U, _logpdf(n, x, y, logweights(n)))
    return llhvals
end

function logpdf!(n::SumNode, x::AbstractMatrix{<:Real}, llhvals::AxisArray{U,2}) where {U<:Real}
    @inbounds Y = mapreduce(c -> llhvals[:, c.id], hcat, children(n))
    @inbounds llhvals[:,n.id] = map(i -> map(U, _logpdf(n, x, view(Y, i, :), logweights(n))), 1:size(x,1))
    return llhvals
end

# ############################################ #
# Concrete implementation for finite sum node. #
# ############################################ #
function (n::FiniteSumNode)(x::AbstractArray{<:Real}, y::AbstractVector{<:Real}; lw::AbstractVector{<:Real}=logweights(n))
    return _logpdf(n, x, y, lw)
end

function _logpdf(n::SumNode, x::AbstractArray{<:Real}, y::AbstractVector{<:Real}, lw::AbstractVector{<:Real})
    l = lw+y
    m = maximum(l)
    lp = log(mapreduce(v -> exp(v-m), +, l)) + m
    return isfinite(lp) ? lp : -Inf
end

function (n::FiniteSumNode)(x::AbstractVector{<:Real}; lw::AbstractVector{<:Real} = logweights(n))
    y = map(c -> c(x), children(n))
    return _logpdf(n, x, y, lw)
end

function (n::FiniteSumNode)(x::AbstractMatrix{<:Real}; lw::AbstractVector{<:Real}=logweights(n))
    @inbounds Y = mapreduce(i -> map(c -> c(view(x, i, :)), children(n)), hcat, 1:size(x,1))
    return @inbounds map(i -> _logpdf(n, x, view(Y, :, i), lw), 1:size(x,1))
end

"""
    logpdf(n::ProductNode, x)

Compute the logpdf of a product node.
"""
logpdf(n::ProductNode, x::AbstractArray{<:Real}) = n(x)
(n::FiniteProductNode)(x::AbstractArray{<:Real}, y::AbstractVector{<:Real}) = _logpdf(n, x, y)
function (n::FiniteProductNode)(x::AbstractVector{<:Real})
    y = map(c -> c(x), children(n))
    return _logpdf(n, x, y)
end

function (n::FiniteProductNode)(x::AbstractMatrix{<:Real})
    @inbounds Y = mapreduce(i -> map(c -> c(view(x, i, :)), children(n)), hcat, 1:size(x,1))
    return _logpdf(n, x, Y)
end

function _logpdf(n::ProductNode, x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    if !hasscope(n)
        return 0.0
    else
        return mapreduce(k -> hasscope(n[k]) ? y[k] : 0.0, +, 1:length(n))
    end
end

function _logpdf(n::ProductNode, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real})
    N = size(x,1)
    if !hasscope(n)
        return zeros(N)
    else
        return mapreduce(k -> hasscope(n[k]) ? y[:,k] : zeros(N), +, 1:length(n))
    end
end

function logpdf!(n::ProductNode, x::AbstractVector{<:Real}, llhvals::AxisArray{U}) where {U<:Real}
    @inbounds y = map(c -> llhvals[c.id], children(n))
    llhvals[n.id] = map(U, _logpdf(n, x, y))
    return llhvals
end

function logpdf!(n::ProductNode, x::AbstractMatrix{<:Real}, llhvals::AxisArray{U}) where {U<:Real}
    @inbounds Y = mapreduce(c -> llhvals[:, c.id], hcat, children(n))
    llhvals[:,n.id] = map(U, _logpdf(n, x, Y))
    return llhvals
end

# ################## #
# Leaf distributions #
# ################## #
"""
    logpdf(n::Leaf, x)

Compute the logpdf of a leaf node.
"""
function logpdf!(n::Leaf, x::AbstractVector{<:Real}, llhvals::AxisArray{U}) where {U<:Real}
    llhvals[n.id] = map(U, n(x))
    return llhvals
end

function logpdf!(n::Leaf, x::AbstractMatrix{<:Real}, llhvals::AxisArray{U}) where {U<:Real}
    llhvals[:,n.id] .= map(U, n(x))
    return llhvals
end

@inline (n::IndicatorNode)(x::AbstractVector{<:Real}) = n(x, n.value)
@inline (n::IndicatorNode)(x::AbstractVector{<:Real}, θ::T) where {T<:Real} = _logpdf(n, x, θ)

function _logpdf(n::IndicatorNode, x::AbstractVector{<:Real}, θ::T) where {T<:Real}
    return isnan(x[scope(n)]) ? 0.0 : (x[scope(n)] == θ ? 0.0 : -Inf)
end

@inline (n::UnivariateNode)(x::AbstractVector{<:Real}) = n(x, Distributions.params(n.dist)...)
@inline (n::UnivariateNode)(x::AbstractVector{<:Real}, θ...) = _logpdf(n, x, θ...)

function _logpdf(n::UnivariateNode, x::AbstractVector{<:Real}, θ...)
    if isnan(x[scope(n)])
        return 0.0
    end

    if θ == params(n)
        return logpdf(n.dist, x[scope(n)])
    else
        return logpdf((typeof(n.dist))(θ...), x[scope(n)])
    end
end

@inline (n::MultivariateNode)(x::AbstractVector{<:Real}) = n(x, Distributions.params(n.dist)...)
@inline (n::MultivariateNode)(x::AbstractVector{<:Real}, θ...) = _logpdf(n, x, θ...)

function _logpdf(n::MultivariateNode, x::AbstractVector{<:Real}, θ...)

    @assert !any(isnan, x[scope(n)]) "NaN values for multivariate leaves is not supported."

    if θ == params(n)
        return logpdf(n.dist, x[scope(n)])
    else
        return logpdf((typeof(n.dist))(θ...), x[scope(n)])
    end
end

(n::IndicatorNode)(x::AbstractMatrix{<:Real}) = @inbounds map(i -> n(view(x,i,:)), 1:size(x,1))
(n::IndicatorNode)(x::AbstractMatrix{<:Real}, θ...) = @inbounds map(i -> n(view(x,i,:), θ...), 1:size(x,1))
(n::UnivariateNode)(x::AbstractMatrix{<:Real}) = @inbounds map(i -> n(view(x,i,:)), 1:size(x,1))
(n::UnivariateNode)(x::AbstractMatrix{<:Real}, θ...) = @inbounds map(i -> n(view(x,i,:), θ...), 1:size(x,1))
(n::MultivariateNode)(x::AbstractMatrix{<:Real}) = @inbounds map(i -> n(view(x,i,:)), 1:size(x,1))
(n::MultivariateNode)(x::AbstractMatrix{<:Real}, θ...) = @inbounds map(i -> n(view(x,i,:), θ...), 1:size(x,1))

@inline logpdf(n::Leaf, x::AbstractArray{<:Real}) = n(x)
@inline logpdf(n::Leaf, x::AbstractArray{<:Real}, y...) = n(x, y...)

rand(spn::SumProductNetwork) = rand(spn.root)

function rand(node::Node)
    @assert isnormalized(node)
    @assert hasscope(node)
    v = rand_(node)
    return map(d -> v[d], sort(scope(node)))
end

rand(node::IndicatorNode) = node.value
rand(node::UnivariateNode) = rand(node.dist)
rand(node::MultivariateNode) = rand(node.dist)

function rand_(node::ProductNode)
    @assert isnormalized(node)
    @assert hasscope(node)
    return mapreduce(c -> rand_(c), merge, filter(c -> hasscope(c), children(node)))
end

function rand_(node::SumNode)
    @assert isnormalized(node)
    @assert hasscope(node)
    w = Float64.(weights(node))
    z = rand(Categorical(w / sum(w))) # Normalisation due to precision errors.
    # Generate observation by drawing from a child.
    return rand_(children(node)[z])
end

rand_(node::IndicatorNode) = Dict(node.scope => node.value)
rand_(node::UnivariateNode) = Dict(node.scope => rand(node.dist))
function rand_(node::MultivariateNode)
	x = rand(node.dist)
	return Dict(d[2] => x[d[1]] for d in enumerate(node.scope))
end


text:
module Interface

using MacroTools: postwalk, isline, inexpr

function collectSymbols!(expression::Expr, symbols::Vector{Symbol})
    for arg in expression.args
        if hasproperty(arg, :args)
            collectSymbols!(arg, symbols)
        elseif !(typeof(arg) <: Number)
            push!(symbols, arg)
        end
    end
end

function isFunction(symbol::Symbol)
    if isdefined(Main, symbol)
        return isa(getfield(Main, symbol), Function)
    end
    return false
end

function replaceVars(expr::Union{Expr,Symbol}, mapping::Dict)
    if expr in keys(mapping)
        return :(u[$(mapping[expr]), :])
    end

    return expr
end

# Handle usage of Integers/Floats within expressions (such as exponents)
replaceVars(expr::Number, mapping::Dict) = return expr

struct ODEFunction <: Function
    func::Function
    dimensions::Int8
end

"""
    ODEFunc!(t, u, [constants...])

Update u based on current position and system defined using @ODE and provided
constants.

# Example
``` julia-repl
julia> lorenz = @ODE quote
    dx = σ*(y-x)
    dy = x*(ρ-z) - y
    dz = x * y - β*z
end

julia> u = [
    [1,2,3]
    [2,3,4]
    [4,5,6]
]

julia> lorenz(0.1, u, σ = 10, ρ = 28, β = 8/3)
```
"""
macro ODE(dx⃗::Expr...)

    dim = length(dx⃗)

    symbols = Vector{Symbol}()

    for expression in dx⃗
        collectSymbols!(expression, symbols)
    end

    # Filter out operators and functions
    filter!(!Meta.isoperator, symbols)
    filter!(!isFunction, symbols)
    filter!(!isequal(:t), symbols)

    diff_symbols = Vector{Symbol}([expr.args[1] for expr in dx⃗ if length(string(expr.args[1])) == 2])

    dependent_vars = Vector{Symbol}([Symbol(string(symbol)[end]) for symbol in diff_symbols])
    mapping_to_indices = Dict(symbol => index for (index, symbol) in enumerate(dependent_vars))

    constants = setdiff(symbols, diff_symbols, dependent_vars, [:t])
    constant_defaults = Dict{Symbol,Expr}([expr.args[1] => expr for expr in dx⃗ if (length(string(expr.args[1])) == 1) & (expr.args[1] in constants)])
    constants = [constant in keys(constant_defaults) ? Expr(:(kw), Expr(:(::), Symbol(constant), :(Real)), constant_defaults[constant].args[2]) : Expr(:(::), constant, :(Real)) for constant in constants]

    filter(expr -> expr.args[1] in keys(constant_defaults), dx⃗)

    diff_var_initialisation = [:($symbol = Vector{Float64}(undef, size(u, 2))) for symbol in diff_symbols]
    system = [var"@__dot__"(LineNumberNode(1), Main, expr).args[1] for expr in dx⃗ if !(expr in values(constant_defaults))] # Apply broadcasting
    system = [postwalk(x -> replaceVars(x, mapping_to_indices), expr) for expr in system] # Replace variables with indices

    if any(inexpr(expr, :t) for expr in dx⃗)
        func = quote
            function ODEFunc(t::Real, u::Matrix{<:Real}; $(constants...))::Matrix{Float64}
                $(diff_var_initialisation...)
                $(system...)
                return hcat($(diff_symbols...))
            end
        end
    else
        func = quote
            function ODEFunc(u::Matrix{<:Real}; $(constants...))::Matrix{Float64}
                $(diff_var_initialisation...)
                $(system...)
                return hcat($(diff_symbols...))
            end
        end
    end

    return func
end

# Method handles quote blocks being used to define ODEs
macro ODE(derivative_block::Expr)
    derivatives = filter(x -> !isline(x), derivative_block.args)
    @eval @ODE($(derivatives...))
end

end