text:
"""
    AbstractSphere{𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType}

An abstract type to represent a unit sphere that is represented isometrically in the embedding.
"""
abstract type AbstractSphere{𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType} end

@doc raw"""
    Sphere{n,𝔽} <: AbstractSphere{𝔽}

The (unit) sphere manifold $𝕊^{n}$ is the set of all unit norm vectors in $𝔽^{n+1}$.
The sphere is represented in the embedding, i.e.

````math
𝕊^{n} := \bigl\{ p \in 𝔽^{n+1}\ \big|\ \lVert p \rVert = 1 \bigr\}
````

where $𝔽\in\{ℝ,ℂ,ℍ\}$. Note that compared to the [`ArraySphere`](@ref), here the
argument `n` of the manifold is the dimension of the manifold, i.e. $𝕊^{n} ⊂ 𝔽^{n+1}$, $n\in ℕ$.

The tangent space at point $p$ is given by

````math
T_p𝕊^{n} := \bigl\{ X ∈ 𝔽^{n+1}\ |\ \Re(⟨p,X⟩) = 0 \bigr \},
````

where $𝔽\in\{ℝ,ℂ,ℍ\}$ and $⟨\cdot,\cdot⟩$ denotes the inner product in the
embedding $𝔽^{n+1}$.

For $𝔽=ℂ$, the manifold is the complex sphere, written $ℂ𝕊^n$, embedded in $ℂ^{n+1}$.
$ℂ𝕊^n$ is the complexification of the real sphere $𝕊^{2n+1}$.
Likewise, the quaternionic sphere $ℍ𝕊^n$ is the quaternionification of the real sphere
$𝕊^{4n+3}$.
Consequently, $ℂ𝕊^0$ is equivalent to $𝕊^1$ and [`Circle`](@ref), while $ℂ𝕊^1$ and $ℍ𝕊^0$
are equivalent to $𝕊^3$, though with different default representations.

This manifold is modeled as a special case of the more general case, i.e. as an embedded
manifold to the [`Euclidean`](@ref), and several functions like the [`inner`](@ref inner(::Euclidean, ::Any...)) product
and the [`zero_vector`](@ref zero_vector(::Euclidean, ::Any...)) are inherited from the embedding.

# Constructor

    Sphere(n[, field=ℝ])

Generate the (real-valued) sphere $𝕊^{n} ⊂ ℝ^{n+1}$, where `field` can also be used to
generate the complex- and quaternionic-valued sphere.
"""
struct Sphere{N,𝔽} <: AbstractSphere{𝔽} end
Sphere(n::Int, field::AbstractNumbers=ℝ) = Sphere{n,field}()

@doc raw"""
    ArraySphere{T<:Tuple,𝔽} <: AbstractSphere{𝔽}

The (unit) sphere manifold $𝕊^{n₁,n₂,...,nᵢ}$ is the set of all unit (Frobenius) norm elements of
$𝔽^{n₁,n₂,...,nᵢ}$, where $𝔽\in\{ℝ,ℂ,ℍ\}. The generalized sphere is
represented in the embedding, and supports arbitrary sized arrays or in other words arbitrary
tensors of unit norm. The set formally reads

````math
𝕊^{n_1, n_2, …, n_i} := \bigl\{ p \in 𝔽^{n_1, n_2, …, n_i}\ \big|\ \lVert p \rVert = 1 \bigr\}
````

where $𝔽\in\{ℝ,ℂ,ℍ\}$. Setting $i=1$ and $𝔽=ℝ$  this  simplifies to unit vectors in $ℝ^n$, see
[`Sphere`](@ref) for this special case. Note that compared to this classical case,
the argument for the generalized case here is given by the dimension of the embedding.
This means that `Sphere(2)` and `ArraySphere(3)` are the same manifold.

The tangent space at point $p$ is given by

````math
T_p 𝕊^{n_1, n_2, …, n_i} := \bigl\{ X ∈ 𝔽^{n_1, n_2, …, n_i}\ |\ \Re(⟨p,X⟩) = 0 \bigr \},
````

where $𝔽\in\{ℝ,ℂ,ℍ\}$ and $⟨\cdot,\cdot⟩$ denotes the (Frobenius) inner product in the
embedding $𝔽^{n_1, n_2, …, n_i}$.

This manifold is modeled as an embedded manifold to the [`Euclidean`](@ref), i.e.
several functions like the [`inner`](@ref inner(::Euclidean, ::Any...)) product and the
[`zero_vector`](@ref zero_vector(::Euclidean, ::Any...)) are inherited from the embedding.

# Constructor

    ArraySphere(n₁,n₂,...,nᵢ; field=ℝ)

Generate sphere in $𝔽^{n_1, n_2, …, n_i}$, where $𝔽$ defaults to the real-valued case $ℝ$.
"""
struct ArraySphere{N,𝔽} <: AbstractSphere{𝔽} where {N<:Tuple} end
function ArraySphere(n::Vararg{Int,I}; field::AbstractNumbers=ℝ) where {I}
    return ArraySphere{Tuple{n...},field}()
end

"""
    check_point(M::AbstractSphere, p; kwargs...)

Check whether `p` is a valid point on the [`AbstractSphere`](@ref) `M`, i.e. is a point in
the embedding of unit length.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_point(M::AbstractSphere, p; kwargs...)
    mpv = invoke(
        check_point,
        Tuple{(typeof(get_embedding(M))),typeof(p)},
        get_embedding(M),
        p;
        kwargs...,
    )
    mpv === nothing || return mpv
    if !isapprox(norm(p), 1.0; kwargs...)
        return DomainError(
            norm(p),
            "The point $(p) does not lie on the $(M) since its norm is not 1.",
        )
    end
    return nothing
end

"""
    check_vector(M::AbstractSphere, p, X; kwargs... )

Check whether `X` is a tangent vector to `p` on the [`AbstractSphere`](@ref) `M`, i.e.
after [`check_point`](@ref)`(M,p)`, `X` has to be of same dimension as `p`
and orthogonal to `p`.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_vector(M::AbstractSphere, p, X; kwargs...)
    mpv = invoke(
        check_vector,
        Tuple{typeof(get_embedding(M)),typeof(p),typeof(X)},
        get_embedding(M),
        p,
        X;
        kwargs...,
    )
    mpv === nothing || return mpv
    if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...)
        return DomainError(
            abs(dot(p, X)),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
        )
    end
    return nothing
end

function decorated_manifold(M::AbstractSphere{𝔽}) where {𝔽}
    return Euclidean(representation_size(M)...; field=𝔽)
end

# Since on every tangent space the Euclidean matric (restricted to this space) is used, this should be fine
default_metric_dispatch(::AbstractSphere, ::EuclideanMetric) = Val(true)

@doc raw"""
    distance(M::AbstractSphere, p, q)

Compute the geodesic distance betweeen `p` and `q` on the [`AbstractSphere`](@ref) `M`.
The formula is given by the (shorter) great arc length on the (or a) great circle
both `p` and `q` lie on.

````math
d_{𝕊}(p,q) = \arccos(\Re(⟨p,q⟩)).
````
"""
distance(::AbstractSphere, p, q) = acos(clamp(real(dot(p, q)), -1, 1))

@doc raw"""
    exp(M::AbstractSphere, p, X)

Compute the exponential map from `p` in the tangent direction `X` on the [`AbstractSphere`](@ref)
`M` by following the great arc eminating from `p` in direction `X`.

````math
\exp_p X = \cos(\lVert X \rVert_p)p + \sin(\lVert X \rVert_p)\frac{X}{\lVert X \rVert_p},
````
where $\lVert X \rVert_p$ is the [`norm`](@ref norm(::AbstractSphere,p,X)) on the
tangent space at `p` of the [`AbstractSphere`](@ref) `M`.
"""
exp(::AbstractSphere, ::Any...)

function exp!(M::AbstractSphere, q, p, X)
    θ = norm(M, p, X)
    q .= cos(θ) .* p .+ usinc(θ) .* X
    return q
end

function get_basis(M::Sphere{n,ℝ}, p, B::DiagonalizingOrthonormalBasis{ℝ}) where {n}
    A = zeros(n + 1, n + 1)
    A[1, :] = transpose(p)
    A[2, :] = transpose(B.frame_direction)
    V = nullspace(A)
    κ = ones(n)
    if !iszero(B.frame_direction)
        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
        V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V)
        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
    end
    T = typeof(similar(B.frame_direction))
    Ξ = [convert(T, V[:, i]) for i in 1:n]
    return CachedBasis(B, κ, Ξ)
end

@doc raw"""
    get_coordinates(M::AbstractSphere{ℝ}, p, X, B::DefaultOrthonormalBasis)

Represent the tangent vector `X` at point `p` from the [`AbstractSphere`](@ref) `M` in
an orthonormal basis by rotating the hyperplane containing `X` to a hyperplane whose
normal is the $x$-axis.

Given $q = p λ + x$, where $λ = \operatorname{sgn}(⟨x, p⟩)$, and $⟨⋅, ⋅⟩_{\mathrm{F}}$
denotes the Frobenius inner product, the formula for $Y$ is
````math
\begin{pmatrix}0 \\ Y\end{pmatrix} = X - q\frac{2 ⟨q, X⟩_{\mathrm{F}}}{⟨q, q⟩_{\mathrm{F}}}.
````
"""
get_coordinates(::AbstractSphere{ℝ}, p, X, ::DefaultOrthonormalBasis)

function get_coordinates!(
    M::AbstractSphere{ℝ},
    Y,
    p,
    X,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
)
    n = manifold_dimension(M)
    p1 = p[1]
    cosθ = abs(p1)
    λ = nzsign(p1, cosθ)
    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
    factor = λ * X[1] / (1 + cosθ)
    Y .= Xend .- pend .* factor
    return Y
end

get_embedding(M::AbstractSphere{𝔽}) where {𝔽} = decorated_manifold(M)

@doc raw"""
    get_vector(M::AbstractSphere{ℝ}, p, X, B::DefaultOrthonormalBasis)

Convert a one-dimensional vector of coefficients `X` in the basis `B` of the tangent space
at `p` on the [`AbstractSphere`](@ref) `M` to a tangent vector `Y` at `p` by rotating the
hyperplane containing `X`, whose normal is the $x$-axis, to the hyperplane whose normal is
`p`.

Given $q = p λ + x$, where $λ = \operatorname{sgn}(⟨x, p⟩)$, and $⟨⋅, ⋅⟩_{\mathrm{F}}$
denotes the Frobenius inner product, the formula for $Y$ is
````math
Y = X - q\frac{2 \left\langle q, \begin{pmatrix}0 \\ X\end{pmatrix}\right\rangle_{\mathrm{F}}}{⟨q, q⟩_{\mathrm{F}}}.
````
"""
get_vector(::AbstractSphere{ℝ}, p, X, ::DefaultOrthonormalBasis)

function get_vector!(
    M::AbstractSphere{ℝ},
    Y,
    p,
    X,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
)
    n = manifold_dimension(M)
    p1 = p[1]
    cosθ = abs(p1)
    λ = nzsign(p1, cosθ)
    pend = view(p, 2:(n + 1))
    pX = dot(pend, X)
    factor = pX / (1 + cosθ)
    Y[1] = -λ * pX
    Y[2:(n + 1)] .= X .- pend .* factor
    return Y
end

@doc raw"""
    injectivity_radius(M::AbstractSphere[, p])

Return the injectivity radius for the [`AbstractSphere`](@ref) `M`, which is globally $π$.

    injectivity_radius(M::Sphere, x, ::ProjectionRetraction)

Return the injectivity radius for the [`ProjectionRetraction`](@ref) on the
[`AbstractSphere`](@ref), which is globally $\frac{π}{2}$.
"""
injectivity_radius(::AbstractSphere) = π
injectivity_radius(::AbstractSphere, ::ExponentialRetraction) = π
injectivity_radius(::AbstractSphere, ::ProjectionRetraction) = π / 2
injectivity_radius(::AbstractSphere, ::Any) = π
injectivity_radius(::AbstractSphere, ::Any, ::ExponentialRetraction) = π
injectivity_radius(::AbstractSphere, ::Any, ::ProjectionRetraction) = π / 2
eval(
    quote
        @invoke_maker 1 AbstractManifold injectivity_radius(
            M::AbstractSphere,
            rm::AbstractRetractionMethod,
        )
    end,
)

@doc raw"""
    inverse_retract(M::AbstractSphere, p, q, ::ProjectionInverseRetraction)

Compute the inverse of the projection based retraction on the [`AbstractSphere`](@ref) `M`,
i.e. rearranging $p+X = q\lVert p+X\rVert_2$ yields
since $\Re(⟨p,X⟩) = 0$ and when $d_{𝕊^2}(p,q) ≤ \frac{π}{2}$ that

````math
\operatorname{retr}_p^{-1}(q) = \frac{q}{\Re(⟨p, q⟩)} - p.
````
"""
inverse_retract(::AbstractSphere, ::Any, ::Any, ::ProjectionInverseRetraction)

function inverse_retract!(::AbstractSphere, X, p, q, ::ProjectionInverseRetraction)
    return (X .= q ./ real(dot(p, q)) .- p)
end

@doc raw"""
    local_metric(M::Sphere{n}, p, ::DefaultOrthonormalBasis)

return the local representation of the metric in a [`DefaultOrthonormalBasis`](@ref), namely
the diagonal matrix of size ``n×n`` with ones on the diagonal, since the metric is obtained
from the embedding by restriction to the tangent space ``T_p\mathcal M`` at ``p``.
"""
function local_metric(::Sphere{n,ℝ}, p, B::DefaultOrthonormalBasis) where {n}
    return Diagonal(ones(SVector{n,eltype(p)}))
end

@doc raw"""
    log(M::AbstractSphere, p, q)

Compute the logarithmic map on the [`AbstractSphere`](@ref) `M`, i.e. the tangent vector,
whose geodesic starting from `p` reaches `q` after time 1.
The formula reads for $x ≠ -y$

````math
\log_p q = d_{𝕊}(p,q) \frac{q-\Re(⟨p,q⟩) p}{\lVert q-\Re(⟨p,q⟩) p \rVert_2},
````

and a deterministic choice from the set of tangent vectors is returned if $x=-y$, i.e. for
opposite points.
"""
log(::AbstractSphere, ::Any...)

function log!(M::AbstractSphere, X, p, q)
    cosθ = clamp(real(dot(p, q)), -1, 1)
    if cosθ ≈ -1 # appr. opposing points, return deterministic choice from set-valued log
        fill!(X, zero(eltype(X)))
        if p[1] ≈ 1
            X[2] = 1
        else
            X[1] = 1
        end
        copyto!(X, X .- real(dot(p, X)) .* p)
        X .*= π / norm(X)
    else
        θ = acos(cosθ)
        X .= (q .- cosθ .* p) ./ usinc(θ)
    end
    return project!(M, X, p, X)
end

@doc raw"""
    manifold_dimension(M::AbstractSphere)

Return the dimension of the [`AbstractSphere`](@ref) `M`, respectively i.e. the
dimension of the embedding -1.
"""
manifold_dimension(M::AbstractSphere) = manifold_dimension(get_embedding(M)) - 1

"""
    mean(
        S::AbstractSphere,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolationWithinRadius(π/2);
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using
[`GeodesicInterpolationWithinRadius`](@ref).
"""
mean(::AbstractSphere, ::Any...)

function Statistics.mean!(
    S::AbstractSphere,
    p,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)
    return mean!(S, p, x, w, GeodesicInterpolationWithinRadius(π / 2); kwargs...)
end

function mid_point!(S::Sphere, q, p1, p2)
    q .= p1 .+ p2
    project!(S, q, q)
    return q
end

@doc raw"""
    project(M::AbstractSphere, p)

Project the point `p` from the embedding onto the [`Sphere`](@ref) `M`.

````math
\operatorname{proj}(p) = \frac{p}{\lVert p \rVert},
````
where $\lVert\cdot\rVert$ denotes the usual 2-norm for vectors if $m=1$ and the Frobenius
norm for the case $m>1$.
"""
project(::AbstractSphere, ::Any)

project!(::AbstractSphere, q, p) = (q .= p ./ norm(p))

@doc raw"""
    project(M::AbstractSphere, p, X)

Project the point `X` onto the tangent space at `p` on the [`Sphere`](@ref) `M`.

````math
\operatorname{proj}_{p}(X) = X - \Re(⟨p, X⟩)p
````
"""
project(::AbstractSphere, ::Any, ::Any)

project!(::AbstractSphere, Y, p, X) = (Y .= X .- real(dot(p, X)) .* p)

@doc raw"""
    representation_size(M::AbstractSphere)

Return the size points on the [`AbstractSphere`](@ref) `M` are represented as, i.e., the
representation size of the embedding.
"""
@generated representation_size(M::ArraySphere{N}) where {N} = size_to_tuple(N)
@generated representation_size(M::Sphere{N}) where {N} = (N + 1,)

@doc raw"""
    retract(M::AbstractSphere, p, X, ::ProjectionRetraction)

Compute the retraction that is based on projection, i.e.

````math
\operatorname{retr}_p(X) = \frac{p+X}{\lVert p+X \rVert_2}
````
"""
retract(::AbstractSphere, ::Any, ::Any, ::ProjectionRetraction)

function retract!(M::AbstractSphere, q, p, X, ::ProjectionRetraction)
    q .= p .+ X
    return project!(M, q, q)
end

Base.show(io::IO, ::Sphere{n,𝔽}) where {n,𝔽} = print(io, "Sphere($(n), $(𝔽))")
function Base.show(io::IO, ::ArraySphere{N,𝔽}) where {N,𝔽}
    return print(io, "ArraySphere($(join(N.parameters, ", ")); field = $(𝔽))")
end

"""
    uniform_distribution(M::Sphere{n,ℝ}, p) where {n}

Uniform distribution on given [`Sphere`](@ref) `M`. Generated points will be of
similar type as `p`.
"""
function uniform_distribution(M::Sphere{n,ℝ}, p) where {n}
    d = Distributions.MvNormal(zero(p), 1.0)
    return ProjectedPointDistribution(M, d, project!, p)
end

@doc raw"""
    vector_transport_to(M::AbstractSphere, p, X, q, ::ParallelTransport)

Compute the parallel transport on the [`Sphere`](@ref) of the tangent vector `X` at `p`
to `q`, provided, the [`geodesic`](@ref) between `p` and `q` is unique. The formula reads

````math
P_{p←q}(X) = X - \frac{\Re(⟨\log_p q,X⟩_p)}{d^2_𝕊(p,q)}
\bigl(\log_p q + \log_q p \bigr).
````
"""
vector_transport_to(::AbstractSphere, ::Any, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(::AbstractSphere, Y, p, X, q, ::ParallelTransport)
    m = p .+ q
    mnorm2 = real(dot(m, m))
    factor = 2 * real(dot(X, q)) / mnorm2
    Y .= X .- m .* factor
    return Y
end

"""
    StereographicAtlas()

The stereographic atlas of ``S^n`` with two charts: one with the singular
point (-1, 0, ..., 0) (called `:north`) and one with the singular
point (1, 0, ..., 0) (called `:south`).
"""
struct StereographicAtlas <: AbstractAtlas{ℝ} end

function get_chart_index(::Sphere{n,ℝ}, ::StereographicAtlas, p) where {n}
    if p[1] < 0
        return :south
    else
        return :north
    end
end

function get_parameters!(::Sphere{n,ℝ}, x, ::StereographicAtlas, i::Symbol, p) where {n}
    if i === :north
        return x .= p[2:end] ./ (1 + p[1])
    else
        return x .= p[2:end] ./ (1 - p[1])
    end
end

function get_point!(::Sphere{n,ℝ}, p, ::StereographicAtlas, i::Symbol, x) where {n}
    xnorm2 = dot(x, x)
    if i === :north
        p[1] = (1 - xnorm2) / (xnorm2 + 1)
    else
        p[1] = (xnorm2 - 1) / (xnorm2 + 1)
    end
    p[2:end] .= 2 * x / (xnorm2 + 1)
    return p
end

function get_coordinates!(
    ::Sphere{n,ℝ},
    Y,
    p,
    X,
    B::InducedBasis{ℝ,TangentSpaceType,<:StereographicAtlas},
) where {n}
    if B.i === :north
        for i in 1:n
            Y[i] = X[i + 1] / (1 + p[1]) - X[1] * p[i + 1] / (1 + p[1])^2
        end
    else
        for i in 1:n
            Y[i] = X[i + 1] / (-1 + p[1]) - X[1] * p[i + 1] / (-1 + p[1])^2
        end
    end
    return Y
end

function get_vector!(
    M::Sphere{n,ℝ},
    Y,
    p,
    X,
    B::InducedBasis{ℝ,TangentSpaceType,<:StereographicAtlas},
) where {n}
    a = get_parameters(M, B.A, B.i, p)
    mult = inv(1 + dot(a, a))^2

    Y[1] = 0
    for j in 1:n
        Y[1] -= 4 * a[j] * mult * X[j]
    end
    for i in 2:(n + 1)
        Y[i] = 0
        for j in 1:n
            if i == j + 1
                Y[i] += 2 * (1 + dot(a, a) - 2 * a[i - 1]^2) * mult * X[j]
            else
                Y[i] -= 4 * a[i - 1] * a[j] * mult * X[j]
            end
        end
        if B.i === :south
            Y[i] *= -1
        end
    end
    return Y
end

function local_metric(
    M::Sphere{n,ℝ},
    p,
    B::InducedBasis{ℝ,TangentSpaceType,StereographicAtlas,Symbol},
) where {n}
    a = get_parameters(M, B.A, B.i, p)
    return (4 / (1 + dot(a, a))^2) * I
end


text:
using EchoviewEvr

using Test

filename = joinpath(dirname(@__FILE__),
                  "data/WCBjr280_1_1_regions.evr")

_regions = collect(regions(filename))

@test length(_regions) == 68
@test _regions[1].classification == "krills"
@test _regions[1].regiontype == "1"

_polygons = polygons(_regions)
@test length(_polygons) == 68

_polygons = polygons(filename)
@test length(_polygons) == 68

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the constructors of the quaternions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/quaternion.jl
# =========================

# Function: Quaternion
# --------------------

@testset "Constructors of the quaternions" begin
    # Constructors with all the components as arguments.
    q = Quaternion(1, 0, 0, 0)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Int

    q = Quaternion(1, 0, 0, 0.0)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === promote_type(Int, typeof(0.0))

    # Constructors with one vector.
    q = Quaternion([0, 0, 0])
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Int

    v = [0, cosd(60), sind(60)]
    q = Quaternion(v)
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == cosd(60)
    @test q.q3 == sind(60)
    @test eltype(q) === eltype(v)

    v = [cosd(60), 0, sind(60), 0]
    q = Quaternion(v)
    @test q.q0 == cosd(60)
    @test q.q1 == 0
    @test q.q2 == sind(60)
    @test q.q3 == 0
    @test eltype(q) === eltype(v)

    # Constructores with one element and one vector.
    r = cosd(60)
    v = [0, sind(60), 0]
    q = Quaternion(r, v)
    @test q.q0 == cosd(60)
    @test q.q1 == 0
    @test q.q2 == sind(60)
    @test q.q3 == 0
    @test eltype(q) === promote_type(eltype(r), eltype(v))

    r = 0
    v = [0, sind(60), cosd(60)]
    q = Quaternion(r, v)
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == sind(60)
    @test q.q3 == cosd(60)
    @test eltype(q) === promote_type(eltype(r), eltype(v))

    r = 0
    v = [0, 0, 1]
    q = Quaternion(r, v)
    @test q.q0 == 0
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 1
    @test eltype(q) === promote_type(eltype(r), eltype(v))

    # Constructors with the uniform scaling.
    q = Quaternion(I)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === eltype(I)

    r = 1.0
    q = Quaternion(r * I)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === eltype(r)

    q = Quaternion{Float32}(I)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === Float32

    q1 = Quaternion{Float64}(I)
    q = Quaternion(I, q1)
    @test q.q0 == 1
    @test q.q1 == 0
    @test q.q2 == 0
    @test q.q3 == 0
    @test eltype(q) === eltype(q1)
end

@testset "Constructors of the quaternions (errors)" begin
    @test_throws ArgumentError Quaternion([1, 1])
    @test_throws ArgumentError Quaternion([1, 1, 2, 3, 4])
end


text:
module GL2017Replication

export hello, domath
export compute_steady_state!,initilize!,aggregate! ,EGM!, compute_distribution!
export calibrate!, calibrate_terminal,load_parameters_iss!
export transition!, EGM_trans!
export ModelGL,TransGL
export describe, print_params,print_status
export plots_figure_1,plots_figure_3,plots_figure_4
# Import scripts
include("model_structures.jl")
include("model_functions.jl")
include("output_functions.jl")

end


text:
using LatSpec
using Test

@testset verbose=true "LatSpec.jl" begin
    include("datapoint.jl")
    include("statistics.jl")
end


text:
using LiteQTL
using CSV
using DelimitedFiles
using DataFrames

pheno_file = "/home/xiaoqihu/git/LiteQTL-G3-supplement/data/tensorqtldata/cleanpheno.csv"
geno_file = "/home/xiaoqihu/git/LiteQTL-G3-supplement/data/tensorqtldata/genotype.csv"
export_matrix = false

@time pheno = LiteQTL.get_pheno_data(pheno_file, Float64, transposed=true)
# 10.865656 seconds (47.21 M allocations: 1.625 GiB, 12.64% gc time)
@time geno = CSV.read(geno_file)
# 109.905806 seconds (41.48 M allocations: 46.943 GiB, 0.58% gc time)
# @time genomat = convert(Matrix{Float32}, geno[:, 3:end]) |> transpose |> collect
# geno = Nothing

# pick out chr1 genos
genobychr = groupby(geno, :Column2)
onechr = genobychr[(Column2="9",)]
# Chr1: 5618644x447;
# Chr2: 1738829;
# Chr3: 919842;
# Chr4: 944044;
# Chr5: 820554;
# Chr6: 849929
# Chr7: 762694
# Chr8: 719376
# Chr9: 558963
# ChrX: 436393

# It takes about 60 seconds to do this step.
datatype = Float64
@time genomat = convert(Matrix{datatype}, onechr[:, 3:end]) |> transpose |> collect
pheno = convert(Matrix{datatype}, pheno)
# gc
genobychr = Nothing;
geno = Nothing;
chr1geno = Nothing;

LiteQTL.set_blas_threads(16);

G = genomat;
Y = pheno;

genomat = Nothing
pheno = Nothing

n = size(Y,1)
m = size(Y,2)
p = size(G,2)
println("******* Indivuduals n: $n, Traits m: $m, Markers p: $p ****************");

@time lodc = scan(Y, G, export_matrix=export_matrix, usegpu=false, lod_or_pval="lod");

@time lodc = LiteQTL.scan(Y, G,n;export_matrix, "lod", false);
@time lodg = LiteQTL.scan(Y, G,n;usegpu=true)


text:
import GLFW

window = GLFW.CreateWindow(800, 600, "GLFW Callback Test")
GLFW.MakeContextCurrent(window)

# Window callbacks
GLFW.SetWindowCloseCallback(window, (_) -> println("window closed"))
GLFW.SetWindowPosCallback(window, (_, x, y) -> println("window position: $x, $y"))
GLFW.SetWindowSizeCallback(window, (_, w, h) -> println("window size: $h x $w"))
GLFW.SetWindowFocusCallback(window, (_, focused) -> println("window focus: $focused"))
GLFW.SetWindowIconifyCallback(window, (_, iconified) -> println("window iconify: $iconified"))

# Input callbacks
actiontext = Dict(
	GLFW.RELEASE => "released",
	GLFW.PRESS   => "pressed",
	GLFW.REPEAT  => "repeats"
)
GLFW.SetKeyCallback(window, (_, key, scancode, action, mods) -> begin
	if '!' <= key && key <= '~'
		key = convert(Char, key)
		key = string("'", key, "'")
	end
	println("key $key ", actiontext[action])
end)

GLFW.SetCharModsCallback(window, (_, c, mods) -> println("char: $c, mods: $mods"))
GLFW.SetMouseButtonCallback(window, (_, button, action, mods) ->
	println("mouse button $button ", actiontext[action]))
GLFW.SetCursorPosCallback(window, (_, x, y) -> println("cursor: $x, $y"))
GLFW.SetScrollCallback(window, (_, xoff, yoff) -> println("scroll: $xoff, $yoff"))

GLFW.SetDropCallback(window, (_, paths) -> println(paths))

glClear() = ccall(@eval(GLFW.GetProcAddress("glClear")), Void, (Cuint,), 0x00004000)

gc() # Force garbage collection so that improper reference management is more apparent via crashes

while !GLFW.WindowShouldClose(window)
	glClear()
	GLFW.SwapBuffers(window)
	GLFW.WaitEvents()
end

GLFW.DestroyWindow(window)


text:
using YaoPlots, Yao
using Test, Compose

@testset "LabelBlock" begin
    x = put(5, (2,3)=>matblock(rand_unitary(4)))
    cb = LabelBlock(x, "x")
    @test mat(copy(cb)) == mat(cb)
    @test isunitary(cb)
    @test ishermitian(cb) == ishermitian(x)
    @test isreflexive(cb) == isreflexive(x)
    @test mat(cb) == mat(x)
    reg = rand_state(5)
    @test apply!(copy(reg), cb) ≈ apply!(copy(reg), x)
    @test cb' isa LabelBlock && mat(cb') ≈ mat(cb)'
    @test (cb').name == "x†" && (cb'').name == "x"

    y = put(5, (3,4)=>matblock(rand_unitary(4)))
    cc = chsubblocks(cb, y)
    @test YaoPlots.is_continuous_chunk([1,2,3]) == true
    @test YaoPlots.is_continuous_chunk([1,2,4]) == false
    @test YaoPlots.is_continuous_chunk([3,2,4]) == true

    c1 = chain(5, [put(5, (2,3)=>label(SWAP, "SWAP")), put(5, 2=>label(I2, "id")), put(5, 2=>label(X, "X")), control(5, (5,3), (2,4,1)=>put(3, (1,3)=>label(SWAP, "SWAP")))])
    c2 = chain(5, [put(5, (2,3)=>label(SWAP, "SWAP")), put(5, 2=>label(I2, "id")), put(5, 2=>label(X, "X")), control(5, (5,3), (2,4,1)=>put(3, (1,2)=>label(SWAP, "SWAP")))])

    @test vizcircuit(c1) isa Compose.Context
    @test_throws ErrorException vizcircuit(c2)
end


text:
using Test, YAAD, YAAD.TestUtils

@testset "cat" begin

@test gradcheck((x, y)->hcat(x, y), Variable(rand(2, 2)), Variable(rand(2)))
@test gradcheck((x, y)->cat(x, y; dims=2), Variable(rand(2, 2)), Variable(rand(2)))

end


text:
module RandomFloat

import Base: frexp, ldexp

export randfloat

# frexp([x,y]) == ([frexp(x)[1], frexp(y)[1]], [frexp(x)[2], frexp(y)[2]]])
# ldexp(frexp([x,y])...) == ldexp(frexp([x,y])[1], frexp([x,y])[2]) == [x,y]
@vectorize_1arg Real frexp
@vectorize_2arg Real ldexp

for (F,I) in ((:Float16, :Int16), (:Float32, :Int32), (:Float64, Int64))
  @eval begin
    function randfloat(r::FloatRange{$F}=sqrt(eps(one($F))):(one($F)/sqrt(eps(one($F)))), dims...)
      firstIsNeg = signbit(first(r))
      lastIsNeg  = signbit(last(r))
      if !firstIsNeg
        lo = reinterpret($I,first(r))
        hi = reinterpret($I,last(r))
        reinterpret($F, rand(lo:hi, dims...)
        #step = reinterpret($I, eps(eps(last(r)))/2)
        #reinterpret($F, rand(lo:step:hi, dims...))
      elseif lastIsNeg
        -randfloat((-last(r)):(-first(r)), dims...)
      else
        throw(ArgumentError("expecting range to be negative or positive, not mixed"))
      end
    end
  end
end

end # module


text:
module PDBenchmark

# Write your package code here.
using PeriDyn
using PDMesh

include("./tests/test.jl")
include("./tests/standardtests.jl")
include("./geom/geom.jl")

end

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark


text:
@everywhere using ArgoData
@everywhere path=joinpath(tempdir(),"Argo_DAC_files")
@everywhere fil=joinpath(path,"Argo_float_files.csv")
#@everywhere ftp="ftp://usgodae.org/pub/outgoing/argo/dac/"

!isdir(path) ? mkdir(path) : nothing

if !isfile(fil)
    files_list=GDAC.files_list()
    path=joinpath(tempdir(),"Argo_DAC_files")
    !isdir(path) ? mkdir(path) : nothing
    GDAC.CSV.write(joinpath(path,"Argo_float_files.csv"),files_list)
end

@everywhere files_list=GDAC.DataFrame(GDAC.CSV.File(fil))

for i in unique(files_list[:,:folder])
    !isdir(joinpath(path,i)) ? mkdir(joinpath(path,i)) : nothing
end

@everywhere n=10
@everywhere N=0
while N<1 #size(files_list,1)
    @sync @distributed for m in 1:nworkers()
        ii=collect(N + (m-1)*n .+ (1:n))
        println(ii[1])
        #[Argo_float_download(files_list[i,:],"prof",ftp) for i in ii];
        [GDAC.download_file(files_list[i,:],"prof") for i in ii];
    end
    @everywhere N=N+n*nworkers()
end


text:
""" Free diagrams in a category.

A [free diagram](https://ncatlab.org/nlab/show/free+diagram) in a category is a
diagram whose shape is a free category. Examples include the empty diagram,
pairs of objects, discrete diagrams, parallel pairs, composable pairs, and spans
and cospans. Limits and colimits are most commonly taken over free diagrams.
"""
module FreeDiagrams
export FreeDiagram, BipartiteFreeDiagram, FixedShapeFreeDiagram,
  DiscreteDiagram, EmptyDiagram, SingletonDiagram, ObjectPair,
  Span, Cospan, Multispan, Multicospan, SMultispan, SMulticospan,
  ParallelPair, ParallelMorphisms, ComposablePair, ComposableMorphisms,
  diagram_type, cone_objects, cocone_objects,
  ob, hom, dom, codom, apex, legs, feet, left, right, bundle_legs,
  nv, ne, src, tgt, vertices, edges, has_vertex, has_edge,
  add_vertex!, add_vertices!, add_edge!, add_edges!,
  ob₁, ob₂, nv₁, nv₂, vertices₁, vertices₂,
  add_vertex₁!, add_vertex₂!, add_vertices₁!, add_vertices₂!

using AutoHashEquals
using StaticArrays: StaticVector, SVector

using ...Present, ...Theories, ...CSetDataStructures, ...Graphs, ..FinCats
import ...Theories: ob, hom, dom, codom, left, right
using ...Graphs.BasicGraphs: TheoryGraph
using ...Graphs.BipartiteGraphs: TheoryUndirectedBipartiteGraph
import ..FinCats: FreeCatGraph, FinDomFunctor, collect_ob, collect_hom

# Diagram interface
###################

""" Given a diagram in a category ``C``, return Julia type of objects and
morphisms in ``C`` as a tuple type of form ``Tuple{Ob,Hom}``.
"""
function diagram_type end

""" Objects in diagram that will have explicit legs in limit cone.

In category theory, it is common practice to elide legs of limit cones that can
be computed from other legs, especially for diagrams of certain fixed shapes.
For example, when it taking a pullback (the limit of a cospan), the limit object
is often treated as having two projections, rather than three. This function
encodes such conventions by listing the objects in the diagram that will have
corresponding legs in the limit object created by Catlab.

See also: [`cocone_objects`](@ref).
"""
cone_objects(diagram) = ob(diagram)

""" Objects in diagram that will have explicit legs in colimit cocone.

See also: [`cone_objects`](@ref).
"""
cocone_objects(diagram) = ob(diagram)

# Diagrams of fixed shape
#########################

""" Abstract type for free diagram of fixed shape.
"""
abstract type FixedShapeFreeDiagram{Ob,Hom} end

diagram_type(::FixedShapeFreeDiagram{Ob,Hom}) where {Ob,Hom} = Tuple{Ob,Hom}

# Discrete diagrams
#------------------

""" Discrete diagram: a diagram with no non-identity morphisms.
"""
@auto_hash_equals struct DiscreteDiagram{Ob,Hom,Obs<:AbstractVector{Ob}} <:
    FixedShapeFreeDiagram{Ob,Hom}
  objects::Obs
end

DiscreteDiagram(objects::Obs, Hom::Type=Any) where {Ob,Obs<:AbstractVector{Ob}} =
  DiscreteDiagram{Ob,Hom,Obs}(objects)

const EmptyDiagram{Ob,Hom} = DiscreteDiagram{Ob,Hom,<:StaticVector{0,Ob}}
const SingletonDiagram{Ob,Hom} = DiscreteDiagram{Ob,Hom,<:StaticVector{1,Ob}}
const ObjectPair{Ob,Hom} = DiscreteDiagram{Ob,Hom,<:StaticVector{2,Ob}}

EmptyDiagram{Ob}(Hom::Type=Any) where Ob = DiscreteDiagram(SVector{0,Ob}(), Hom)
SingletonDiagram(ob, Hom::Type=Any) = DiscreteDiagram(SVector(ob), Hom)
ObjectPair(first, second, Hom::Type=Any) =
  DiscreteDiagram(SVector(first, second), Hom)

ob(d::DiscreteDiagram) = d.objects

Base.iterate(d::DiscreteDiagram, args...) = iterate(d.objects, args...)
Base.eltype(d::DiscreteDiagram) = eltype(d.objects)
Base.length(d::DiscreteDiagram) = length(d.objects)
Base.getindex(d::DiscreteDiagram, i) = d.objects[i]
Base.firstindex(d::DiscreteDiagram) = firstindex(d.objects)
Base.lastindex(d::DiscreteDiagram) = lastindex(d.objects)

# Spans and cospans
#------------------

""" Multispan of morphisms in a category.

A [multispan](https://ncatlab.org/nlab/show/multispan) is like a [`Span`](@ref)
except that it may have a number of legs different than two. A colimit of this
shape is a pushout.
"""
@auto_hash_equals struct Multispan{Ob,Hom,Legs<:AbstractVector{Hom}} <:
    FixedShapeFreeDiagram{Ob,Hom}
  apex::Ob
  legs::Legs
end

function Multispan(legs::AbstractVector)
  !isempty(legs) || error("Empty list of legs but no apex given")
  allequal(dom.(legs)) || error("Legs $legs do not have common domain")
  Multispan(dom(first(legs)), legs)
end

const SMultispan{N,Ob,Hom} = Multispan{Ob,Hom,<:StaticVector{N,Hom}}

SMultispan{N}(apex, legs::Vararg{Any,N}) where N =
  Multispan(apex, SVector{N}(legs...))
SMultispan{0}(apex) = Multispan(apex, SVector{0,typeof(id(apex))}())
SMultispan{N}(legs::Vararg{Any,N}) where N = Multispan(SVector{N}(legs...))

""" Span of morphims in a category.

A common special case of [`Multispan`](@ref). See also [`Cospan`](@ref).
"""
const Span{Ob,Hom} = SMultispan{2,Ob,Hom}

""" Apex of multispan or multicospan.

The apex of a multi(co)span is the object that is the (co)domain of all the
[`legs`](@ref).
"""
apex(span::Multispan) = span.apex

""" Legs of multispan or multicospan.

The legs are the morphisms comprising the multi(co)span.
"""
legs(span::Multispan) = span.legs

""" Feet of multispan or multicospan.

The feet of a multispan are the codomains of the [`legs`](@ref).
"""
feet(span::Multispan) = map(codom, span.legs)

""" Left leg of span or cospan.
"""
left(span::Span) = span.legs[1]

""" Right leg of span or cospan.
"""
right(span::Span) = span.legs[2]

cocone_objects(span::Multispan) = feet(span)

Base.iterate(span::Multispan, args...) = iterate(span.legs, args...)
Base.eltype(span::Multispan) = eltype(span.legs)
Base.length(span::Multispan) = length(span.legs)

""" Multicospan of morphisms in a category.

A multicospan is like a [`Cospan`](@ref) except that it may have a number of
legs different than two. A limit of this shape is a pullback.
"""
@auto_hash_equals struct Multicospan{Ob,Hom,Legs<:AbstractVector{Hom}} <:
    FixedShapeFreeDiagram{Ob,Hom}
  apex::Ob
  legs::Legs
end

function Multicospan(legs::AbstractVector)
  !isempty(legs) || error("Empty list of legs but no apex given")
  allequal(codom.(legs)) || error("Legs $legs do not have common codomain")
  Multicospan(codom(first(legs)), legs)
end

const SMulticospan{N,Ob,Hom} = Multicospan{Ob,Hom,<:StaticVector{N,Hom}}

SMulticospan{N}(apex, legs::Vararg{Any,N}) where N =
  Multicospan(apex, SVector{N}(legs...))
SMulticospan{0}(apex) = Multicospan(apex, SVector{0,typeof(id(apex))}())
SMulticospan{N}(legs::Vararg{Any,N}) where N = Multicospan(SVector{N}(legs...))

""" Cospan of morphisms in a category.

A common special case of [`Multicospan`](@ref). See also [`Span`](@ref).
"""
const Cospan{Ob,Hom} = SMulticospan{2,Ob,Hom}

apex(cospan::Multicospan) = cospan.apex
legs(cospan::Multicospan) = cospan.legs
feet(cospan::Multicospan) = map(dom, cospan.legs)
left(cospan::Cospan) = cospan.legs[1]
right(cospan::Cospan) = cospan.legs[2]

cone_objects(cospan::Multicospan) = feet(cospan)

Base.iterate(cospan::Multicospan, args...) = iterate(cospan.legs, args...)
Base.eltype(cospan::Multicospan) = eltype(cospan.legs)
Base.length(cospan::Multicospan) = length(cospan.legs)

""" Bundle together legs of a multi(co)span.

For example, calling `bundle_legs(span, SVector((1,2),(3,4)))` on a multispan
with four legs gives a span whose left leg bundles legs 1 and 2 and whose right
leg bundles legs 3 and 4. Note that in addition to bundling, this function can
also permute legs and discard them.

The bundling is performed using the universal property of (co)products, which
assumes that these (co)limits exist.
"""
bundle_legs(span::Multispan, indices) =
  Multispan(apex(span), map(i -> bundle_leg(span, i), indices))
bundle_legs(cospan::Multicospan, indices) =
  Multicospan(apex(cospan), map(i -> bundle_leg(cospan, i), indices))

bundle_leg(x::Union{Multispan,Multicospan}, i::Int) = legs(x)[i]
bundle_leg(x::Union{Multispan,Multicospan}, i::Tuple) = bundle_leg(x, SVector(i))
bundle_leg(span::Multispan, i::AbstractVector{Int}) = pair(legs(span)[i])
bundle_leg(cospan::Multicospan, i::AbstractVector{Int}) = copair(legs(cospan)[i])

# Parallel morphisms
#-------------------

""" Parallel morphims in a category.

[Parallel morphisms](https://ncatlab.org/nlab/show/parallel+morphisms) are just
morphisms with the same domain and codomain. A (co)limit of this shape is a
(co)equalizer.

For the common special case of two morphisms, see [`ParallelPair`](@ref).
"""
@auto_hash_equals struct ParallelMorphisms{Dom,Codom,Hom,Homs<:AbstractVector{Hom}} <:
    FixedShapeFreeDiagram{Union{Dom,Codom},Hom}
  dom::Dom
  codom::Codom
  homs::Homs
end

function ParallelMorphisms(homs::AbstractVector)
  @assert !isempty(homs) && allequal(dom.(homs)) && allequal(codom.(homs))
  ParallelMorphisms(dom(first(homs)), codom(first(homs)), homs)
end

""" Pair of parallel morphisms in a category.

A common special case of [`ParallelMorphisms`](@ref).
"""
const ParallelPair{Dom,Codom,Hom} =
  ParallelMorphisms{Dom,Codom,Hom,<:StaticVector{2,Hom}}

function ParallelPair(first, last)
  dom(first) == dom(last) ||
    error("Domains of parallel pair do not match: $first vs $last")
  codom(first) == codom(last) ||
    error("Codomains of parallel pair do not match: $first vs $last")
  ParallelMorphisms(dom(first), codom(first), SVector(first, last))
end

dom(para::ParallelMorphisms) = para.dom
codom(para::ParallelMorphisms) = para.codom
hom(para::ParallelMorphisms) = para.homs

cone_objects(para::ParallelMorphisms) = SVector(dom(para))
cocone_objects(para::ParallelMorphisms) = SVector(codom(para))

Base.iterate(para::ParallelMorphisms, args...) = iterate(para.homs, args...)
Base.eltype(para::ParallelMorphisms) = eltype(para.homs)
Base.length(para::ParallelMorphisms) = length(para.homs)
Base.getindex(para::ParallelMorphisms, i) = para.homs[i]
Base.firstindex(para::ParallelMorphisms) = firstindex(para.homs)
Base.lastindex(para::ParallelMorphisms) = lastindex(para.homs)

allequal(xs::AbstractVector) = isempty(xs) || all(==(xs[1]), xs)

# Composable morphisms
#---------------------

""" Composable morphisms in a category.

Composable morphisms are a sequence of morphisms in a category that form a path
in the underlying graph of the category.

For the common special case of two morphisms, see [`ComposablePair`](@ref).
"""
@auto_hash_equals struct ComposableMorphisms{Ob,Hom,Homs<:AbstractVector{Hom}} <:
    FixedShapeFreeDiagram{Ob,Hom}
  homs::Homs
end

function ComposableMorphisms(homs::Homs) where {Hom, Homs<:AbstractVector{Hom}}
  doms, codoms = dom.(homs), codom.(homs)
  all(c == d for (c,d) in zip(codoms[1:end-1], doms[2:end])) ||
    error("Domain mismatch in composable sequence: $homs")
  ComposableMorphisms{Union{eltype(doms),eltype(codoms)},Hom,Homs}(homs)
end

""" Pair of composable morphisms in a category.

[Composable pairs](https://ncatlab.org/nlab/show/composable+pair) are a common
special case of [`ComposableMorphisms`](@ref).
"""
const ComposablePair{Ob,Hom} = ComposableMorphisms{Ob,Hom,<:StaticVector{2,Hom}}

ComposablePair(first, last) = ComposableMorphisms(SVector(first, last))

dom(comp::ComposableMorphisms) = dom(first(comp))
codom(comp::ComposableMorphisms) = codom(last(comp))
hom(comp::ComposableMorphisms) = comp.homs

Base.iterate(comp::ComposableMorphisms, args...) = iterate(comp.homs, args...)
Base.eltype(comp::ComposableMorphisms) = eltype(comp.homs)
Base.length(comp::ComposableMorphisms) = length(comp.homs)
Base.getindex(comp::ComposableMorphisms, i) = comp.homs[i]
Base.firstindex(comp::ComposableMorphisms) = firstindex(comp.homs)
Base.lastindex(comp::ComposableMorphisms) = lastindex(comp.homs)

# Diagrams of flexible shape
############################

# Bipartite free diagrams
#------------------------

@abstract_acset_type _BipartiteFreeDiagram <: AbstractUndirectedBipartiteGraph

""" A free diagram with a bipartite structure.

Such diagrams include most of the fixed shapes, such as spans, cospans, and
parallel morphisms. They are also the generic shape of diagrams for limits and
colimits arising from undirected wiring diagrams. For limits, the boxes
correspond to vertices in ``V₁`` and the junctions to vertices in ``V₂``.
Colimits are dual.
"""
const BipartiteFreeDiagram{Ob,Hom,S} = _BipartiteFreeDiagram{S,Tuple{Ob,Hom}}

@present TheoryBipartiteFreeDiagram <: TheoryUndirectedBipartiteGraph begin
  Ob::AttrType
  Hom::AttrType
  ob₁::Attr(V₁, Ob)
  ob₂::Attr(V₂, Ob)
  hom::Attr(E, Hom)
end

""" The default concrete type for bipartite free diagrams.
"""
@acset_type BasicBipartiteFreeDiagram(TheoryBipartiteFreeDiagram,
  index=[:src, :tgt]) <: _BipartiteFreeDiagram

@present TheoryFreeDiagramAsBipartite <: TheoryBipartiteFreeDiagram begin
  V::Ob
  orig_vert₁::Hom(V₁, V)
  orig_vert₂::Hom(V₂, V)
end

""" A free diagram that has been converted to a bipartite free diagram.
"""
@acset_type FreeDiagramAsBipartite(TheoryFreeDiagramAsBipartite,
  index=[:src, :tgt], unique_index=[:orig_vert₁, :orig_vert₂]) <: _BipartiteFreeDiagram

ob₁(d::BipartiteFreeDiagram, args...) = subpart(d, args..., :ob₁)
ob₂(d::BipartiteFreeDiagram, args...) = subpart(d, args..., :ob₂)
hom(d::BipartiteFreeDiagram, args...) = subpart(d, args..., :hom)

diagram_type(::BipartiteFreeDiagram{Ob,Hom}) where {Ob,Hom} = Tuple{Ob,Hom}
cone_objects(diagram::BipartiteFreeDiagram) = ob₁(diagram)
cocone_objects(diagram::BipartiteFreeDiagram) = ob₂(diagram)

BipartiteFreeDiagram{Ob,Hom}() where {Ob,Hom} =
  BasicBipartiteFreeDiagram{Ob,Hom}()

BipartiteFreeDiagram(obs₁::AbstractVector{Ob₁}, obs₂::AbstractVector{Ob₂},
                     homs::AbstractVector{Tuple{Hom,Int,Int}}) where {Ob₁,Ob₂,Hom} =
  BipartiteFreeDiagram{Union{Ob₁,Ob₂},Hom}(obs₁, obs₂, homs)

function BipartiteFreeDiagram{Ob,Hom}(obs₁::AbstractVector, obs₂::AbstractVector,
                                      homs::AbstractVector) where {Ob,Hom}
  @assert all(obs₁[s] == dom(f) && obs₂[t] == codom(f) for (f,s,t) in homs)
  d = BipartiteFreeDiagram{Ob,Hom}()
  add_vertices₁!(d, length(obs₁), ob₁=obs₁)
  add_vertices₂!(d, length(obs₂), ob₂=obs₂)
  add_edges!(d, getindex.(homs,2), getindex.(homs,3), hom=first.(homs))
  return d
end

BipartiteFreeDiagram(d::FixedShapeFreeDiagram{Ob,Hom}; kw...) where {Ob,Hom} =
  BipartiteFreeDiagram{Ob,Hom}(d; kw...)

function BipartiteFreeDiagram{Ob,Hom}(discrete::DiscreteDiagram;
                                      colimit::Bool=false) where {Ob,Hom}
  d = BipartiteFreeDiagram{Ob,Hom}()
  if colimit
    add_vertices₂!(d, length(discrete), ob₂=collect(discrete))
  else
    add_vertices₁!(d, length(discrete), ob₁=collect(discrete))
  end
  return d
end

function BipartiteFreeDiagram{Ob,Hom}(span::Multispan) where {Ob,Hom}
  d = BipartiteFreeDiagram{Ob,Hom}()
  v₀ = add_vertex₁!(d, ob₁=apex(span))
  vs = add_vertices₂!(d, length(span), ob₂=feet(span))
  add_edges!(d, fill(v₀, length(span)), vs, hom=legs(span))
  return d
end

function BipartiteFreeDiagram{Ob,Hom}(cospan::Multicospan) where {Ob,Hom}
  d = BipartiteFreeDiagram{Ob,Hom}()
  v₀ = add_vertex₂!(d, ob₂=apex(cospan))
  vs = add_vertices₁!(d, length(cospan), ob₁=feet(cospan))
  add_edges!(d, vs, fill(v₀, length(cospan)), hom=legs(cospan))
  return d
end

function BipartiteFreeDiagram{Ob,Hom}(para::ParallelMorphisms) where {Ob,Hom}
  d = BipartiteFreeDiagram{Ob,Hom}()
  v₁ = add_vertex₁!(d, ob₁=dom(para))
  v₂ = add_vertex₂!(d, ob₂=codom(para))
  add_edges!(d, fill(v₁,length(para)), fill(v₂,length(para)), hom=hom(para))
  return d
end

""" Convert a free diagram to a bipartite free diagram.

Reduce a free diagram to a free bipartite diagram with the same limit (the
default, `colimit=false`) or the same colimit (`colimit=true`). The reduction is
essentially the same in both cases, except for the choice of where to put
isolated vertices, where we follow the conventions described at
[`cone_objects`](@ref) and [`cocone_objects`](@ref). The resulting object is a
bipartite free diagram equipped with maps from the vertices of the bipartite
diagram to the vertices of the original diagram.
"""
function BipartiteFreeDiagram{Ob,Hom}(F::Functor{<:FinCat{Int}};
                                      colimit::Bool=false) where {Ob,Hom}
  d = FreeDiagramAsBipartite{Ob,Hom}()
  g = graph(dom(F))
  add_parts!(d, :V, nv(g))
  for v in vertices(g)
    x = ob_map(F, v)
    out_edges, in_edges = incident(g, v, :src), incident(g, v, :tgt)
    v₁ = if !isempty(out_edges) || (!colimit && isempty(in_edges))
      add_vertex₁!(d, ob₁=x, orig_vert₁=v)
    end
    v₂ = if !isempty(in_edges) || (colimit && isempty(out_edges))
      add_vertex₂!(d, ob₂=x, orig_vert₂=v)
    end
    if !(isnothing(v₁) || isnothing(v₂))
      add_edge!(d, v₁, v₂, hom=id(x))
    end
  end
  for e in edges(g)
    v₁ = incident(d, src(g, e), :orig_vert₁)
    v₂ = incident(d, tgt(g, e), :orig_vert₂)
    add_edge!(d, v₁, v₂, hom=hom_map(F, e))
  end
  return d
end
BipartiteFreeDiagram(F::Functor{<:FinCat{Int},<:Cat{Ob,Hom}}; kw...) where {Ob,Hom} =
  BipartiteFreeDiagram{Ob,Hom}(F; kw...)

# Free diagrams
#--------------

@present TheoryFreeDiagram <: TheoryGraph begin
  Ob::AttrType
  Hom::AttrType
  ob::Attr(V,Ob)
  hom::Attr(E,Hom)
end

@abstract_acset_type AbstractFreeDiagram <: AbstractGraph

@acset_type FreeDiagram(TheoryFreeDiagram, index=[:src,:tgt]) <:
  AbstractFreeDiagram

ob(d::AbstractFreeDiagram, args...) = subpart(d, args..., :ob)
hom(d::AbstractFreeDiagram, args...) = subpart(d, args..., :hom)

diagram_type(::AbstractFreeDiagram{S,T}) where {S,T<:Tuple} = T

FreeDiagram(obs::AbstractVector{Ob},
            homs::AbstractVector{Tuple{Hom,Int,Int}}) where {Ob,Hom} =
  FreeDiagram{Ob,Hom}(obs, homs)

function FreeDiagram{Ob,Hom}(obs::AbstractVector,
                             homs::AbstractVector) where {Ob,Hom}
  @assert all(obs[s] == dom(f) && obs[t] == codom(f) for (f,s,t) in homs)
  d = FreeDiagram{Ob,Hom}()
  add_vertices!(d, length(obs), ob=obs)
  length(homs) > 0 && add_edges!(d, getindex.(homs,2), getindex.(homs,3), hom=first.(homs))
  return d
end

FreeDiagram(d::FixedShapeFreeDiagram{Ob,Hom}; kw...) where {Ob,Hom} =
  FreeDiagram{Ob,Hom}(d; kw...)

function FreeDiagram{Ob,Hom}(discrete::DiscreteDiagram) where {Ob,Hom}
  d = FreeDiagram{Ob,Hom}()
  add_vertices!(d, length(discrete), ob=collect(discrete))
  return d
end

function FreeDiagram{Ob,Hom}(span::Multispan) where {Ob,Hom}
  d = FreeDiagram{Ob,Hom}()
  v₀ = add_vertex!(d, ob=apex(span))
  vs = add_vertices!(d, length(span), ob=feet(span))
  add_edges!(d, fill(v₀, length(span)), vs, hom=legs(span))
  return d
end

function FreeDiagram{Ob,Hom}(cospan::Multicospan) where {Ob,Hom}
  d = FreeDiagram{Ob,Hom}()
  vs = add_vertices!(d, length(cospan), ob=feet(cospan))
  v₀ = add_vertex!(d, ob=apex(cospan))
  add_edges!(d, vs, fill(v₀, length(cospan)), hom=legs(cospan))
  return d
end

function FreeDiagram{Ob,Hom}(para::ParallelMorphisms) where {Ob,Hom}
  d = FreeDiagram{Ob,Hom}()
  add_vertices!(d, 2, ob=[dom(para), codom(para)])
  add_edges!(d, fill(1,length(para)), fill(2,length(para)), hom=hom(para))
  return d
end

function FreeDiagram{Ob,Hom}(comp::ComposableMorphisms) where {Ob,Hom}
  d = FreeDiagram{Ob,Hom}()
  n = length(comp)
  add_vertices!(d, n+1, ob=[dom.(comp); codom(comp)])
  add_edges!(d, 1:n, 2:(n+1), hom=hom(comp))
  return d
end

function FreeDiagram{Ob,Hom}(diagram::BipartiteFreeDiagram) where {Ob,Hom}
  d = FreeDiagram{Ob,Hom}()
  vs₁ = add_vertices!(d, nv₁(diagram), ob=ob₁(diagram))
  vs₂ = add_vertices!(d, nv₂(diagram), ob=ob₂(diagram))
  add_edges!(d, vs₁[src(diagram)], vs₂[tgt(diagram)], hom=hom(diagram))
  return d
end
FreeDiagram(diagram::BipartiteFreeDiagram{Ob,Hom}) where {Ob,Hom} =
  FreeDiagram{Ob,Hom}(diagram)

function FreeDiagram{Ob,Hom}(F::Functor{<:FinCat{Int}}) where {Ob,Hom}
  diagram = FreeDiagram{Ob,Hom}()
  copy_parts!(diagram, graph(dom(F)))
  diagram[:ob] = collect_ob(F)
  diagram[:hom] = collect_hom(F)
  diagram
end
FreeDiagram(F::Functor{<:FinCat{Int},<:Cat{Ob,Hom}}) where {Ob,Hom} =
  FreeDiagram{Ob,Hom}(F)

(::Type{BFD})(diagram::FreeDiagram; kw...) where BFD <: BipartiteFreeDiagram =
  BFD(FinDomFunctor(diagram); kw...)

# FinDomFunctors as diagrams
#---------------------------

diagram_type(F::FinDomFunctor{Dom,Codom}) where {Ob,Hom,Dom,Codom<:Cat{Ob,Hom}} =
  Tuple{Ob,Hom}
cone_objects(F::FinDomFunctor) = collect_ob(F)
cocone_objects(F::FinDomFunctor) = collect_ob(F)

""" Wrapper type to interpret `FreeDiagram` as a `FinDomFunctor`.
"""
@auto_hash_equals struct FreeDiagramFunctor{Ob,Hom,Codom} <:
    FinDomFunctor{FreeCatGraph{FreeDiagram{Ob,Hom}},Codom}
  diagram::FreeDiagram{Ob,Hom}
  codom::Codom
end
FinDomFunctor(diagram::FreeDiagram, codom::Cat) =
  FreeDiagramFunctor(diagram, codom)
FinDomFunctor(diagram::FreeDiagram{Ob,Hom}) where {Ob,Hom} =
  FreeDiagramFunctor(diagram, TypeCat(Ob, Hom))

dom(F::FreeDiagramFunctor) = FreeCatGraph(F.diagram)

Categories.do_ob_map(F::FreeDiagramFunctor, x) = ob(F.diagram, x)
Categories.do_hom_map(F::FreeDiagramFunctor, f) = hom(F.diagram, f)

collect_ob(F::FreeDiagramFunctor) = ob(F.diagram)
collect_hom(F::FreeDiagramFunctor) = hom(F.diagram)

end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

const __increment_identifier = 0
const __decrement_identifier = 1
const __goto_identifier = 2
const __ifzero_goto_identifier = 3
const __halt_identifier = (4, 0)

@doc raw"""
```julia
struct Instruction <: ProgrammeCompoment
```

An instruction is a command spanning a single line in a goto programme.

This struct has the fields:
  - `I::Integer`
  - `first::Integer`
  - `second::Integer`
  - `third::Union{Integer, Nothing}`
  - `instruction::Tuple`

---

```julia
Instruction(I::Integer)
```

Given an integer, this constructor function decodes the instruction.  An instruction can be one of five things, defined as follows:
  - The code for "``Rn := Rn + 1``" is ``\left\langle 0, n\right\rangle``;
  - The code for "``Rn := Rn - 1``" is ``\left\langle 1, n\right\rangle``;
  - The code for "``\texttt{goto } k``" is ``\left\langle 2, k\right\rangle``;
  - The code for "``\texttt{if } Rn = 0 \texttt{ goto } k``" is ``\left\langle 3, \left\langle n, k\right\rangle\right\rangle``; and
  - The code for "``\texttt{halt}``" is ``\left\langle 4, 0\right\rangle``.

```julia
Instruction(instruction::Tuple)
Instruction(i::Integer, j::Integer...) = Instruction((i, j...))
```

Given a tuple or a list of values, the value of the instruction is the pair of all inputs.

---

### Examples

```julia
julia> Instruction((3, (1, 2))).I
58
```
"""
struct Instruction <: ProgrammeCompoment
	I::Integer
	first::Integer
	second::Integer
	third::Union{Integer, Nothing}
	instruction::Tuple

	function Instruction(I::Integer)
		first, second = π(I)
		instruction = nothing
		third = nothing

		if isequal(first, 3)
			second, third = π(second)
			instruction = (first, (second, third))
		else
			instruction = (first, second)
		end

		new(I, first, second, third, instruction)
	end # end constructor function

	function Instruction(instruction::Tuple)
		instruction_error = "Instructions whose coding tuple is more than three is not a valid instruction."
		length(instruction) > 3 && throw(error("$instruction_error"))
		first, second = Base.first(instruction), instruction[2]
		third = I = nothing

		if isequal(length(instruction), 3)
			third = instruction[3]
			I = pair_tuple(first, pair_tuple(second, third))
		else
			if instruction[2] isa Tuple
				second_tuple = second
				second = second_tuple[1]
				third = second_tuple[2]
				I = pair_tuple(first, pair_tuple(second, third))
			else
				I = pair_tuple(first, second)
			end
		end

		new(I, first, second, third, instruction)
	end # end constructor function

	Instruction(i::Integer, j::Integer...) = Instruction((i, j...))
end # end struct

@doc raw"""
```julia
struct Sequence <: ProgrammeCompoment
```

For any ``d \in \mathbb{N}``, the sequence code of a sequence ``(a_0,a_1,\ldots,a_{d-1}) \in \mathbb{N}^d``, denoted by ``[a_0,a_1,\ldots,a_{d-1}]``, is

```math
\left\langle d, \left\langle a a_0,a_1,\ldots,a_{d-1}\right\rangle\right\rangle
```

The struct has the following fields:
  - `q::Integer`
  - `seq_length::Integer`
  - `instructions::Tuple`

```julia
Sequence(q::Integer)
Sequence(t::Tuple)
Sequence(i::Integer, j::Tuple)
Sequence(i::Integer, j::Integer...)
```

The constructor methods for this struct decode the given value(s) into the sequence.

---

### Examples

```julia
julia> Sequence(972292871301644916468488152875266508938968846389326007980307063346008398713128885682044504108288931767348821063618087715644933567266540511345568504718733339523678538338052787779884557674350959673597803113281693069940562881722205193604550737455583875504348606989700013337656597740101535).instructions
(328, 4, 531, 4, 5, 0, 14)
```
"""
struct Sequence <: ProgrammeCompoment
	q::Integer
	seq_length::Integer
	instructions::Tuple

	function Sequence(q::Integer)
		seq_length, instructions = π(q)
		instructions = isone(seq_length) ? (instructions,) : π(instructions, seq_length)

		new(q, seq_length, instructions)
	end

	function Sequence(t::Tuple)
		q = nothing
		seq_length, instructions = Base.first(t), t[2]

		if eltype(instructions) <: Integer
			q = pair_tuple(t[1], t[2]...)
		end

		if eltype(instructions) <: Tuple
			# q = pair_tuple(Base.first(t), pair_tuple(pair_tuple.(t[2])))
			q = pair_tuple(Base.first(t), pair_tuple(t[2]))
		end

		sequence_length_error = "The first number in your sequence should match the length of your sequence."
		seq_length ≠ length(instructions) && throw(error("$(sequence_length_error)"))

		new(q, seq_length, instructions)
	end

	Sequence(i::Integer, j::Tuple) = Sequence(pair_tuple(i, pair_tuple(j...)))
	Sequence(i::Integer, j::Integer...) = Sequence((i, tuple(j...)))
end # end struct

@doc raw"""
```julia
increment(n::Integer)
```

Constructs an `Instruction` object for incrementing the ``n\textsuperscript{th}`` register.

The code for "``Rn := Rn + 1``" is ``\left\langle 0, n\right\rangle``.
"""
increment(n::Integer) = Instruction(__increment_identifier, n)

@doc raw"""
```julia
decrement(n::Integer)
```

Constructs an `Instruction` object for decrementing the ``n\textsuperscript{th}`` register.

The code for "``Rn := Rn - 1``" is ``\left\langle 1, n\right\rangle``.
"""
decrement(n::Integer) = Instruction(__decrement_identifier, n)

@doc raw"""
```julia
goto(n::Integer)
```

Constructs an `Instruction` object for going to line ``k``.

The code for "``\texttt{goto } k``" is ``\left\langle 2, k\right\rangle``.
"""
goto(k::Integer) = Instruction(__goto_identifier, k)

@doc raw"""
```julia
ifzero_goto(t::NTuple{2, Integer})
ifzero_goto(n::Integer, k::Integer)
```

Constructs an `Instruction` object for going to line ``k`` if and only  if the ``n\textsuperscript{th}`` register is zero.

The code for "``\texttt{if } Rn = 0 \texttt{ goto } k``" is ``\left\langle 3, \left\langle n, k\right\rangle\right\rangle``.
"""
ifzero_goto(t::NTuple{2, Integer}) = Instruction(__ifzero_goto_identifier, t...)
ifzero_goto(n::Integer, k::Integer) = Instruction(__ifzero_goto_identifier, (n, k)...)

@doc raw"""
```julia
halt()
```

Constructs an `Instruction` object for the final line of all goto programmes: `halt`.

The code for "``\texttt{halt}``" is ``\left\langle 4, 0\right\rangle``.
"""
halt() = Instruction(__halt_identifier...)

@doc raw"""
```julia
struct GoToProgramme <: Programme
```

This struct has fields:
  - `P::Integer`
  - `programme_length::Integer`
  - `instructions::Vector{<:Tuple}`
  - `max_line::Integer`

---

```julia
GoToProgramme(P::Integer)
GoToProgramme(S::Sequence)
```

A goto programme can be coded by a single number!  This number in the struct is `P`.

---

For any ``d \in \mathbb{N}``, the sequence code of a sequence ``(a_0,a_1,\ldots,a_{d-1}) \in \mathbb{N}^d``, denoted by ``[a_0,a_1,\ldots,a_{d-1}]``, is

```math
\left\langle d, \left\langle a a_0,a_1,\ldots,a_{d-1}\right\rangle\right\rangle
```

Then, P is ``[a_0,a_1,\ldots,a_{d-1}]``.  So, given ``P``, `π(P, 2, 0)` is the number of lines of the programme, d;
and `π(P, 2, 1)` is the programme itself.  Depairing that ``d`` times, we get a tuple of integers: the
instructions of the programme (one integer codes one line in the programme).
We define each line of the programme as being coded as follows:
  - The code for "``Rn := Rn + 1``" is ``\left\langle 0, n\right\rangle``;
  - The code for "``Rn := Rn - 1``" is ``\left\langle 1, n\right\rangle``;
  - The code for "``\texttt{goto } k``" is ``\left\langle 2, k\right\rangle``;
  - The code for "``\texttt{if } Rn = 0 \texttt{ goto } k``" is ``\left\langle 3, \left\langle n, k\right\rangle\right\rangle``; and
  - The code for "``\texttt{halt}``" is ``\left\langle 4, 0\right\rangle``.

It should also be noted that the sequence code for some base cases is defined as follows:
  - if ``d = 0``, the sequence code ``[]`` for the empty sequence is the number 0; and
  - if ``d = 1``, for ``a \in \mathbb{N}``, the sequence code ``[a]`` for the sequence ``(a)`` of length 1, is ``\left\langle 1, a\right\rangle``.

The constructor function for this struct will ensure that the given integer `P` is a valid goto programme.

---

### Examples

```julia
julia> GoToProgramme(121).length
1

julia> GoToProgramme(121).instructions
1-element Array{Tuple{BigInt,BigInt},1}:
 (4, 0)
```
"""
struct GoToProgramme <: Programme
    P::Integer
    programme_length::Integer
    instructions::Vector{<:Tuple}
    max_line::Integer

    # declare constructor function
    function GoToProgramme(P::Sequence) # P is the code for a sequence of a programme
        # Ensure the programme P is at least the nothing programme
		smallest_programme = Sequence(1, halt().I).q
        P.q < smallest_programme && throw(error("The smallest possible programme is coded by $(smallest_programme)."))

        # generate the snapshot of programme P
		snapshot = Instruction.(P.instructions)
        programme_length = P.seq_length

        # construct list of codes for each instruction
        instructions = iszero(programme_length) ? 0 : Tuple[i.instruction for i in snapshot]
        # check that the programme halts at the end
    	instructions[end] ≠ halt().instruction && throw(error("Goto programmes neccesarily have a halting instruction."))

		max_line = programme_length - 1
		lower_bound, upper_bound = extrema(Int[__increment_identifier, __decrement_identifier, __goto_identifier, __ifzero_goto_identifier, __halt_identifier...])
        row_counter = 0

        # check that all programme instruction codes are valid instructions
        for instruction in instructions
            primary_identifier, secondary_identifier = instruction

			if primary_identifier < lower_bound || primary_identifier > upper_bound || (isequal(primary_identifier, halt().first) && ! isequal(secondary_identifier, halt().second))
                throw(error("No known instruction for code ⟨$(join(instruction, ","))⟩"))
            end

            if isequal(primary_identifier, halt().first) && ! isequal(secondary_identifier, halt().second) && ! isequal(max_line, row_counter)
                throw(error("You must have exactly one halting instruction at the end of the programme."))
            end

            if isequal(primary_identifier, __goto_identifier) && secondary_identifier > max_line
                throw(error("I cannot go to line $secondary_identifier of a programme which has $max_line instructions."))
            elseif isequal(primary_identifier, __goto_identifier) && isequal(secondary_identifier, row_counter)
                throw(error("I am told to go to my own line (at line $secondary_identifier, goto line $secondary_identifier), and so I am stuck in an infinite loop.  The only way to escape is to tell you.  Please help me."))
            end

			row_counter += 1
        end

        # construct fields
        new(P.q, programme_length, instructions, max_line)
    end # end constructor (GoToProgramme) function

	# Add method for taking in an integer coding a sequence
	GoToProgramme(P::Integer) = GoToProgramme(Sequence(P))
end # end struct

"""
```julia
show_programme(io::IO, P::GoToProgramme)
show_programme(P::GoToProgramme)
show_programme(P::Integer)
```

Given a goto programme, this function will decode it into its constituent components.  It will default to `stdout`.

---

### Examples

```julia
julia> show_programme(121)
0    halt
```
"""
function show_programme(io::IO, P::GoToProgramme)
	# println("\033[1;38mThe number for $(P.P) pertains to the following programme:\033[0;38m\n")
    row_counter = 0
    message = string() # initalise message with empty string

    for instruction in P.instructions
        primary_identifier, secondary_identifier = instruction

		if isequal(primary_identifier, __increment_identifier)
            n = secondary_identifier
            message = "R$n := R$n + 1"
        elseif isequal(primary_identifier, __decrement_identifier)
            n = secondary_identifier
            message = "R$n := R$n - 1"
        elseif isequal(primary_identifier, __goto_identifier)
            k = secondary_identifier
            message = "goto $k"
        elseif isequal(primary_identifier, __ifzero_goto_identifier)
            n, k = secondary_identifier
            message = "if R$n = 0 goto $k"
        elseif isequal(primary_identifier, halt().first) && isequal(secondary_identifier, halt().second)
            message = "halt"
        else
            message = "No known instruction for code ⟨$(join(instruction, ","))⟩"
        end

        @printf(io, "%-3.3s  %-60.60s\n", "$row_counter", "$message")
		row_counter += 1
    end

    return nothing
end # end show_programme function

# Given an integer, show_programme assumes it is a programme
show_programme(io::IO, P::Integer) = show_programme(io::IO, GoToProgramme(P))

# Fall back to standard output
show_programme(P::GoToProgramme) = show_programme(stdout, P)
show_programme(P::Integer) = show_programme(stdout, GoToProgramme(P))

@doc raw"""
```julia
rand_unsafe(::Type{GoToProgramme}, d::Integer; upper_bound::Integer = 200)
```

Finds a random go-to programme.

`rand_unsafe` takes in the type, `GoToProgramme`, a number of lines of the programme, `d`, and an upper bound for the main instructions (coded; not including the `halt` instruction) which defaults to 200.  (Recall how we code instructions, using `pair_tuple`, so numbers between 1 and 200 will produce reasonably small codes)

!!! note

	*This sometimes fails, as not every random number will deconstruct into a nice go-to programme.  Programmatically, I have done as much as I can to try to avoid this failing.  It is more likely to fail as you increase the upper bound.  To avoid this, try using `rand`, though be warned it might be a little slower as it will keep trying till it finds something.

---

### Examples

```julia
julia> show_programme(ComputabilityTheory.rand_unsafe(GoToProgramme, 3, upper_bound = 200)) # a reasonably small random programme with 3 lines
0    R3 := R3 + 1
1    if R0 = 0 goto 1
2    halt
```
"""
rand_unsafe(::Type{GoToProgramme}, d::Integer; upper_bound::Integer = 200) =
	GoToProgramme(Sequence(d, (rand(1:upper_bound), halt().I)))

@doc raw"""
```julia
rand(::Type{GoToProgramme}, d::Integer; upper_bound::Integer = 200)
```

Finds a random go-to programme.

`rand` takes in the type, `GoToProgramme`, a number of lines of the programme, `d`, and an upper bound for the main instructions (coded; not including the `halt` instruction) which defaults to 200.  (Recall how we code instructions, using `pair_tuple`, so numbers between 1 and 200 will produce reasonably small codes)

!!! note

	*This may be slow*.  Depending on how large your upper bound is, it may take a while to find a valid code.

---

### Examples

```julia
julia> show_programme(rand(GoToProgramme, 3, upper_bound = 200)) # a reasonably small random programme with 3 lines
0    R3 := R3 + 1
1    if R0 = 0 goto 1
2    halt
```
"""
function Base.rand(::Type{GoToProgramme}, d::Integer; upper_bound::Integer = 200)
	try
		return rand_unsafe(GoToProgramme, d; upper_bound = upper_bound)
	catch
		return rand(GoToProgramme, d; upper_bound = upper_bound)
	end

	# should never get here
	return nothing
end


text:
using Descartes
using GeometryBasics: Mesh

p = Piping(1.0, [[0,0,0],
                 [10,0,0],
                 [10,10,0],
                 [10,10,10],
                 [5,5,5]])

m = Mesh(p)

#Descartes.visualize(m)

#save("piping.ply", m)


text:
# Advent of Code 2021, Day 5: Hydrothermal Venture
# (c) 2022 TillW

struct Point
    x::Int; y::Int
end

struct Line
    first::Point; last::Point
end

function getCoords(txt::AbstractString)::Point
    parts = split(txt, ",")
    Point(parse(Int,parts[1]), parse(Int,parts[2]))
end

function readInput(filename::String)
    lines = Line[]
    input = readlines(filename)
    for txtLine in input
        txtPoints = split(txtLine, " -> ");
        line = Line(getCoords(txtPoints[1]), getCoords(txtPoints[2]))
        lines = [lines;line]
    end
    lines
end

function addLine(area::Array{Int}, line::Line, ignoreDiagonal::Bool)
    if (line.first.x == line.last.x || line.first.y == line.last.y)
        for x = min(line.first.x, line.last.x):max(line.first.x, line.last.x)
            for y = min(line.first.y, line.last.y):max(line.first.y, line.last.y)
                area[x+1,y+1] += 1
            end
        end
    elseif !ignoreDiagonal
        left = line.first.x < line.last.x ? line.first : line.last
        right = left == line.last ? line.first : line.last
        slope = left.y < right.y ? 1 : -1
        for i = 0:(right.x - left.x)
            area[left.x+1+i, left.y+1+(i*slope)] += 1
        end
    end
end

function printArea(area::Array{Int})
    for i = 1:length(area)
        print(area[i] == 0 ? "." : area[i])
        if (i % size(area, 1) == 0) println() end
    end
end

function solvePuzzle(inputFilename::String, ignoreDiagonalLines::Bool)
    lines = readInput(inputFilename)
    maxx = maximum([map((e) -> e.first.x, lines); map((e) -> e.last.x, lines)])
    maxy = maximum([map((e) -> e.first.y, lines); map((e) -> e.last.y, lines)])
    area = zeros(Int, maxx+1, maxy+1)
    for line in lines
        addLine(area, line, ignoreDiagonalLines)
    end
    count((x) -> x > 1, area)
end

function solvePuzzle1(inputFileName::String) solvePuzzle(inputFileName, true) end
function solvePuzzle2(inputFileName::String) solvePuzzle(inputFileName, false) end

print("Number of overlaps: ")
print(solvePuzzle1("input_05.txt"))
print("/")
println(solvePuzzle2("input_05.txt"))


text:
function Molecule(;atoms::Array{Atom,1}=Array{Atom,1}())
    Molecule(atoms)
end


text:
creshape(x) = reshape(x, :, size(x)[end])
rreshape(x) = reshape(x, size(x, 1), :)

csize(a) = (ndims(a) == 1 ? size(a) : size(a)[1:end-1])
csize(a, n) = tuple(csize(a)..., n) # size if you had n columns
rsize(a) = (ndims(a) == 1 ? size(a) : size(a)[2:end])
rsize(a, n) = tuple(n, rsize(a)...) # size if you had n columns

ccount(a) = (ndims(a) == 1 ? length(a) : size(a, ndims(a)))
rcount(a) = (ndims(a) == 1 ? length(a) : size(a, 1))

@generated function subslice(x::AbstractArray{T, N}) where {T, N}
    inds = ntuple(i -> (:), N - 1)
    :($inds)
end
subslice(x) = ntuple(i -> (:), ndims(x) - 1)

cview(a, i) = view(a, subslice(a)..., i)
rview(a, i) = view(a, i, subslice(a)...)
cget(a, i) = getindex(a, subslice(a)..., i)
rget(a, i) = getindex(a, i, subslice(a)...)
cset!(a, x, i) = setindex!(a, x, subslice(a)..., i)
rset!(a, x, i) = setindex!(a, x, i, subslice(a)...)

unsqueeze(x; dims) = reshape(x, ntuple(d -> d == dims ? 1 : d < dims ? size(x, d) : size(x, d - 1), Val(ndims(x) + 1)))

function concatenate(xs::Vector{<:AbstractArray{T, N}}; dims = -1) where {T, N}
    shape = ntuple(Val(N)) do d
        d != dims ? size(xs[1], d) :
        sum(size.(xs, dims))
    end
    xc, i = zeros(T, shape), 0
    for x in xs
        is = ntuple(Val(N)) do d
            d != dims ? (1:shape[d]) :
            (i + 1):(i + size(x, dims))
        end
        copyto!(view(xc, is...), x)
        i += size(x, dims)
    end
    return xc
end

cconcatenate(xs) = concatenate(xs, dims = ndims(xs[1]))
rconcatenate(xs) = concatenate(xs, dims = 1)

stack(xs; dims) = concatenate(unsqueeze.(xs, dims = dims), dims = dims)
unstack(xs; dims) = [selectdim(xs, dims, i) for i in 1:size(xs, dims)]

cstack(xs) = stack(xs, ndims(first(xs)) + 1)
rstack(xs) = stack(xs, 1)

indbatch(x, b, offset = 0) = (C = ccount(x); min(i + offset, C):min(i + offset + b -1, C) for i in 1:b:C)
minibatch(x, batchsize) = [cview(x, ind) for ind in indbatch(x, batchsize)]

function Base.split(x::AbstractArray, n::Integer)
    cview(x, 1:n), cview(x, (n + 1):ccount(x))
end

function Base.split(x::AbstractArray, is::AbstractArray{<:Integer})
    [cview(x, i:i′) for (i, i′) in zip([0; is[1:end - 1]] .+ 1, is)]
end

colvec(x) = reshape(x, length(x), 1)
rowvec(x) = reshape(x, 1, length(x))

splat(list) = [item for sublist in list for item in (isa(sublist, AbstractArray) ? sublist : [sublist])]


text:
using Distributions
using ForwardDiff
using LinearAlgebra
using PERK
using Random
using StableRNGs
using Statistics
using Test

@testset "PERK.jl" begin
    include("kernels.jl")
    include("krr.jl")
    include("estimation.jl")
    include("complex.jl")
    include("holdout.jl")
    include("forwarddiff.jl")
end


text:
# traceradvdiff.jl's test function

import FourierFlows.TracerAdvDiff

"""
    test_constvel(; kwargs...)

Advects a gaussian concentration c0(x, y, t) with a constant velocity flow
u(x, y, t) = uvel and v(x, y, t) = vvel and compares the final state with
cfinal = c0(x-uvel*tfinal, y-vvel*tfinal)
"""
function test_constvel(stepper, dt, nsteps)

  nx, Lx = 128, 2π
  uvel, vvel = 0.2, 0.1
  u(x, y) = uvel
  v(x, y) = vvel

  prob = TracerAdvDiff.Problem(; nx=nx, Lx=Lx, kap=0.0, u=u, v=v, dt=dt, stepper=stepper)
  s, v, p, g = prob.state, prob.vars, prob.params, prob.grid

  σ = 0.1
  c0func(x, y) = 0.1*exp.(-(x.^2+y.^2)/(2σ^2))

  c0 = c0func.(g.X, g.Y)
  tfinal = nsteps*dt
  cfinal = c0func.(g.X .- uvel*tfinal, g.Y .- vvel*tfinal)

  TracerAdvDiff.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, v.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_timedependenttvel(; kwargs...)

Advects a gaussian concentration c0(x, y, t) with a time-varying velocity flow
u(x, y, t) = uvel and v(x, y, t) = vvel*sign(-t+tfinal/2) and compares the final
state with cfinal = c0(x-uvel*tfinal, y)
"""
function test_timedependentvel(stepper, dt, tfinal)

  nx, Lx = 128, 2π
  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=1e-12)
      error("tfinal is not multiple of dt")
  end

  uvel, vvel = 0.2, 0.1
  u(x, y, t) = uvel
  v(x, y, t) = t <= tfinal/2 ? vvel : -vvel

  prob = TracerAdvDiff.ConstDiffProblem(; nx=nx, Lx=Lx, kap=0.0, u=u, v=v, dt=dt, stepper=stepper)
  s, v, p, g = prob.state, prob.vars, prob.params, prob.grid

  σ = 0.1
  c0func(x, y) = 0.1*exp.(-(x.^2+y.^2)/(2σ^2))

  c0 = c0func.(g.X, g.Y)
  tfinal = nsteps*dt
  cfinal = c0func.(g.X .- uvel*tfinal, g.Y)

  TracerAdvDiff.set_c!(prob, c0func)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, v.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_diffusion(; kwargs...)

Diffuses a gaussian concentration c0(x, y, t) and compares the final state with
the analytic solution of the heat equation, cfinal
"""
function test_diffusion(stepper, dt, tfinal; steadyflow = true)

  nx  = 128
  Lx  = 2π
  kap = 0.01
  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=1e-12)
      error("tfinal is not multiple of dt")
  end

  grid = TwoDGrid(nx, Lx)
  prob = TracerAdvDiff.ConstDiffProblem(; steadyflow=steadyflow, grid=grid, nx=nx,
    Lx=Lx, kap=kap, dt=dt, stepper=stepper)
  s, v, p, g = prob.state, prob.vars, prob.params, prob.grid

  c0ampl, σ = 0.1, 0.1
  c0func(x, y) = c0ampl*exp.(-(x.^2+y.^2)/(2σ^2))

  c0 = c0func.(g.X, g.Y)
  tfinal = nsteps*dt
  σt = sqrt(2*kap*tfinal + σ^2)
  cfinal = c0ampl*σ^2/σt^2 * exp.(-(g.X.^2 + g.Y.^2)/(2*σt^2))

  TracerAdvDiff.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, v.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_hyperdiffusion(; kwargs...)

Diffuses a gaussian concentration c0(x, y, t) using hyperdiffusivity and
compares the final state with the analytic solution of the heat equation, cfinal
"""
function test_hyperdiffusion(stepper, dt, tfinal; steadyflow = true)

    nx  = 128
    Lx  = 2π
    kap = 0.0   # no diffusivity
    eta = kap   # no diffusivity
   kaph = 0.01  # hyperdiffusivity coeff
  nkaph = 1     # nkaph=1 converts hyperdiffusivity to plain diffusivity
                # so we can compare with the analytic solution of heat equation

  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=1e-12)
      error("tfinal is not multiple of dt")
  end

   g = TwoDGrid(nx, Lx)

  u, v = 0*g.X, 0*g.X

  vs = TracerAdvDiff.Vars(g)
  pr = TracerAdvDiff.ConstDiffSteadyFlowParams(eta, kap, kaph, nkaph, u, v, g)
  eq = TracerAdvDiff.Equation(pr, g)
  ts = FourierFlows.autoconstructtimestepper(stepper, dt, eq.LC, g)
  prob = FourierFlows.Problem(g, vs, pr, eq, ts)

  c0ampl, σ = 0.1, 0.1
  c0func(x, y) = c0ampl*exp.(-(x.^2+y.^2)/(2σ^2))

  c0 = c0func.(g.X, g.Y)
  tfinal = nsteps*dt
  σt = sqrt(2*kaph*tfinal + σ^2)
  cfinal = c0ampl*σ^2/σt^2 * exp.(-(g.X.^2 + g.Y.^2)/(2*σt^2))

  TracerAdvDiff.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, vs.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

# --
# Run tests
# --

stepper = "RK4"

dt, nsteps  = 1e-2, 40
@test test_constvel(stepper, dt, nsteps)

dt, tfinal  = 0.002, 0.1
@test test_timedependentvel(stepper, dt, tfinal)

dt, tfinal  = 0.005, 0.1
@test test_diffusion(stepper, dt, tfinal; steadyflow=true)

dt, tfinal  = 0.005, 0.1
@test test_diffusion(stepper, dt, tfinal; steadyflow=false)

dt, tfinal  = 0.005, 0.1
@test test_hyperdiffusion(stepper, dt, tfinal)


text:
export one_step_em,
    component_weights_per_example,
    initial_weights,
    clustering,
    log_likelihood_per_instance_per_component,
    estimate_parameters_cached,
    learn_circuit_mixture,
    learn_strudel

using Statistics: mean
using LinearAlgebra: normalize!
using Clustering: kmeans, nclusters, assignments
using DataFrames

function one_step_em(spc, data, values, flows, node2id, component_weights; pseudocount)
    # E step
    lls = log_likelihood_per_instance_per_component(spc, data, values, flows, node2id)
    lls .+= log.(component_weights)

    example_weights = component_weights_per_example(lls)
    component_weights .= sum(example_weights, dims=1)
    normalize!(component_weights, 1.0)

    # M step
    estimate_parameters_cached(spc, example_weights, values, flows, node2id; pseudocount=pseudocount)
    logsumexp(lls, 2), component_weights
end

function component_weights_per_example(log_p_of_x_and_c)
    log_p_of_x = logsumexp(log_p_of_x_and_c, 2) # marginalize out components
    log_p_of_given_x_query_c = mapslices(col -> col .- log_p_of_x, log_p_of_x_and_c, dims=[1])
    p_of_given_x_query_c = exp.(log_p_of_given_x_query_c) # no more risk of underflow, so go to linear space
    @assert sum(p_of_given_x_query_c) ≈ size(log_p_of_x_and_c, 1) "$(sum(p_of_given_x_query_c)) != $(size(log_p_of_x_and_c))"# each row has proability 1
    Matrix(p_of_given_x_query_c)
end

function initial_weights(data, mix_num::Int64; alg="cluster")::Vector{Float64}
    if alg == "cluster"
        clustered = clustering(data, mix_num)
        counting = Float64.(num_examples.(clustered))
        return normalize!(counting, 1)
    elseif alg == "random"
        return normalize!(rand(Float64, mix_num), 1)
    else
        error("Initialize weights algorithm is $undefined")
    end
end

function clustering(data, mix_num::Int64; maxiter=200)::Vector
    n = num_examples(data)
    if mix_num == 1
        return [data]
    end
    data = Matrix(data)
    R = kmeans(data', mix_num; maxiter=maxiter)
    @assert nclusters(R) == mix_num
    a = assignments(R)

    clustered_data = Vector()
    for k in 1 : mix_num
        push!(clustered_data, DataFrame(data[findall(x -> x == k, a), :]))
    end

    return clustered_data
end

function log_likelihood_per_instance_per_component(pc::SharedProbCircuit, data::DataFrame, values::Matrix{UInt64}, flows::Matrix{UInt64}, node2id)
    @assert isbinarydata(data) "Can only calculate EVI on Bool data"

    N = num_examples(data)
    num_mix = num_components(pc)
    log_likelihoods = zeros(Float64, N, num_mix)
    indices = init_array(Bool, N)::BitVector

    ll(n::SharedProbCircuit) = ()
    ll(n::SharedSumNode) = begin
        if num_children(n) != 1 # other nodes have no effect on likelihood
            for i in 1 : num_children(n)
                c = children(n)[i]
                log_theta = reshape(n.log_probs[i, :], 1, num_mix)
                indices = downflow_all(values, flows, N, n, c, node2id)
                view(log_likelihoods, indices::BitVector, :) .+=  log_theta
            end
         end
    end

    foreach(ll, pc)
    log_likelihoods
end

function estimate_parameters_cached(pc::SharedProbCircuit, example_weights::Matrix{Float64},
        values::Matrix{UInt64}, flows::Matrix{UInt64}, node2id; pseudocount::Float64)
    N = size(example_weights, 1)
    foreach(pc) do pn
        if is⋁gate(pn)
            if num_children(pn) == 1
                pn.log_probs .= 0.0
            else
                smoothed_flow = Float64.(sum(example_weights[downflow_all(values, flows, N, pn, node2id), :], dims=1)) .+ pseudocount
                uniform_pseudocount = pseudocount / num_children(pn)
                children_flows = vcat(map(c -> sum(example_weights[downflow_all(values, flows, N, pn, c, node2id), :], dims=1), children(pn))...)
                @. pn.log_probs = log((children_flows + uniform_pseudocount) / smoothed_flow)
                @assert all(sum(exp.(pn.log_probs), dims=1) .≈ 1.0) "Parameters do not sum to one locally"
                # normalize away any leftover error
                pn.log_probs .-= logsumexp(pn.log_probs, dims=1)
            end
        end
    end
end

"""
    learn_strudel(train_x; num_mix = 5, init_maxiter = 10, em_maxiter=20)

Learn a mixture of circuits
See "Strudel: Learning Structured-Decomposable Probabilistic Circuits. [arxiv.org/abs/2007.09331](https://arxiv.org/abs/2007.09331)
"""
function learn_strudel(train_x; num_mix = 5,
    pseudocount=1.0,
    init_maxiter = 10,
    em_maxiter = 20,
    verbose = true)

    pc = learn_circuit(train_x; maxiter = init_maxiter, verbose = verbose)
    learn_circuit_mixture(pc, train_x; num_mix = num_mix, pseudocount = pseudocount,
                          em_maxiter = em_maxiter, verbose = verbose)
end

"""
Given a circuit, learns a mixture of structure decomposable circuits based on that circuit

    learn_circuit_mixture(pc, data; num_mix=5, pseudocount=1.0, maxiter=20)
"""
function learn_circuit_mixture(pc, data;
        num_mix=5,
        pseudocount=1.0,
        em_maxiter=20,
        verbose = true)

    spc = compile(SharedProbCircuit, pc, num_mix)
    values, flows, node2id = satisfies_flows(spc, data)
    component_weights = reshape(initial_weights(data, num_mix), 1, num_mix)
    estimate_parameters_cached(spc, ones(Float64, num_examples(data), num_mix) ./ num_mix, values, flows, node2id; pseudocount=pseudocount)

    lls = nothing
    for iter in 1 : em_maxiter
        @assert isapprox(sum(component_weights), 1.0; atol=1e-10)
        lls, component_weights = one_step_em(spc, data, values, flows, node2id, component_weights; pseudocount=pseudocount)
        verbose && println("EM Iteration $iter/$em_maxiter. Log likelihood $(mean(lls))")
    end
    spc, component_weights, lls
end


text:
"""
    LaneFollowingAccel
Longitudinal acceleration

# Fields
- `a::Float64` longitudinal acceleration [m/s^2]
"""
struct LaneFollowingAccel
    a::Float64
end

function propagate(veh::Vehicle1D, action::LaneFollowingAccel, roadway::StraightRoadway, Δt::Float64)

    a = action.a
    s, v = veh.state.s, veh.state.v

    s′ = s + v*Δt + a*Δt*Δt/2
    v′ = v + a*Δt

    s′ = mod_position_to_roadway(s′, roadway)

    return State1D(s′, v′)
end

function propagate(veh::Entity{VehicleState,D,I}, action::LaneFollowingAccel, roadway::Roadway, ΔT::Float64) where {D,I}

    a_lon = action.a

    ds = veh.state.v

    ΔT² = ΔT*ΔT
    Δs = ds*ΔT + 0.5*a_lon*ΔT²

    v₂ = ds + a_lon*ΔT

    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    posG = roadway[roadind].pos
    posF = Frenet(roadind, roadway, t=veh.state.posF.t, ϕ=veh.state.posF.ϕ)
    VehicleState(posG, posF, v₂)
end

#XXX these should probably be removed
# Base.show(io::IO, a::LaneFollowingAccel) = @printf(io, "LaneFollowingAccel(%6.3f)", a.a)
# Base.length(::Type{LaneFollowingAccel}) = 1
# Base.convert(::Type{LaneFollowingAccel}, v::Vector{Float64}) = LaneFollowingAccel(v[1])

# function Base.copyto!(v::Vector{Float64}, a::LaneFollowingAccel)
#     v[1] = a.a
#     v
# end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo

timestep=1.0e-2
gravity=-9.81

T=Float64

# Parameters
height = 1.0

# Bodies
origin = Origin{T}()

bodies = [Sphere(0.1, 1.0, name=Symbol("ball"))]

# Joints
joints = [JointConstraint(Floating(origin, bodies[1]))]

# Contacts
contacts = [contact_constraint(bodies[1], SA[0.0; 0.0; 1.0],
            friction_coefficient=0.5,
            contact_origin=SA[0.0; 0.0; 0.0],
            contact_radius=0.1,
            contact_type=:nonlinear)]

mech = Mechanism(origin, bodies, joints, contacts,
    gravity=gravity,
    timestep=timestep)

mech.bodies[1].state.x2 = [0.0; 0.0; 0.5]

storage = simulate!(mech, 1.0,
    record=true,
    opts=SolverOptions(rtol=1.0e-6, btol=1.0e-6, verbose=true))

vis = Visualizer()
render(vis)
visualize(mech, storage,
    vis=vis)


text:
using RTMidi
using Test

@testset "RTMidi.jl" begin
    mi = MidiIn("RTMidi")
    mo = MidiOut("RTMidi")
end


text:
module FloorLayout

using JuMP, MathProgBase

export Problem,
       Solution,
       yalparser,
       load_yal,
       verify_and_print

type Problem
    N::Int
    W::Real
    H::Real
    wlb
    hlb
    wub
    hub
    area
    aspect
    c
end

type Solution
    prob::Problem
    cˣ
    cʸ
    ℓˣ
    ℓʸ
    cost::Float64
    form
    cuts
    model
    time::Float64
    nvar::Int
    ncons::Int
    stat::Symbol
    nodes::Int
    gap::Float64
    lowerbound::Float64
end

type Variables{Formulation}
    cˣ
    cʸ
    ℓˣ
    ℓʸ
    dˣ
    dʸ
    v  # "raw" binary variables, sizes/interpretations depend on Formulation
    zˣ # zˢ[i,j] is 1 iff box i sits to the left of box j in direction s,
    zʸ # and 0 otherwise
end
unpack(a::Variables) = (a.cˣ,a.cʸ,a.ℓˣ,a.ℓʸ,a.dˣ,a.dʸ,a.v,a.zˣ,a.zʸ)

include("yalparser.jl")

function validate_solution(prob::Problem, sol::Solution)
    N, W, H = prob.N, prob.W, prob.H
    x, y, w, h = sol.cˣ, sol.cʸ, sol.ℓˣ, sol.ℓʸ
    tol = 1e-3*max(W,H)
    for i in 1:N
        if !(prob.wlb[i] - tol <= w[i] <= prob.wub[i] + tol) ||
           !(prob.hlb[i] - tol <= h[i] <= prob.hub[i] + tol)
            return false
        end
        if !(w[i] * h[i] >= prob.area[i] - tol)
            return false
        end
        if x[i] + 0.5w[i] >= W + tol ||
           x[i] - 0.5w[i] <= 0 - tol ||
           y[i] + 0.5h[i] >= H + tol ||
           y[i] - 0.5h[i] <= 0 - tol
            return false
        end
    end
    for i in 1:N
        for j in (i+1):N
            if abs(x[i] - x[j]) <= 0.5*(w[i] + w[j]) - tol &&
               abs(y[i] - y[j]) <= 0.5*(h[i] + h[j]) - tol
                return false
            end
        end
    end
    return true
end

function verify_and_print(fp, prob::Problem, sol::Solution)
    valid = validate_solution(prob, sol)
    println(fp, "$(sol.form): ")
    println(fp, "cuts = $(sol.cuts)")
    sol.stat == :Optimal || print(fp, "NOT OPTIMAL ")
    valid || print(fp, "NOT VALID")
    println(fp)
    println(fp, """
        $(sol.time) seconds
        $(sol.nodes) nodes
        $(sol.cost) cost
        $(sol.lowerbound) lowerbound
        $(sol.gap) gap
        $(sol.nvar) variables
        $(sol.ncons) (linear) constraints
    """)
    flush(fp)
end

export @unpack
macro unpack(prob)
    return esc(quote
        n   = $(prob).N
        wlb = $(prob).wlb
        hlb = $(prob).hlb
        wub = $(prob).wub
        hub = $(prob).hub
        c   = $(prob).c
        area = $(prob).area
        W = $(prob).W
        H = $(prob).H
    end)
end

include("get_problem_data.jl")

include("formulations/base.jl")
include("formulations/base_cuts.jl")

include("formulations/sdp.jl")

end


text:
import PWS
using PWS: SMCEstimate, DirectMCEstimate, marginal_configuration, ConditionalEnsemble, MarginalEnsemble, gene_expression_system, generate_configuration, log_marginal, logpdf, simulate

smc = SMCEstimate(16)
dmc = DirectMCEstimate(4)

# using Distributed
# addprocs(8)

# @everywhere begin
#     import Pkg
#     Pkg.activate(".")
# end
# @everywhere import PWS

system = PWS.chemotaxis_system(mean_L=50, dtimes=0:0.1:10.0)
@time conf = generate_configuration(system)

cens = ConditionalEnsemble(system)
mens = MarginalEnsemble(system)
mconf = marginal_configuration(conf)

cr = simulate(smc, conf, cens)
mr = simulate(smc, mconf, mens)
log_marginal(cr) - log_marginal(mr)

result = PWS.mutual_information(system, smc, num_responses=10)

# @time result = PWS.run_parallel(system_fn, smc, 80)

using Plots
# using Statistics
plot(system.dtimes, result.MutualInformation, color=:gray, legend=false)
# plot!(system.dtimes, dresult.MutualInformation, color=:pink, label="")

# plot(system.dtimes, mean(result.MutualInformation), ribbon=sqrt.(var(result.MutualInformation)./size(result, 1)), label="SMC", ylabel="Path mutual information", xlabel="Trajectory length")
# plot!(system.dtimes, mean(dresult.MutualInformation), ribbon=sqrt.(var(dresult.MutualInformation)./size(dresult, 1)), label="DMC")


text:
################################################################################

module MindReader

################################################################################

# dependencies
using Parameters
using DataFrames
using Dates
using EDF
using FreqTables
using DelimitedFiles
using NamedArrays
using OrderedCollections
using StatsBase
using CSV
using FFTW

# Flux
import Flux: Chain, Dense
import Flux: mse, throttle, ADAM
import Flux: throttle, params, train!
import Flux: @epochs
import Flux.Data: DataLoader

# CairoMakie
import CairoMakie: Figure, Axis, Colorbar, Relative
import CairoMakie: heatmap!, save

################################################################################

# readEDF
export getSignals

# writeCSV
export writeHMM, writePerformance

# screening
export sensitivitySpecificity, predictiveValue

# signalBin
export extractSignalBin

# fastFourierTransform
export extractFFT

# architect
export buildAutoencoder

# shapeShifter
export shifter, reshifter

# autoencoder
export modelTrain!

# graphics
export mindGraphics

################################################################################

# declare tool directories
begin
  montageDir = "Montage/"                                  # EEG montage
  utilDir    = "Utilities/"                                # IO utilities
  performDir = "Performance/"                              # performance
  signalDir  = "SignalProcessing/"                         # preprocessing
  arqDir     = "Architect/"                                # neural network
  graphDir   = "Graphics/"                                 # graphic rendering
end;

################################################################################

# load functions
begin

  # EEG montage
  include( string(montageDir, "electrodeID.jl") )

  # IO utilities
  include( string(utilDir,    "readEDF.jl") )
  include( string(utilDir,    "writeCSV.jl") )

  # performance
  include( string(performDir, "screening.jl") )

  # preprocessing
  include( string(signalDir,  "signalBin.jl") )
  include( string(signalDir,  "fastFourierTransform.jl") )

  # neural network
  include( string(arqDir,     "architect.jl") )
  include( string(arqDir,     "shapeShifter.jl") )
  include( string(arqDir,     "autoencoder.jl") )

  # graphic rendering
  include( string(graphDir,   "statesHeatMap.jl") )

end;

################################################################################

end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
println("Running tests for Jags-j1.0.0-v2.0.0:")

using Compat, Jags
using Test

code_tests = ["test_cmd.jl"]

println("Run execution_tests only if Jags is available")
execution_tests = [
  "test_line1.jl",
  "test_line2.jl",
  "test_line3.jl",
  "test_line4.jl",
  "test_rats.jl",
  "test_bones1.jl",
  "test_bones2.jl",
  "test_dyes.jl"
]

for my_test in code_tests
    println("\n  * $(my_test) *")
    include(my_test)
end

if isdefined(Main, :JAGS_HOME) && length(JAGS_HOME) > 0
  @testset "Jags.jl" begin

    for my_test in code_tests
        println("\n\n\n  * $(my_test) *")
        include(my_test)
    end

    for my_test in execution_tests
        println("\n\n\n  * $(my_test) *")
        include(my_test)
    end

    println("\n")

  end
else
  println("\n\nJAGS_HOME not found.")
  println("Skipping all tests that depend on the Jags executable!\n")
end


text:
#!/usr/bin/env julia

using Luxor, Colors

function drawbbox(apoly)
    gsave()
    setline(0.3)
    setdash("dotted")
    box(polybbox(apoly), :stroke)
    grestore()
end

function test1(p, x, y)
    prettypoly(p, :fill,
              # all these commands are executed for each vertex of the polygon
              :(
              scale(0.5, 0.5);
              randomhue();
              opacity = currentdrawing.alpha;
              setopacity(rand(5:10)/10);
              ngon(0, 0, 15, rand(3:12), 0, :fill);
              setopacity(opacity);
              sethue("black");
              # we can pass the origin for the whole polygon into the expression like this
              # but the current vertex is 0/0
              textcentred(string($(x)) * "/" *  string($(y)))))
end

function test2(p)
    prettypoly(p, :fill,
                  :(
                  ## all these commands are executed for each vertex of the polygon
                  randomhue();
                  scale(0.5, 0.5);
                  rotate(pi/2);
                  prettypoly($(p), :fill,
                    :(
                    # and all these commands are executed for each vertex of that polygon
                    randomhue();
                    scale(0.25, 0.25);
                    rotate(pi/2);
                    prettypoly($($(p)), :fill)))))
end

function get_png_files(folder)
    cd(folder)
    imagelist = filter(f -> !startswith(f, ".") && endswith(f, "png"), readdir(folder))
    imagelist = filter(f -> !startswith(f, "tiled-images"), imagelist) # don't recurse... :)
    return map(realpath, imagelist)
end

function test3(p)
    imagelist = get_png_files(dirname(@__FILE__))
    shuffle!(imagelist)
    img = readpng(imagelist[1])
    w = img.width
    h = img.height
    prettypoly(p, :fill,
                  :(
                  ## all these commands are executed for each vertex of the polygon
                  # to make it "through" the expression to the base drawing functions, these have to be made literal using $
                  placeimage($readpng($imagelist[rand(1:end)]), -$w/2, -$h/2)
                  ))
end

function draw_lots_of_polys(pagewidth, pageheight)
    setopacity(0.75)
    pagetiles = Tiler(width, height, 9, 9, margin=20)
    for (pos, n) in pagetiles
        if rand(Bool)
            p = randompointarray(rand(-140:-10), rand(-140:-10), rand(10:140), rand(10:140), rand(5:12))
        else
            p = ngon(0, 0, rand(50:100), rand(3:12), vertices=true)
        end
        gsave()
        translate(pos)
        setline(1)
        randomhue()
        @eval [test1($p, $pos.x, $pos.y), test2($p), test3($p)][rand(1:end)] # choose a test at random
        drawbbox(p)
        grestore()
    end
end

width, height = 3000, 3000
fname = "/tmp/pretty-poly-test.pdf"
Drawing(width, height, fname)
origin()
background("ivory")
draw_lots_of_polys(width, height)
finish()
println("finished test: output in $(fname)")