text:
const LessThanIndicatorSetOne{T} = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE, MOI.LessThan{T}}
const LessThanIndicatorSetZero{T} = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO, MOI.LessThan{T}}

# Needed by test spread over several files, defining it here make it easier to comment out tests
# Model supporting every MOI functions and sets

MOIU.@model(Model,
            (MOI.ZeroOne, MOI.Integer),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval,
             MOI.Semicontinuous, MOI.Semiinteger),
            (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
             MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
             MOI.GeometricMeanCone, MOI.ExponentialCone, MOI.DualExponentialCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.PositiveSemidefiniteConeSquare,
             MOI.RootDetConeTriangle, MOI.RootDetConeSquare, MOI.LogDetConeTriangle,
             MOI.LogDetConeSquare),
            (MOI.PowerCone, MOI.DualPowerCone, MOI.SOS1, MOI.SOS2, LessThanIndicatorSetOne, LessThanIndicatorSetZero),
            (), (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))


text:
using ArrayLayouts, BlockBandedMatrices, BandedMatrices, BlockArrays, LinearAlgebra, Test
import Base: getindex, size
import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedStyle, bandeddata, BandedColumns

@testset "Diagonal interface" begin
    n = 10
    D = Diagonal(randn(n^2))
    @test blockbandwidths(D) == subblockbandwidths(D) == (0,0)

    PD = PseudoBlockArray(D, Fill(n,n), Fill(n,n))
    @test blockbandwidths(PD) == bandwidths(PD) == (0,0)
    @test MemoryLayout(typeof(PD)) isa DiagonalLayout{DenseColumnMajor}
    @test bandeddata(PD) == bandeddata(D)

    V = view(PD, Block(1,1))
    @test bandwidths(V) == (0,0)
    @test Broadcast.BroadcastStyle(typeof(V)) == BandedStyle()
    @test MemoryLayout(typeof(V)) isa BandedColumns{DenseColumnMajor}
    @test bandeddata(V) == bandeddata(PD)[:,1:n]
    @test BandedMatrix(V) == V
end

@testset "Block Diagonal" begin
    A = BlockBandedMatrices.BlockDiagonal(fill([1 2],3))
    @test blockbandwidths(A) == (0,0)
    @test isblockbanded(A)
    @test A[Block(1,1)] == [1 2]
    @test @inferred(getblock(A,1,2)) == @inferred(A[Block(1,2)]) == [0 0]
    @test_throws DimensionMismatch A+I
    A = BlockBandedMatrices.BlockDiagonal(fill([1 2; 1 2],3))
    @test A+I == I+A == mortar(Diagonal(fill([2 2; 1 3],3))) == Matrix(A) + I
end

@testset "Block Bidiagonal" begin
    Bu = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :U)
    Bl = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :L)
    @test blockbandwidths(Bu) == (0,1)
    @test blockbandwidths(Bl) == (1,0)
    @test isblockbanded(Bu)
    @test isblockbanded(Bl)
    @test Bu[Block(1,1)] == Bl[Block(1,1)] == [1 2]
    @test @inferred(Bu[Block(1,2)]) == @inferred(Bl[Block(2,1)]) == [3 4]
    @test @inferred(getblock(Bu,1,3)) == @inferred(Bu[Block(1,3)]) == [0 0]
    @test_throws DimensionMismatch Bu+I
    Bu = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :U)
    Bl = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :L)
    @test Bu+I == I+Bu == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) + I
    @test Bl+I == I+Bl == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :L)) == Matrix(Bl) + I
    @test Bu-I == mortar(Bidiagonal(fill([0 2; 1 1],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) - I
    @test I-Bu == mortar(Bidiagonal(fill([0 -2; -1 -1],4), fill(-[3 4; 3 4],3), :U)) == I - Matrix(Bu)
end

@testset "Block Tridiagonal" begin
    A = BlockTridiagonal(fill([1 2],3), fill([3 4],4), fill([4 5],3))
    @test blockbandwidths(A) == (1,1)
    @test isblockbanded(A)
    @test A[Block(1,1)] == [3 4]
    @test @inferred(A[Block(1,2)]) == [4 5]
    @test @inferred(getblock(A,1,3)) == @inferred(A[Block(1,3)]) == [0 0]
    @test_throws DimensionMismatch A+I
    A = BlockTridiagonal(fill([1 2; 1 2],3), fill([3 4; 3 4],4), fill([4 5; 4 5],3))
    @test A+I == I+A == mortar(Tridiagonal(fill([1 2; 1 2],3), fill([4 4; 3 5],4), fill([4 5; 4 5],3))) == Matrix(A) + I
    @test A-I == mortar(Tridiagonal(fill([1 2; 1 2],3), fill([2 4; 3 3],4), fill([4 5; 4 5],3))) == Matrix(A) - I
    @test I-A == mortar(Tridiagonal(fill(-[1 2; 1 2],3), fill([-2 -4; -3 -3],4), fill(-[4 5; 4 5],3))) == I - Matrix(A)
end

ext: jl


text:
# Type for specifying paths through nested dictionaries.
struct KeyChain{T}
    keys::T
end
KeyChain(args...) = KeyChain(args)

function Base.getindex(d::Dict, k::KeyChain)
    # Using the "pass-by-sharing" semantics, argument `d` will not be changed.
    for i in k.keys
        d = d[i]
    end
    return d
end

function Base.setindex!(d::Dict, x, k::KeyChain)
    for i in 1:length(k.keys)-1
        d = d[i]
    end
    d[endof(k.keys)] = x
end

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

type_sanitize(::Type{T}, v::T) where T = v
function type_sanitize(::Type{T}, v::U) where {T,U}
    throw(TypeError(:type_sanitize, "Unexpected type for link definitions",T,U))
end

"""
    getkeys(d::T, keys, required = true) where T <: Dict

Return a dictionary `r` of type `T` with just the requested keys and
corresponding values from `d`. If `required = true`, throw `KeyError` if a key
`k` is not found. Otherwise, set `r[k] = missing`.
"""
function getkeys(d::T, keys) where T <: Dict
    r = T()
    for k in keys
        !haskey(d,k) && throw(KeyError(k))
        r[k] = d[k]
    end
    return r
end

################################################################################
stripped_contents(dir::String) = [first(gzsplitext(i)) for i in readdir(dir)]
number_regex(str) = Regex("(?<=$(str)_)\\d+")

function append_suffix(iter, key::String)
    rgx = number_regex(key)

    matches = Int[]
    for k in iter
        attach!(matches, match(rgx, k))
    end

    return length(matches) == 0 ? "$(key)_1" : "$(key)_$(1+maximum(matches))"
end

attach!(a, m) = nothing
function attach!(a::Vector{T}, m::RegexMatch) where T
    val = tryparse(T, m.match)
    if !isnull(val)
        push!(a, val.value)
    end
end

# Helpful for dealing with files with an additional .gz extension
function gzsplitext(s)
    y,z = splitext(s)
    if z == ".gz"
        x,y = splitext(y)
        return x, y*z
    end
    return y,z
end

"""
    augment(dir::String, new::String)

Add a numeric suffix to `new` so it does not conflict with anything in directory
`dir`. Create `dir` if it does not exist.
"""
function augment(dir::String, new::String)
    dir = isempty(dir) ? "." : dir
    ispath(dir) || mkpath(dir)

    prefix, ext = gzsplitext(new)
    newprefix = append_suffix(stripped_contents(dir), prefix)

    return joinpath(dir, newprefix*ext)
end


text:
#
#      V(k,n) - The Stiefel-Manifold
#  A point on this manifold is a n×k matrix.
#  A point in an tangential space is a n×k matrix.
#
import LinearAlgebra: norm, dot, nullspace, det, tr, qr, triu, rank, svd, diag, sylvester, lyap, Diagonal
import Base: exp, log, show, cat, rand, Matrix, real, atan
export Stiefel, StPoint, StTVector, getValue
export dot, exp, log, manifoldDimension, norm, parallelTransport, randomTVector, randomMPoint, retractionQR, retractionPolar, inverseRetractionPolar, inverseRetractionQR, projection, retraction, inverseRetraction
export zeroTVector
#
# Type definitions
#

@doc doc"""
    Stiefel{T<:Union{U, Complex{U}} <: Manifold

The manifold $\mathcal M = \mathrm{St}(k,n)$ represented by $n\times k$
orthogonal matrices, that are either real- or complex-valued.

# Abbreviation
  `St`

# Constructor
     Stiefel(k, n[, d=Float64])

generate the manifold $\mathrm{St}(k,n)$ where the integer `n` is the
number of rows and `k` is the number of columns of the matrices and the optional
parameter `d` sets the `DataType` of the matrix entries.
"""
struct Stiefel{T<:Union{U, Complex{U}} where U<:AbstractFloat} <: Manifold
  name::String
  dimensionlines::Int
  dimensioncolumns::Int
  abbreviation::String
  function Stiefel{T}(dimensioncolumns::Int, dimensionlines::Int) where T<:Union{U, Complex{U}} where U<:AbstractFloat
    if dimensioncolumns > dimensionlines
      throw(ErrorException("dimensioncolumns can't be bigger than dimensionlines: $dimensioncolumns > $dimensionlines"))
    else
     new("Stiefel-Manifold St($dimensioncolumns,$dimensionlines) in $T", dimensionlines, dimensioncolumns,"V($dimensioncolumns,$dimensionlines)")
    end
  end
end
Stiefel(dimensionlines::Int, dimensioncolumns::Int, D::DataType=Float64) = Stiefel{D}(dimensionlines::Int, dimensioncolumns::Int)

@doc doc"""
    StPoint <: MPoint

A point $x$ on the manifold $\mathcal M = \mathrm{St}(k,n)$ is represented by an
orthogonal matrix from $\mathbb{K}^{n\times k}.$

# Constructor
    StPoint(Matrix)

where Matrix is an orthogonal matrix of dimension $n×k$.
"""
struct StPoint{T<:Union{U, Complex{U}} where U<:AbstractFloat} <: MPoint
  value::Matrix{T}
  StPoint{T}(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = new(value)
end
getValue(x::StPoint) = x.value;
StPoint(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = StPoint{T}(value)

@doc doc"""
    StTVector <: TVector

A tangent vector $\xi \in T_x\mathcal M$ on the manifold
$\mathcal M = \mathrm{St}(k,n)$. The tangent space is given by as

$T_x\mathrm{St}(k,n) = \bigl\{\xi \in \mathbb{K}^{n\times k} \big| x^{\mathrm{T}}ξ+ξ^{\mathrm{T}}x=0 \bigr\}$.

# Constructor

    StTVector(ξ)

where `ξ` is an $n\times k$ `Matrix` that satisfies the above.
"""
struct StTVector{T<:Union{U, Complex{U}} where U<:AbstractFloat} <: TVector
  value::Matrix{T}
  StTVector{T}(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = new(value)
end
getValue(ξ::StTVector) = ξ.value;
StTVector(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = StTVector{T}(value)

#
# Traits
#
@traitimpl IsMatrixM{Stiefel}
@traitimpl IsMatrixP{StPoint}
@traitimpl IsMatrixTV{StTVector}

#
# Functions
#
@doc doc"""
    dot(M,x,ξ,ν)

compute the Riemannian inner product for two [`StTVector`](@ref)s `ξ` and `ν`
from $T_x\mathcal M$ of the [`Stiefel`](@ref) manifold `M` given by

$\langle \xi, \nu \rangle_x = \operatorname{trace}({\bar \xi}^{\mathrm{T}}\nu).$
"""
function dot(M::Stiefel{T}, x::StPoint{T}, ξ::StTVector{T}, ν::StTVector{T})::Float64 where T<:Union{U, Complex{U}} where U<:AbstractFloat
  return real(tr(getValue(ξ)'*getValue(ν)))
end

@doc doc"""
    exp(M,x,ξ [,t=1.0])

compute the exponential map on the [`Stiefel`](@ref) manifold `M`$=\mathrm{SO}(n)$ with
respect to the [`StPoint`](@ref) `x` and the [`StTVector`](@ref) `ξ`, which can
be shortened with `t` to `tξ`. The formula reads

$\operatorname{exp}_{x} tξ = \begin{pmatrix}
   x\\t\xi
 \end{pmatrix}
 \operatorname{Exp}
 \left(
 \begin{pmatrix} {\bar x}^{\mathrm{T}}\xi & -{\bar \xi}^{\mathrm{T}}\xi\\
 I_{k×k} & {\bar x}^{\mathrm{T}}\xi\end{pmatrix}
 \right)
\begin{pmatrix}  \operatorname{Exp}( -{\bar x}^{\mathrm{T}}\xi) \\ 0_{k×k}\end{pmatrix}$

where $\operatorname{Exp}$ denotes matrix exponential, and $I_{k×k}$ and
$0_{k×k} are the identity matrix and the zero matrix of dimension $k×k$,
respectively.
"""
function exp(M::Stiefel{T},x::StPoint{T},ξ::StTVector{T},t::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  Z = zeros(M.dimensioncolumns,M.dimensioncolumns)
  I = one(Z)
  Ξ = t*getValue(ξ)
  Y = [getValue(x) Ξ] * exp([getValue(x)'*Ξ -Ξ'*Ξ; I getValue(x)'*Ξ]) * [exp(-getValue(x)'*Ξ); Z]
  StPoint{T}(Y)
end

@doc doc"""
    inverseRetractionPolar(M,x,y)

return a [`StTVector`](@ref) `ξ` of the tagent space $T_x\mathrm{SO}(n)$
of the [`StPoint`](@ref) `x` on the [`Stiefel`](@ref) manifold `M`
with which the [`StPoint`](@ref) `y` can be reached by the
[`retractionPolar`](@ref)  after time 1.
The formula reads

$ξ = ys-x$

where $s$ is the solution to the Lyapunov equation

$x^{\mathrm{T}}ys + s(x^{\mathrm{T}}y)^{\mathrm{T}} + 2\mathrm{I}_k = 0.$

This function is implemented only for the case $\mathbb{K}=\mathbb{R}$.
"""
function inverseRetractionPolar(M::Stiefel{T},x::StPoint{T},y::StPoint{T}) where T<:AbstractFloat
  A = transpose(getValue(x)) * getValue(y)
  H = -2 * one(transpose(getValue(x))*getValue(x))
  B = lyap(A, H)
  C = getValue(y)*B - getValue(x)
  StTVector{T}( C )
end

@doc doc"""
    inverseRetractionQR(M,x,y)

return a [`StTVector`](@ref) `ξ` of the tagent space $T_x\mathrm{SO}(n)$
of the [`StPoint`](@ref) `x` on the [`Stiefel`](@ref) manifold `M`
with which the [`StPoint`](@ref) `y` can be reached by the
[`retractionQR`](@ref) after time 1.
This function is implemented only for the case $\mathbb{K}=\mathbb{R}$.
This is also the standard retraction.
"""
function inverseRetractionQR(M::Stiefel{T},x::StPoint{T},y::StPoint{T}) where T<:AbstractFloat
  A = transpose(getValue(x)) * getValue(y)
  R = zeros(M.dimensioncolumns, M.dimensioncolumns)
  for i = 1:M.dimensioncolumns
    b = zeros(i)
    b[end] = 1
    b[1:(end-1)] = - transpose(R[1:(i-1), 1:(i-1)]) * A[i, 1:(i-1)]
    R[1:i, i] = A[1:i, 1:i] \ b
  end
  C =  getValue(y)*R-getValue(x)
  StTVector{T}( C )
end

function inverseRetraction(M::Stiefel{T},x::StPoint{T},y::StPoint{T}) where T<:AbstractFloat
  return inverseRetractionQR(M,x,y)
end

@doc doc"""
    manifoldDimension(x)

return the dimension of the [`Stiefel`](@ref) manifold `M`$= \mathrm{St}(k,n)$, the
[`StPoint`](@ref) `x`, itself embedded in $\mathbb R^{n\times k}$, belongs to.
The dimension for $\mathbb{K}=\mathbb{R}$ is given by

$nk - \frac{1}{2}k(k+1)$

and for $\mathbb{K}=\mathbb{C}$

$2nk - k^2.$
"""
function manifoldDimension(x::StPoint{T}) where T<:AbstractFloat
  return Int(size(getValue(x),1) * size(getValue(x),2) - 0.5 * size(getValue(x),2) * (size(getValue(x),2) + 1))
end

function manifoldDimension(x::StPoint{T}) where T<:Complex{U} where U<:AbstractFloat
  return Int(2 * size(getValue(x),1) * size(getValue(x),2) - size(getValue(x),2)^2)
end

@doc doc"""
    manifoldDimension(M)

return the dimension of the [`Stiefel`](@ref) manifold `M`.
The dimension for $\mathbb{K}=\mathbb{R}$ is given by

$nk - \frac{1}{2}k(k+1)$

and for $\mathbb{K}=\mathbb{C}$

$2nk - k^2.$
"""
function manifoldDimension(M::Stiefel{T}) where T<:AbstractFloat
  return Int(M.dimensionlines * M.dimensioncolumns - 0.5 * M.dimensioncolumns * (M.dimensioncolumns + 1))
end
function manifoldDimension(M::Stiefel{T}) where T<:Complex{U} where U<:AbstractFloat
  return Int(2 * M.dimensionlines * M.dimensioncolumns - (M.dimensioncolumns)^2)
end

@doc doc"""
    norm(M,x,ξ)

compute the norm of the [`StTVector`](@ref) `ξ` in the tangent space
$T_x\mathcal M$ at [`StPoint`](@ref) `x` of the [`Stiefel`](@ref) manifold `M`.

$\lVert \xi \rVert_x = \sqrt{\sum_{i,j=0}^n \xi_{ij}^2}$

where $\xi_{ij}$ are the entries of the matrix `ξ`, i.e. the norm is given by
the Frobenius norm of `ξ`.
"""
function norm(M::Stiefel{T}, x::StPoint{T}, ξ::StTVector{T})::Float64 where T<:Union{U, Complex{U}} where U<:AbstractFloat
  norm(getValue(ξ))
end

@doc doc"""
    parallelTransport(M,x,y,ξ)

compute the paralllel transport of the [`StTVector`](@ref) `ξ` from
the tangent space $T_x\mathcal M$ at [`StPoint`](@ref) `x` to
$T_y\mathcal M$ at [`StPoint`](@ref) `y` on the [`Stiefel`](@ref) manifold `M`
provided that the corresponding [`geodesic`](@ref) $g(\cdot;x,y)$ is unique.
The formula reads

$P_{x\to y}(\xi) = \operatorname{proj}_{\mathcal M}(y,\xi).$

where $\operatorname{proj}_{\mathcal M}$ is the projection onto the
tangent space $T_y\mathcal M$.
"""
function parallelTransport(M::Stiefel{T}, x::StPoint{T}, y::StPoint{T}, ξ::StTVector{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
    projection(M,y,getValue(ξ))
end

@doc doc"""
    projection(M,x,q)

project a `Matrix q` orthogonally on the tangent space of the
[`StPoint`](@ref) `x` on the [`Stiefel`](@ref) manifold `M`. The formula reads

$\operatorname{proj}_{\mathcal M}(x,q) = q-xB$

where

$B=\frac{1}{2} (x^{\mathrm{T}}{\bar q})^{\mathrm{T}} {\bar x}^{\mathrm{T}}q.$

# see also
[`parallelTransport`](@ref), [`randomTVector`](@ref)
"""
function projection(M::Stiefel{T}, x::StPoint{T}, q::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
    A = getValue(x)'*q
    B = 0.5 * (A' + A)
    return StTVector{T}(q - getValue(x) * B)
end

@doc doc"""
    randomMPoint(M [,:Gaussian, σ=1.0])

return a random (Gaussian) [`StPoint`](@ref) `x` on the manifold
[`Stiefel`](@ref) manifold `M` by generating a (Gaussian) matrix with standard deviation
`σ` and return the orthogonalized version, i.e. return the Q
component of the QR decomposition of the random matrix of size $n×k$.
"""
function randomMPoint(M::Stiefel{T}, ::Val{:Gaussian}, σ::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
    A = σ*randn(T, (M.dimensionlines, M.dimensioncolumns))
    return StPoint{T}( Matrix(qr(A).Q) )
end

@doc doc"""
    randomTVector(M,x [,:Gaussian, σ=1.0])

return a random vector [`StTVector`](@ref) in the tangential space
$T_x\mathrm{St}(k,n)$ by generating a random matrix of size $n×k$ and projecting
it onto [`StPoint`](@ref) `x` with [`projection`](@ref).
"""
function randomTVector(M::Stiefel{T}, x::StPoint{T}, ::Val{:Gaussian}, σ::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  y = σ * randn(T, (M.dimensionlines,M.dimensioncolumns))
  Y = projection(M, x, y)
  return 1/(norm(getValue(Y))) * Y
end

@doc doc"""
    retractionPolar(M,x,ξ,[t=1.0])

move the [`StPoint`](@ref) `x` in the direction of the [`StTVector`](@ref) `ξ`
on the  [`Stiefel`](@ref) manifold `M`. This SVD-based retraction is a second-order
approximation of the exponential map [`exp`](@ref). Let

$USV = x + tξ$

be the singular value decomposition, then the formula reads

$\operatorname{retr}_x\xi = UV^\mathrm{T}.$

This function is implemented only for the case $\mathbb{K}=\mathbb{R}$.
"""
function retractionPolar(M::Stiefel{T},x::StPoint{T},ξ::StTVector{T},t::Float64=1.0) where T<:AbstractFloat
  y = getValue(x) + t * getValue(ξ)
  S = svd(y)
  A = S.U * transpose(S.V)
  StPoint{T}(A)
end

@doc doc"""
    retractionQR(M,x,ξ,[t=1.0])

move the [`StPoint`](@ref) `x` in the direction of the [`StTVector`](@ref) `ξ`
on the  [`Stiefel`](@ref) manifold `M`. This QR-based retraction is a
first-order approximation of the exponential map [`exp`](@ref). Let

$QR = x + tξ$

be the QR decomposition, then the formula reads

$\operatorname{retr}_x\xi = QD$

where D is a $n×k$ matrix with the signs of the diagonal entries of $R$ plus
$0.5$ on the upper diagonal, i.e.

$D_{ij}=\begin{cases} \operatorname{sgn}(R_{ij}+0,5) & \text{if} \; i=j \\ 0 & \, \text{otherwise} \end{cases}.$

This is also the standard retraction.
"""
function retractionQR(M::Stiefel{T},x::StPoint{T},ξ::StTVector{T},t::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  y = getValue(x) + t * getValue(ξ)
  QRdecomp = qr(y)
  d = diag(QRdecomp.R)
  D = Diagonal( sign.( sign.(d .+ 0.5) ) )
  B = zeros(M.dimensionlines,M.dimensioncolumns)
  B[1:M.dimensioncolumns,1:M.dimensioncolumns] = D
  A = QRdecomp.Q * B
  StPoint{T}(A)
end

function retraction(M::Stiefel{T},x::StPoint{T},ξ::StTVector{T},t::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  return retractionQR(M,x,ξ,t)
end

@doc doc"""
    ξ = zeroTVector(M,x)

return a zero vector in the tangent space $T_x\mathcal M$ of the
[`StPoint`](@ref) $x\in\mathrm{St}(k,n)$ on the [`Stiefel`](@ref) manifold `M`.
"""
function zeroTVector(M::Stiefel{T},x::StPoint{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  return StTVector{T}( zero(getValue(x)) )
end

@doc doc"""
    validateMPoint(M,x)

validate that the [`StPoint`](@ref) `x` is a valid point on the
[`Stiefel`](@ref) manifold `M`, i.e. that both dimensions and the rank are correct
as well as that all columns are orthonormal.
"""
function validateMPoint(M::Stiefel{T}, x::StPoint{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  if size(getValue(x), 1) ≠ M.dimensionlines
    throw( ErrorException("The dimension of $x must be $(M.dimensionlines) × $(M.dimensioncolumns) but it is $(size(getValue(x), 1)) × $(size(getValue(x), 2)).") )
  end
  if size(getValue(x), 2) ≠ M.dimensioncolumns
    throw( ErrorException("The dimension of $x must be $(M.dimensionlines) × $(M.dimensioncolumns) but it is $(size(getValue(x), 1)) × $(size(getValue(x), 2)).") )
  end
if rank(getValue(x)) ≠ M.dimensioncolumns #is this necessary?
    throw( ErrorException("$x must have  full rank.") )
  end

  z=getValue(x)'*getValue(x)
  if z ≉ one(z)
    throw( ErrorException("$x has to be orthonormal but it's not.") )
  end
end
@doc doc"""
    validateMPoint(M,x,ξ)

validate that the [`StTVector`](@ref) `ξ` is a valid tangent vector to
[`StPoint`](@ref) `x` on the [`Stiefel`](@ref) manifold `M`, i.e. that both
dimensions agree and ${\bar x}^{\mathrm{T}}$ is skew symmetric.
"""
function validateTVector(M::Stiefel{T}, x::StPoint{T}, ξ::StTVector{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  if size(getValue(ξ), 1) ≠ M.dimensionlines
    throw( ErrorException("The dimension of $ξ must be $(M.dimensionlines) × $(M.dimensioncolumns) but it is $(size(getValue(ξ), 1)) × $(size(getValue(ξ), 2))") )
  end
  if size(getValue(ξ), 2) ≠ M.dimensioncolumns
    throw( ErrorException("The dimension of $ξ must be $(M.dimensionlines) × $(M.dimensioncolumns) but it is $(size(getValue(ξ), 1)) × $(size(getValue(ξ), 2))") )
  end
  if norm(getValue(x)'*getValue(ξ) + getValue(ξ)'*getValue(x)) > 10^(-15)
    throw( ErrorException("The matrix $x'$ξ must be skew-symmetric!") )
  end
end


text:
module TEST_CIRULARSYS
using Samsara, Test

@testset "circularsys initiation" begin
    CircularSys = Samsara.CircularSys
    current_index = Samsara.current_index
    current_state = Samsara.current_state
    step_up! = Samsara._step_up!
    step_down! = Samsara._step_down!

    case = CircularSys()
    @test isa(case, AbstractSystem)

    # default antall noder: 3 zeroes
    @test length(case._all_nodes) == 3
    @test length(case) == 3
    @test case._all_nodes[1] == :A  # default list of nodes: [:A, :B, :C]

    # current state: The current node in a Tuple
    case = CircularSys()
    @test system_state(case) == (:A, )
    " system_state(case) returns the value of the current node "

    case._current_index = 1
    @test step_up!(case) == 2
    " step_up! steps one 'up' "

    case._current_index = 3
    @test current_index(case) == 3
    @test step_up!(case) == 1
    @test current_index(case) == 1
    " step_up! steps around to beginning when rounding the end => Circular! "

    case = CircularSys( nodes=(1,2,3,4) )
    @test current_index(case) == 1
    @test step_up!(case) == 2
    @test step_up!(case) == 3
    @test step_up!(case) == 4
    @test step_up!(case) == 1
    " integration test: Step up works for a full circle "

    case = CircularSys( nodes=(:state1, :state2) )
    case._current_index = 2
    @test step_down!(case) == 1
    step_down!(case)
    @test current_index(case) == 2
    @test system_state(case) == (:state2, )
    @test current_state(case) == :state2
    " Full round integration test for step_down! (on length(case) == 2) "
end

@testset "System testing: circular system with deterministic actions" begin
    CircularSys    = Samsara.CircularSys
    current_state  = Samsara.current_state
    current_action = Samsara.current_action
    set_action_in!  = Samsara.set_action_in!

    case_system = CircularSys( nodes=(:s1, :s2, :s3) )
    @assert current_state(case_system) == :s1
    @assert isnothing(case_system._latent_variables)
    @test isnothing(current_action(case_system))
    " Sircular system have _latent_variables of type action "

    @assert current_state(case_system) == :s1
    case_system._latent_variables = :up
    step_system_mechanics!(case_system)
    @test current_state(case_system) == :s2
    case_system._latent_variables = :down
    step_system_mechanics!(case_system)
    @test current_state(case_system) == :s1
    " System mechanics works as dictated by _latent_variables i.e. action "

    set_action_in!(case_system, :up)
    @test current_action(case_system) == :up
    set_action_in!(case_system, :down)
    @test current_action(case_system) == :down
    @test_throws ArgumentError set_action_in!(case_system, :everything_else)
    #@test current_action(case_system) == nothing
    " helper-funciton set_action_in!(System, action) "

    set_action_in!(case_system, :up)
    @assert current_state(case_system) == :s1
    @test step_system_mechanics!(case_system) == :s2
    set_action_in!(case_system, :down)
    @test step_system_mechanics!(case_system) == :s1
    set_action_in!(case_system, nothing)
    @test step_system_mechanics!(case_system) == :s1
    " System works as expected for actions :up and :down, and nothing "
end

@testset "Action space : configuable to be whatever" begin
    set_action_in! = Samsara.set_action_in!
    current_action = Samsara.current_action
    current_state = Samsara.current_state

    case = Samsara.CircularSys( nodes=(:s1, :s2, :s3), actions = (:auke, :minke, :drikke_kaffe) )
    @test case._all_actions == (:auke, :minke, :drikke_kaffe)
    " The action set, i.e. what represents the 3 actions [:up, :down, nothing] can be set by ctor arg "

    #too_few_actions = (:første, :andre)
    #@test_throws ArgumentError Samsara.CircularSys( nodes=(:s1, :s2, :s3), actions = too_few_actions)
    #too_many_actions = (:første, :andre, :a3, :4, :a5)
    #@test_throws ArgumentError Samsara.CircularSys( nodes=(:s1, :s2, :s3), actions = too_many_actions)
    " The set of actions have to be exactly 2: representing :up, :down, and :notning: "

    case = Samsara.CircularSys( nodes=(:s1, :s2, :s3), actions = ("auke", :minke, 42) )
    set_action_in!(case, "auke")
    current_action(case) == "auke"
    @assert current_state(case) == :s1
    @test step_system_mechanics!(case) == :s2
    @assert current_state(case) == :s2
    set_action_in!(case, :minke)
    @test step_system_mechanics!(case) == :s1
    #set_action_in!(case, "Trallalalala, drikke kaffe!")  # (ikkje lista som gyldige actions)
    #@test step_system_mechanics!(case) == :s1
    " Actions kan være kva faen det vil, så lenge det følger (opp, ned, ingenting)-konvensjonen ! "
end

end#module TEST_CIRULARSYS

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
using LAR
include("src/inters.jl")
include("test/boundary.jl")

function boxes3lar(boxes)
	V = Array{Float64,1}[]
	EV = Array{Int64,1}[]
	FV = Array{Int64,1}[]
	for k=1:size(boxes,2)
		xm,ym,zm,xM,yM,zM = boxes[:,k]
		push!(V,[xm,ym,zm])
		push!(V,[xM,yM,zm])
		push!(V,[xm,yM,zm])
		push!(V,[xM,ym,zm])
		push!(V,[xm,ym,zM])
		push!(V,[xM,yM,zM])
		push!(V,[xm,yM,zM])
		push!(V,[xM,ym,zM])

		push!(EV,[ 8(k-1)+1, 8(k-1)+3 ])
		push!(EV,[ 8(k-1)+1, 8(k-1)+4 ])
		push!(EV,[ 8(k-1)+2, 8(k-1)+3 ])
		push!(EV,[ 8(k-1)+2, 8(k-1)+4 ])
		push!(EV,[ 8(k-1)+4+1, 8(k-1)+4+3 ])
		push!(EV,[ 8(k-1)+4+1, 8(k-1)+4+4 ])
		push!(EV,[ 8(k-1)+4+2, 8(k-1)+4+3 ])
		push!(EV,[ 8(k-1)+4+2, 8(k-1)+4+4 ])
		push!(EV,[ 8(k-1)+1, 8(k-1)+4+1 ])
		push!(EV,[ 8(k-1)+2, 8(k-1)+4+2 ])
		push!(EV,[ 8(k-1)+3, 8(k-1)+4+3 ])
		push!(EV,[ 8(k-1)+4, 8(k-1)+4+4 ])

		push!(FV,[ 8(k-1)+1, 8(k-1)+2, 8(k-1)+3, 8(k-1)+4])
		push!(FV,[ 8(k-1)+4+1, 8(k-1)+4+2, 8(k-1)+4+3, 8(k-1)+4+4])
		push!(FV,[ 8(k-1)+1, 8(k-1)+3, 8(k-1)+4+1, 8(k-1)+4+3])
		push!(FV,[ 8(k-1)+1, 8(k-1)+4, 8(k-1)+4+1, 8(k-1)+4+4])
		push!(FV,[ 8(k-1)+2, 8(k-1)+3, 8(k-1)+4+2, 8(k-1)+4+3])
		push!(FV,[ 8(k-1)+2, 8(k-1)+4, 8(k-1)+4+2, 8(k-1)+4+4])
	end
	return hcat(V...), hcat(EV...), hcat(FV...)
end

function boxes2lar(boxes)
	V = Array{Float64,1}[]
	EV = Array{Int64,1}[]
	FV = Array{Int64,1}[]
	for k=1:size(boxes,2)
		xm,ym,xM,yM = boxes[:,k]
		push!(V,[xm,ym])
		push!(V,[xM,yM])
		push!(V,[xm,yM])
		push!(V,[xM,ym])
		push!(EV,[ 4(k-1)+1, 4(k-1)+3 ])
		push!(EV,[ 4(k-1)+1, 4(k-1)+4 ])
		push!(EV,[ 4(k-1)+2, 4(k-1)+3 ])
		push!(EV,[ 4(k-1)+2, 4(k-1)+4 ])
		push!(FV,[ 4(k-1)+1, 4(k-1)+2, 4(k-1)+3, 4(k-1)+4])
	end
	return hcat(V...), hcat(EV...), hcat(FV...)
end

function splitting(boxes,nboxes,ncoords,bucket, finalBuckets,splittingStack,batch=14)
	k, theBucket = bucket
	k += 1
	k = k % ncoords == 0 ? ncoords : k % ncoords
	threshold = ( sum(boxes[k,:])+sum(boxes[k+ncoords,:]) )/(2*nboxes)
	#@show threshold
	pre,into,post = [],[],[]
	for h in theBucket
		themin, themax = boxes[k,h], boxes[k+ncoords,h]
		if themin <= threshold <= themax
			push!(into, h)
		elseif themax < threshold
			push!(pre, h)
		elseif themin > threshold
			push!(post, h)
		end
	end
	current = Set(theBucket)
	#@show p.AA(length)([pre,into,post])

	appended = 0
	if length(pre) <= batch
		push!(pre, copy(into))
		append!(finalBuckets, pre)
		appended += 1
	end
	if length(post) <= batch
		push!(post, copy(into))
		append!(finalBuckets, post)
		appended += 2
	end
	if appended==1
		append!(post, copy(into))
		newblock = (k,post)
		if noloop(splittingStack, ncoords, newblock,current)
			if newblock <= batch
				push!(splittingStack, newblock)
			else
				push!(finalBuckets, newblock[2])
			end
		else
			push!(finalBuckets, newblock[2])
		end
		return
	elseif appended==2
		append!(pre, copy(into))
		newblock = (k,pre)
		if noloop(splittingStack, ncoords, newblock,current)
			if newblock <= batch
				push!(splittingStack, newblock)
			else
				push!(finalBuckets, newblock[2])
			end
		else
			push!(finalBuckets, newblock[2])
		end
		return
	elseif appended==3
		return
	elseif appended==0
		append!(pre, copy(into))
		push!(splittingStack, (k,pre))
		append!(post, copy(into))
		push!(splittingStack, (k,post))
	end
end

function noloop(splittingStack, ncoords, newblock,current)
	j,newterm = 0,Set(newblock[2])
	if newterm == current
		return false
	end
	len = length(splittingStack)
	while j < ncoords
		stored = Set(splittingStack[len-j][2])
		if stored == newterm
			return false
		end
		j += 1
	end
	return true
end

function boxIndexing(boxes)
	m = size(boxes,1)
	nboxes = size(boxes,2)
	index = [k for k=1:nboxes]
	ncoords = Int(m/2)

    bucket = (2,collect(1:size(boxes,2)))
    splittingStack = Tuple{Int64,Array{Any,1}}[bucket]
    finalBuckets = Any[]
    while splittingStack != []
        bucket = pop!(splittingStack)
		splitting(boxes,nboxes,ncoords,bucket, finalBuckets,splittingStack)
	end
	return finalBuckets
end

function boxBuckets2(boxes)
	boxinstances = [ Set{Int64}([]) for k=1:size(boxes,2) ]
	finalBuckets = boxIndexing(boxes)
	for bucket in finalBuckets
		println(bucket)
		for b in bucket
			println(b)
			boxinstances[b] = union(boxinstances[b], Set{Int64}(bucket))
		end
	end
	return [sort(collect(item)) for item in boxinstances]
end

bb = boxBuckets2(boxes)

finalBuckets = boxIndexing(boxes)

boxgroups = [hcat([boxes[:,b] for b in bucket]...) for bucket in finalBuckets]
colors = [p.RED,p.GREEN,p.BLUE,p.CYAN,p.MAGENTA,p.YELLOW,p.BLACK,p.WHITE,p.GRAY,p.BROWN,p.ORANGE,p.PURPLE]
p.VIEW(p.STRUCT([ p.COLOR(colors[k%12+1])(lar2hpc(boxes3lar(group)[1:2]...)) for (k,group) in enumerate(boxgroups) ]))

pairs = collect(zip(boxBuckets(boxes),boxBuckets2(boxes)))

p.AA(p.AA(length))(pairs)

for pair in pairs println(p.AA(length)(pair)) end


text:
## ------------------------------------------------------------------
let
    @info("Testing PropertyChangedEvent")
    global __glob = rand()
    e = PropertyChangedEvent(Main)
    update!(e, :__glob)

    events_count = 0
    trigger_at = [3, 4, 5]
    for it in 1:10

        if has_event!(e, :__glob)
            println(:__glob, " changed!!!")
            events_count += 1
        end

        if it in trigger_at
            __glob = rand()
        end
        sleep(0.1)
    end
    @test events_count == length(trigger_at)
end


text:
using DiffEqBayes

using Distributions
using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions, DiffEqUncertainity
using Plots
using DiffEqMonteCarlo

gr(fmt=:png)

f = @ode_def LotkaVolterraTest begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d

u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1,0]

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(1,stop=10,length=10))
sig = 0.49
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))

scatter(t, data[1,:], lab="#prey (data)")
scatter!(t, data[2,:], lab="#predator (data)")
plot!(sol)

priors = [Truncated(Normal(1.5,0.5),0.5,2.5),Truncated(Normal(1.2,0.5),0,2),Truncated(Normal(3.0,0.5),1,4),Truncated(Normal(1.0,0.5),0,2)]

cb = AdaptiveProbIntsUncertainty(5)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=100,callback=cb,reltol=1e-5,abstol=1e-5)
plot(sim,vars=(0,1),linealpha=0.4)

@time bayesian_result_stan = stan_inference(prob,t,data,priors;reltol=1e-5,abstol=1e-5,vars =(StanODEData(),InverseGamma(3,2)))

plot_chain(bayesian_result_stan)

@time bayesian_result_turing = turing_inference(prob,Tsit5(),t,data,priors)

plot_chain(bayesian_result_turing)


text:
const unary_operators = [:(+), :(-), :(!), :(*)]

const numeric_unary_operators = [:(+), :(-)]

const logical_unary_operators = [:(!)]

const elementary_functions = [:abs, :sign, :acos, :acosh, :asin,
                              :asinh, :atan, :atanh, :sin, :sinh,
                              :cos, :cosh, :tan, :tanh, :ceil, :floor,
                              :round, :trunc, :exp, :exp2, :expm1, :log, :log10, :log1p,
                              :log2, :exponent, :sqrt, :gamma, :lgamma, :digamma,
                              :erf, :erfc]

const two_argument_elementary_functions = [:round, :ceil, :floor, :trunc]

const special_comparison_operators = [:isless]

const comparison_operators = [:(==), :(.==), :(!=), :(.!=),
                              :(>), :(.>), :(>=), :(.>=), :(<), :(.<),
                              :(<=), :(.<=)]

const scalar_comparison_operators = [:(==), :(!=), :(>), :(>=),
                                     :(<), :(<=)]

const array_comparison_operators = [:(.==), :(.!=), :(.>), :(.>=), :(.<), :(.<=)]

const vectorized_comparison_operators = [(:(.==), :(==)), (:(.!=), :(!=)),
                                         (:(.>), :(>)), (:(.>=), :(>=)),
                                         (:(.<), :(<)), (:(.<=), :(<=))]

const binary_operators = [:(+), :(.+), :(-), :(.-), :(*), :(.*), :(/), :(./),
                          :(.^), :(div), :(mod), :(fld), :(rem)]

const induced_binary_operators = [:(^)]

const arithmetic_operators = [:(+), :(.+), :(-), :(.-), :(*), :(.*), :(/), :(./),
                              :(.^), :(div), :(mod), :(fld), :(rem)]

const induced_arithmetic_operators = [:(^)]

const biscalar_operators = [:(max), :(min)]

const scalar_arithmetic_operators = [:(+), :(-), :(*), :(/),
                                     :(div), :(mod), :(fld), :(rem)]

const induced_scalar_arithmetic_operators = [:(^)]

const array_arithmetic_operators = [:(+), :(.+), :(-), :(.-), :(.*), :(.^)]

const bit_operators = [:(&), :(|), :($)]

const unary_vector_operators = [:minimum, :maximum, :prod, :sum, :mean, :median, :std,
                                :var, :mad, :norm, :skewness, :kurtosis]

# TODO: dist, iqr, rle, inverse_rle

const pairwise_vector_operators = [:diff, :reldiff]#, :percent_change]

const cumulative_vector_operators = [:cumprod, :cumsum, :cumsum_kbn, :cummin, :cummax]

const ffts = [:fft]

const binary_vector_operators = [:dot, :cor, :cov, :cor_spearman]

const rowwise_operators = [:rowmins, :rowmaxs, :rowprods, :rowsums,
                           :rowmeans, :rowmedians, :rowstds, :rowvars,
                           :rowffts, :rownorms]

const columnar_operators = [:colmins, :colmaxs, :colprods, :colsums,
                            :colmeans, :colmedians, :colstds, :colvars,
                            :colffts, :colnorms]

const boolean_operators = [:any, :all]

# Swap arguments to fname() anywhere in AST. Returns the number of
# arguments swapped
function swapargs(ast::Expr, fname::Symbol)
    if ast.head == :call &&
       (ast.args[1] == fname ||
        (isa(ast.args[1], Expr) && ast.args[1].head == :curly &&
         ast.args[1].args[1] == fname)) &&
       length(ast.args) == 3

        ast.args[2], ast.args[3] = ast.args[3], ast.args[2]
        1
    else
        n = 0
        for arg in ast.args
            n += swapargs(arg, fname)
        end
        n
    end
end
function swapargs(ast, fname::Symbol)
    ast
    0
end

# Return a block consisting of both the given function and a copy of
# the function in which arguments to the function itself and any
# 2-argument calls to a function of the same name are swapped
macro swappable(func, syms...)
    if (func.head != :function && func.head != :(=)) ||
       func.args[1].head != :call || length(func.args[1].args) != 3
        error("@swappable may only be applied to functions of two arguments")
    end

    func2 = deepcopy(func)
    fname = func2.args[1].args[1]
    if isa(fname, Expr)
        if fname.head == :curly
            fname = fname.args[1]
        else
            error("Unexpected function name $fname")
        end
    end

    for s in unique([fname, syms...])
        if swapargs(func2, s) < 1
            error("No argument swapped")
        end
    end
    esc(Expr(:block, func, func2))
end

#
# Unary operator macros for DataFrames and DataArrays
#

macro dataframe_unary(f)
    esc(:($(f)(d::DataFrame) = DataFrame([$(f)(d[i]) for i=1:size(d, 2)], deepcopy(index(d)))))
end

#
# Binary operator macros for DataFrames and DataArrays
#

macro dataframe_binary(f)
    esc(quote
        function $(f)(a::DataFrame, b::DataFrame)
            if size(a) != size(b); error("argument dimensions must match"); end
            DataFrame([$(f)(a[i], b[i]) for i=1:size(a, 2)], deepcopy(index(a)))
        end
        @swappable $(f)(a::DataFrame, b::Union(Number, String)) =
            DataFrame([$(f)(a[i], b) for i=1:size(a, 2)], deepcopy(index(a)))
        @swappable $(f)(a::DataFrame, b::NAtype) =
            DataFrame([$(f)(a[i], b) for i=1:size(a, 2)], deepcopy(index(a)))
    end)
end

# Unary operators, DataFrames
@dataframe_unary !
@dataframe_unary -
# As in Base, these are identity operators
for f in (:(+), :(*))
    @eval $(f)(d::DataFrame) = d
end

#
# Elementary functions
#
# XXX: The below should be revisited once we have a way to infer what
# the proper return type of an array should be.

# One-argument elementary functions that return the same type as their
# inputs
for f in (:abs, :sign)
    @eval begin
        @dataframe_unary $(f)
    end
end

# One-argument elementary functions that always return floating points
for f in (:acos, :acosh, :asin, :asinh, :atan, :atanh, :sin, :sinh, :cos,
          :cosh, :tan, :tanh, :exp, :exp2, :expm1, :log, :log10, :log1p,
          :log2, :exponent, :sqrt, :gamma, :lgamma, :digamma, :erf, :erfc)
    @eval begin
        @dataframe_unary $(f)
    end
end

# Elementary functions that take varargs
for f in (:round, :ceil, :floor, :trunc)
    @eval begin
        $(f)(d::DataFrame, args::Integer...) =
            DataFrame([$(f)(d[i], args...) for i=1:size(d, 2)], deepcopy(index(d)))
    end
end

#
# Bit operators
#

for f in (:&, :|, :$)
    @eval begin
        # DataFrame
        @dataframe_binary $(f)
    end
end

function isequal(df1::AbstractDataFrame, df2::AbstractDataFrame)
    if size(df1, 2) != size(df2, 2)
        return false
    end
    for idx in 1:size(df1, 2)
        if !isequal(df1[idx], df2[idx])
            return false
        end
    end
    return true
end

for sf in scalar_comparison_operators
    vf = symbol(".$sf")
    @eval begin
        @dataframe_binary $vf
    end
end

#
# Binary operators
#

for f in arithmetic_operators
    @eval begin
        @dataframe_binary $f
    end
end

for f in (:minimum, :maximum, :prod, :sum, :mean, :median, :std, :var, :norm)
    colf = symbol("col$(f)s")
    rowf = symbol("row$(f)s")
    @eval begin
        function ($colf)(df::AbstractDataFrame)
            p = ncol(df)
            res = DataFrame()
            for j in 1:p
                res[j] = DataArray(($f)(df[j]))
            end
            names!(res, names(df))
            return res
        end
    end
end

#
# Boolean operators
#

function all(df::AbstractDataFrame)
    for i in 1:size(df, 2)
        x = all(df[i])
        if isna(x)
            return NA
        end
        if !x
            return false
        end
    end
    true
end

function any(df::AbstractDataFrame)
    has_na = false
    for i in 1:size(df, 2)
        x = any(df[i])
        if !isna(x)
            if x
                return true
            end
        else
            has_na = true
        end
    end
    has_na ? NA : false
end


text:
abstract SeqRead

type FastaRead <: SeqRead
	name::ASCIIString
	sequence::Sequence
end

type FastqRead <: SeqRead
	name::ASCIIString
	sequence::Sequence
	strand::ASCIIString
	quality::Quality
end

"""
read1 and read2 of a FastqPair are considered with different directions
"""
type FastqPair
	read1::FastqRead
	read2::FastqRead

	# template length, if it's not calculated, it is set as -1
	template_length::Int
	# overlap length, if it's not calculated, it is set as -1
	overlap_length::Int
end

function FastqPair(read1::FastqRead, read2::FastqRead)
	return FastqPair(read1, read2, -1, -1)
end

Base.length(read::SeqRead) = length(read.sequence)

Base.getindex(read::FastqRead, i::Int64) = (getindex(read.sequence, i), getindex(read.quality, i))
Base.getindex(read::FastqRead, r::UnitRange{Int64}) = FastqRead(read.name, getindex(read.sequence, r), read.strand, getindex(read.quality, r))
Base.getindex(read::FastqRead, indx::AbstractArray{Int64,1}) = FastqRead(read.name, getindex(read.sequence, indx), read.strand, getindex(read.quality, indx))
Base.reverse(read::FastqRead) = FastqRead(read.name, reverse(read.sequence), read.strand, reverse(read.quality))
complement(read::FastqRead) = FastqRead(read.name, complement(read.sequence), read.strand, read.quality)
reverse_complement(read::FastqRead) = FastqRead(read.name, reverse_complement(read.sequence), read.strand, reverse(read.quality))
Base.display(read::FastqRead) = display_fastq(read)

Base.getindex(read::FastaRead, i::Int64) = getindex(read.sequence, i)
Base.getindex(read::FastaRead, r::UnitRange{Int64}) = FastaRead(read.name, getindex(read.sequence, r))
Base.getindex(read::FastaRead, indx::AbstractArray{Int64,1}) = FastaRead(read.name, getindex(read.sequence, indx))
Base.reverse(read::FastaRead) = FastaRead(read.name, reverse(read.sequence))
complement(read::FastaRead) = FastaRead(read.name, complement(read.sequence))
reverse_complement(read::FastaRead) = FastaRead(read.name, reverse_complement(read.sequence))
Base.display(read::FastaRead) = display_fasta(read)

Base.display(pair::FastqPair) = display_fastq_pair(pair)

-(read::FastqRead) = reverse(read)
!(read::FastqRead) = complement(read)
~(read::FastqRead) = reverse_complement(read)
-(read::FastaRead) = reverse(read)
!(read::FastaRead) = complement(read)
~(read::FastaRead) = reverse_complement(read)

function display_fasta(read::FastaRead)
	println(">" * read.name)
	display_sequence(read.sequence)
end

function display_fastq(read::FastqRead)
	println(read.name)
	println(read.sequence.seq)
	println(read.strand)
	println(read.quality.qual)
end

function display_fastq_pair(pair::FastqPair)
	println("----Read #1----")
	display_fastq(pair.read1)
	println("----Read #2----")
	display_fastq(pair.read2)
end


text:
get_language(config::I18nConfig = GlobalI18nConfig) = config.current_language
get_language(ctx::I18nContext) = get_language(get_config(ctx))

function set_language!(
    config::I18nConfig,
    locale_name::AbstractString = get_system_language(),
    fallback::Union{AbstractString, AbstractVector{<:AbstractString}} = GlobalI18nConfig.fallback
)
    if fallback isa AbstractString
        fallback = [fallback]
    end
    previous = config.current_language
    locale_name = parse_locale_name(locale_name)
    config.current_language = locale_name
    config.fallback = parse_locale_name.(fallback)
    for func in get(OnLanguageChange, config, [])
        func(locale_name, previous)
    end
    config
end

set_language(
    locale_name::AbstractString = get_system_language(),
    fallback::Union{AbstractString, AbstractVector{<:AbstractString}} = GlobalI18nConfig.fallback
) = set_language!(GlobalI18nConfig, locale_name, fallback)

function set_language!(
    ctx::I18nContext,
    locale_name::Union{AbstractString, Nothing} = nothing,
    fallback::Union{AbstractString, AbstractVector{<:AbstractString}} = get_config(ctx).fallback
)
    config = get_config(ctx)
    if isnothing(locale_name)
        isnothing(config) && return ctx
        ctx.override_config = nothing
        pop!(OnLanguageChange, config, nothing)
        return ctx
    end
    if isnothing(config)
        config = ctx.config = I18nConfig()
    end
    set_language!(config, locale_name, fallback)
end

function get_entry(data::Union{I18nData, String}, path)
    for x in path
        data isa I18nData && haskey(data, x) || return nothing
        data = data[x]
    end
    data isa String ? data : nothing
end

function try_languages(ctx::I18nContext, path, languages)
    data = ctx.data
    for lang in languages
        haskey(data, lang) || continue
        entry = get_entry(data[lang], path)
        isnothing(entry) || return entry
    end
end

function i18n(ctx::I18nContext, key, default = key; language = nothing)
    language = if isnothing(language)
        get_language(ctx)
    else
        parse_locale_name(language)
    end
    path = split(key, '.')
    entry = try_languages(ctx, path, [language])
    isnothing(entry) || return entry
    entry = try_languages(ctx, path, ctx.fallback(language))
    isnothing(entry) ? default : entry
end

function i18n(current_module::Module, key, default = key; language = nothing)
    ctx = get(I18nContexts, current_module, nothing)
    isnothing(ctx) && return default
    i18n(ctx, key, default, language = language)
end

@noinline function i18n(key, default = key; language = nothing)
    current_module = @caller_module(5)
    i18n(current_module, key, default, language = language)
end

macro i_str(s)
    :(i18n(@caller_module(2), $s))
end


text:
function process_session(filename::String)
    ongoing = MAT.matread(filename)
    MouseID = match(r"[a,b,c,d,p]{1,2}\d{1,2}",filename).match #ongoing["saved_history"]["SavingSection_ratname"][3]
    preday = match(r"\d{6}",filename).match
    println(preday," ",MouseID)
    Day = Date("20"*replace(preday,r"a.mat" => ""),"yyyymmdd")
    Session = MouseID * "_" * string(Day)
    Phase = get(Phase_Calendar,Day,"training")
    if string(MouseID[1]) == "p"
        Group = string(MouseID[end]) in ["1","3","5","7","9"] ? "Group C" : "Group D"
    else
        Group = string(MouseID[1]) in ["a","d"] ? "Group A" : "Group B"
    end
    Treatment = get_treatment(Phase,Group,Day) #get_injection(Phase,Group,Day)
    Injection = occursin("Veh",Treatment) ? "Vehicle" : Treatment
    ExpType = occursin("pc",MouseID) ? "Optogenetic" : "Pharmacology"
    pokes = process_pokes(ongoing)
    sort!(pokes,:PokeIn)
    streaks = process_streaks(pokes)
    for df in (pokes,streaks)
        for (n,v) in zip([:MouseID,:Day,:Phase,:Group,:Treatment, :Injection, :ExpType],[MouseID,Day,Phase,Group,Treatment, Injection, ExpType])
            df[!,n] .= v
        end
    end
    return pokes, streaks
end

function process_all()
    source = files_dict()
    full_pokes = DataFrame()
    full_streaks = DataFrame()

    for idx in 1:length(source["Files"])
        filepath = get_rawfile(idx)
        ongoing_pokes, ongoing_streaks = process_session(filepath)
        println(filepath)
        if isempty(full_pokes)
            full_pokes = ongoing_pokes
            full_streaks = ongoing_streaks
        else
            append!(full_pokes,ongoing_pokes)
            append!(full_streaks,ongoing_streaks)
        end
    end
    full_pokes[!,:ExpDay] = exp_calendar(full_pokes)
    full_pokes[!,:ProtocolDay] = exp_calendar(full_pokes,:Phase)
    full_streaks[!,:ExpDay] = exp_calendar(full_streaks)
    full_streaks[!,:ProtocolDay] = exp_calendar(full_streaks,:Phase)
    return full_pokes,full_streaks
end


text:
function tunningexpand(ex)
end

function tunninginput(ex)
end

function tunningvar(i)
end

shouldexpand(ex) = (isa(ex, Expr) && (ex.head == :call) && (ex.args[1] == :tunningexpand))

isinput(ex) = (isa(ex, Expr) && (ex.head == :call) && (ex.args[1] == :tunninginput))

isvar(ex) = (isa(ex, Expr) && (ex.head == :call) && (ex.args[1] == :tunningvar))

"""
    getsymbols!(ex, s)

This functions search recursively in `ex` by variables that are marked
with `tunninginput`.

It returns a new expression with this mark removed and also append to
Set `s` the name of the variables found during the search.

"""
function getsymbols!(ex, s::Set{Symbol})

    # If is the `isinput` command, remove it and return only the
    # argument
    if isinput(ex)

        push!(s, ex.args[2])

        return ex.args[2]

    end

    # If is the `tunningvar` command, add :__x to the set of symbols
    # (as usual), but associate with this variable the provided
    # argument.
    if isvar(ex)

        push!(s, :__x)

        return Expr(:ref, :__x, ex.args[2])

    end

    isa(ex, Expr) &&
        (return Expr(ex.head, (getsymbols!(exx, s) for exx in ex.args)...))

    return ex

end

"""
    expand_args(args, pn=0)

This function receives a list `args` or arguments to a function call
and, possibly, the number `pn` of optimization variables known before
a call to this function. It recursively walks through the AST of each
expression in the arguments and perform some actions

  - If symbol `:_` is found, then this argument is a new optimization
    argument

  - If a call to `tunningexpand` is performed (usually passing an
    Iterable as argument), then the call is removed and the associated
    expression is replicated, replacing the argument by each value
    from the evaluation of the argument of `tunningexpand`

  - If a call to `tunninginput` is found, then a new user-needed
    variable is found, but this variable is not an optimization
    variable. Simply remove the call and add the variable to the list
    of symbols.

It returns a `list` where, at each position there is a proper list of
the arguments to the function call. It also returns the number `n` of
optimization variables and a `Set` with all the symbols (including the
optimization variable) that will be needed to call the function.

"""
function expand_args(args, pn=0)

    nargs = length(args)

    args_list = [Vector{Any}(undef, nargs)]

    n = pn

    symbols = Set{Symbol}()

    for (i, ex) in enumerate(args)

        if ex == :_

            # In this case, we add a new variable for the optimization
            # problem. Unlike the `tunningvar(i)` option, in this
            # case, we are free to set the order of variables

            # Check if there is no use of :_ and :tunningvar
            (n == 0) && (:__x in symbols) && error("It is forbidden to use :_ and :tunningvar in the same call")

            n += 1

            push!(symbols, :__x)

            for vargs in args_list

                # It was _, now it is __x[n]
                vargs[i] = Expr(:ref, :__x, n)

            end

        elseif isa(ex, Expr) && (ex.head in [:kw, :parameters] || shouldexpand(ex))

            if shouldexpand(ex)

                # If it was indicated that the current expression
                # should be expanded to generate different function
                # calls

                isfirst = true

                local tmplist = []

                for j in eval(ex.args[2])

                    for vargs in args_list
                        tmpargs = (isfirst) ? vargs : copy(vargs)
                        tmpargs[i] = j
                        (!isfirst) && (push!(tmplist, tmpargs))
                    end
                    # After the first loop, need to create copies
                    isfirst = false

                end

                append!(args_list, tmplist)

            elseif ex.head == :kw

                local tmplist = []

                isfirst = true

                # When dealing with keywords, only have to expand the
                # second argument, if necessary
                expanded_args, n, new_symbols = expand_args(ex.args[2:end], n)

                union!(symbols, new_symbols)

                for j in expanded_args

                    for vargs in args_list
                        tmpargs = (isfirst) ? vargs : copy(vargs)
                        tmpargs[i] = Expr(:kw, ex.args[1], j...)
                        (!isfirst) && (push!(tmplist, tmpargs))
                    end
                    # After the first loop, need to create copies
                    isfirst = false

                end

                append!(args_list, tmplist)

            elseif ex.head == :parameters

                local tmplist = []

                isfirst = true

                # Recursively expands all the keywords associated with
                # the parameters and update the number of variables
                expanded_args, n, new_symbols = expand_args(ex.args, n)

                union!(symbols, new_symbols)

                for j in expanded_args

                    for vargs in args_list
                        tmpargs = (isfirst) ? vargs : copy(vargs)
                        # Here we only have to add a new expression
                        # and "unroll" the arguments
                        tmpargs[i] = Expr(:parameters, j...)
                        (!isfirst) && (push!(tmplist, tmpargs))
                    end
                    # After the first loop, need to create copies
                    isfirst = false

                end

                append!(args_list, tmplist)

            end

        else

            # Default case: simply add the current expression in all
            # vectors, after removing the symbols

            new_symbols = Set{Symbol}()

            exx = getsymbols!(ex, new_symbols)

            # Check if both :_ and :tunningvar are being used
            (n > 0) && (:__x in new_symbols) &&
                error("It is forbidden to use :_ and :tunningvar in the same call")

            union!(symbols, new_symbols)

            for vargs in args_list

                vargs[i] = exx

            end

        end

    end

    return args_list, n, symbols

end

macro prepare_tunning(f_call, comparefunction)

    (f_call.head != :call) && error("Must be a function call.")

    # Expand, if necessary, all the arguments after the function name
    expanded_args, n, symbols = expand_args(f_call.args[2:end])

    # Add the function call to each different configuration of the
    # algorithm
    fexprs = map((args) -> Expr(:call, f_call.args[1], args...),
                 expanded_args)

    # Now create the functions to measure each different configuration
    # of the algorithm call

    fex = Vector{Any}(undef, length(fexprs))
    flist = Vector{Any}(undef, length(fexprs))

    for (i, ex) in enumerate(fexprs)

        local fi = Symbol("f", i)
        fex[i] = quote
            function $(fi)(__x)
                # Call the solver in a series of problems
                ff($(symbols...)) = $ex
                return $comparefunction(__x, ff)
            end
        end

        flist[i] = :($(fi))
    end

    # Create the functions and the list

    return Expr(:block, fex..., Expr(:vect, flist...))

end

function rcopycode(fname, code)

    !isa(code, Expr) && return code, []

    # If found the given function, then create copies, etc...
    if (code.head == :call) && (code.args[1] == fname)

        # Expand, if necessary, all the arguments after the function name
        expanded_args, n, symbols = expand_args(code.args[2:end])

        # Add the function call to each different configuration of the
        # algorithm
        fexprs = map((args) -> Expr(:call, code.args[1], args...),
                     expanded_args)

        return code, fexprs

    end

    # Otherwise, try to find it further
    for (i, ex) in enumerate(code.args)

        c, l = rcopycode(fname, ex)

        !isempty(l) && return c, l

    end

    return code, []

end

macro copycode(fname, code)

    # Something like this
    target, copies = rcopycode(fname, code)

    # Now create the functions to measure each different configuration
    # of the algorithm call

    fex = Vector{Any}(undef, length(copies))
    flist = Vector{Any}(undef, length(copies))

    for (i, ex) in enumerate(copies)

        # Put the correct arguments
        target.args .= ex.args

        local fi = Symbol("f", i)
        fex[i] = quote
            function $(fi)(__x)
                # Add the code
                $(copy(code))
            end
        end

        flist[i] = :($(fi))
    end

    # Create the functions and the list

    return Expr(:block, fex..., Expr(:vect, flist...))

end


text:
###########
# Effects #
###########

abstract type Effect end # What's the virtue of having an abstract type for this?

struct NoEffect <: Effect end
struct Gravity <: Effect; specific_force::Vector{Float64}; end
struct BodyForce <: Effect; force::Vector{Float64}; position::Vector{Float64}; end
struct BodyTorque <: Effect; torque::Vector{Float64}; end

mutable struct BodyStateEffect{T <: Real} <: Effect
    r_be_I::Vector{T} # Position of center of mass wrt center of mass of earth in ICRF's coordinates
    v_be_I::Vector{T} # Rate of change of above (as viewed from the inertia frame)
    q_BI::Vector{T}   # Attitude quaternion of body frame wrt ICRF
    ω_BI_B::Vector{T} # Rotation rate of body frame wrt ICRF, in body coordinates
end
BodyStateEffect() = BodyStateEffect(zeros(3), zeros(3), [0.; 0.; 0.; 1.], zeros(3))

#############
# Functions #
#############

# Find a single instance of an effect in a dictionary of effects.
# If the dictionary of effects were always the same, then we could
# make a _generated_ function to just return the one we want. That would
# be super fast. In fact, it would work even if the list changed sometimes,
# because we'd just get a new generated function, but that would cause weird
# delays at runtime, and it would be slow if the list were frequently different.
function find_effect(effects_bus::Dict{String, Tuple}, null_effect::T) where T
    for effects in values(effects_bus) # effects contains tuple of Effects from each component
        index = findfirst(e->isa(e, T), effects)
        if !isnothing(index)
            return (effects[index], true)
        end
    end
    return (null_effect, false)
end

# Find a single instance of an effect in a vector of effects.
function find_effect(effects::Vector, null_effect::T) where T
    index = findfirst(e->isa(e, T), effects)
    if !isnothing(index)
        return (effects[index], index)
    end
    return (null_effect, 0)
end

# function find_effect(all_effects::Dict{String, Tuple}, effect_type::Type)
#     for effects in values(all_effects)
#         index = findfirst(e->isa(e, effect_type), effects)
#         if !iszero(index)
#             return effects[index]
#         end
#     end
#     return NoEffect()
# end
#
# function find_effects(all_effects::Dict{String, Tuple}, effect_type::Type)
#     for effects in values(all_effects)
#         index = findfirst(e->isa(e, effect_type), effects)
#         if !iszero(index)
#             return effects[index]
#         end
#     end
#     return NoEffect()
# end

# # Assumes effects is an iterable array and that the + method exist for T.
# function accumulate(x::T, effects) where T
#     for effect in Iterators.filter(e -> isa(e, T), effects)
#         x += effect
#     end
#     return x
# end
#
# # Allows a user to define a function to accumulate the type into another type
# # (e.g., to extract a vector).
# function accumulate(x, effects, effect_type::Type, f = e->e, args...)
#     for effect in Iterators.filter(e -> isa(e, effect_type), effects)
#         x += f(effect, args...)
#     end
#     return x
# end
#
# # Acts in place and prevents unnecessary creation of vectors.
# function accumulate!(x::Vector, effects, effect_type::Type, f = e->e, args...)
#     for effect in Iterators.filter(e -> isa(e, effect_type), effects)
#         x[:] += f(effect, args...)
#     end
# end

# Acts in place and prevents unnecessary creation of vectors.
# We could also make this a generated function.
function accumulate!(x::Vector, effects_bus::Dict, effect_type::Type, f = e->e, args...)
    for effects in values(effects_bus)
        for effect in effects
            if isa(effect, effect_type)
                x[:] += f(effect, args...)
            end
        end
    end
end

# function accumulate(x, effects_bus::Dict{String, Tuple}, effect_type::Type, f = e->e, args...)
#     return accumulate(x, flatten(values(effects_bus)), effect_type, f, args...)
# end


text:
# # Quantum Circuit Born Machine
using Yao, YaoExtensions
import Yao: probs
using QuAlgorithmZoo: Adam, update!

struct QCBM{BT<:AbstractBlock, MT<:MMD}
    circuit::BT
    mmd::MT
end

"""generated output probability distribution"""
function probs(qcbm::QCBM)
    zero_state(qcbm.circuit |> nqubits) |> qcbm.circuit |> probs
end

function loss(qcbm::QCBM)
    expect(qcbm.mmd, probs(qcbm) |> as_weights)
end

function getgrad(qcbm::QCBM)
    expect'(qcbm.mmd, zero_state(nqubits(qcbm.circuit))=>qcbm.circuit).second
end

# ## DATA: Target Probability Distribution
# The gaussian probability disctribution in phase space of 2^6
nbit = 6
N = 1<<nbit

function gaussian_pdf(x, μ::Real, σ::Real)
    pl = @. 1 / sqrt(2pi * σ^2) * exp(-(x - μ)^2 / (2 * σ^2))
    pl / sum(pl)
end
pg = gaussian_pdf(1:N, N/2-0.5, N/4);

# ## MODEL: Quantum Circuit and Loss
# Using a random differentiable circuit of depth 6 for training, the kernel function is universal RBF kernel
depth = 6
kernel = rbf_kernel(0.25)
c = variational_circuit(nbit, depth, pair_ring(nbit))
dispatch!(c, :random)
qcbm = QCBM(c, MMD(kernel, pg))

# ## TRAINING: Adam Optimizer
# We probide the QCBMGo! iterative interface for training
niter = 100
optim = Adam(lr=0.1)

params = parameters(qcbm.circuit)
for i=1:niter
    # initialize the parameters
    update!(params, getgrad(qcbm), optim)
    dispatch!(qcbm.circuit, params)
    println("Step = $i, Loss = $(loss(qcbm))")
end


text:
module Fretboards

using ArgCheck: @argcheck
using Crayons: Crayon
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using OffsetArrays: OffsetMatrix
using UnPack: @unpack

include("pitch.jl")
include("fretboard.jl")

end # module


text:
using Colors

abstract type ColoredObject <: Object end
abstract type Pigment <: JR end

export FT
struct FT <: Pigment
    filter::Float64
    transmit::Float64

    function FT(filter,transmit)
        if filter<0 || transmit<0 || filter+transmit>1
            error("だめです")
        else
            new(filter,transmit)
        end
    end
end

export rgbColor
struct rgbColor <: ColoredObject
    object :: Object
    color :: Color

    function rgbColor(object,color)
        if object == Empty()
            Empty()
        else
            new(object,color)
        end
    end
end

export rgbftColor
struct rgbftColor <:ColoredObject
    object :: Object
    color :: Color
    transparence :: FT

    function rgbftColor(object,color,ft)
        if object == Empty()
            Empty()
        else
            new(object,color,ft)
        end
    end
end

export rgbftColor2
struct rgbftColor2 <:ColoredObject
    object :: Object
    color :: Color
    transparence :: FT

    function rgbftColor2(object,color,ft)
        if object == Empty()
            Empty()
        else
            new(object,color,ft)
        end
    end
end

# Color
function povray_script(color :: Color)
    r = string(Float64(color.r))
    g = string(Float64(color.g))
    b = string(Float64(color.b))
    return "rgb<"*r*","*g*","*b*">"
end

function povray_script(color :: Color, ft:: FT)
    r = string(Float64(color.r))
    g = string(Float64(color.g))
    b = string(Float64(color.b))
    f = string(Float64(ft.filter))
    t = string(Float64(ft.transmit))
    return "rgbft<"*r*","*g*","*b*","*f*","*t*">"
end

function povray_script(rgbcolor :: rgbColor; preindent = 0)
    script = "object{\n" * "  "^(preindent+1)
    script *= povray_script(rgbcolor.object, preindent=preindent+1) * "\n" * "  "^(preindent+1)
    script *= "pigment{" * povray_script(rgbcolor.color) * "}" * "\n" * "  "^preindent
    script *= "}"
    return script
end

function povray_script(rgbftcolor :: rgbftColor; preindent = 0)
    r = string(Float64(rgbftcolor.color.r))
    g = string(Float64(rgbftcolor.color.g))
    b = string(Float64(rgbftcolor.color.b))
    f = string(Float64(rgbftcolor.transparence.filter))
    t = string(Float64(rgbftcolor.transparence.transmit))

    script = "object{\n" * "  "^(preindent+1)
    script *= povray_script(rgbftcolor.object, preindent=preindent+1) * "\n" * "  "^(preindent+1)
    script *= "pigment{rgbft<"*r*","*g*","*b*","*f*","*t*">}" * "\n" * "  "^preindent
    script *= "}"
    return script
end

function povray_script(rgbftcolor :: rgbftColor2; preindent = 0)
    r = string(Float64(rgbftcolor.color.r))
    g = string(Float64(rgbftcolor.color.g))
    b = string(Float64(rgbftcolor.color.b))
    f = string(Float64(rgbftcolor.transparence.filter))
    t = string(Float64(rgbftcolor.transparence.transmit))

    script = "object{\n" * "  "^(preindent+1)
    script *= povray_script(rgbftcolor.object, preindent=preindent+1) * "\n" * "  "^(preindent+1)
    script *= "texture{pigment{rgbft<"*r*","*g*","*b*","*f*","*t*">}}" * "\n" * "  "^preindent
    script *= "texture{pigment{rgbft<"*r*","*g*","*b*","*f*","*t*">}}" * "\n" * "  "^preindent
    script *= "}"
    return script
end

export Transparent
function Transparent(object::PrimitiveObject,ft::FT)
    return object
end
function Transparent(csg::csgObject,ft::FT)
    csgtype = typeof(csg)
    return csgtype(map(obj->Transparent(obj,ft),csg.objects))
end

function Transparent(rgbcolor::rgbColor,ft::FT)
    return rgbftColor(rgbcolor.object,rgbcolor.color,ft)
end
function Transparent(rgbftcolor::rgbftColor,ft::FT)
    f0 = rgbftcolor.transparence.filter
    t0 = rgbftcolor.transparence.transmit
    f1 = ft.filter
    t1 = ft.transmit
    f2 = 1-(1-f0)*(1-f1)
    t2 = 1-(1-t0)*(1-t1)
    return rgbftColor(rgbftcolor.object,rgbftcolor.color,FT(f2,t2))
end
function Transparent(affinetransform::AffineTransform,ft::FT)
    return AffineTransform(Transparent(affinetransform.object,ft),affinetransform.matrix, affinetransform.vector)
end
function Transparent(paralleltranslation::ParallelTranslation,ft::FT)
    return ParallelTranslation(Transparent(paralleltranslation.object,ft),paralleltranslation.vector)
end
function Transparent(scaling::Scaling,ft::FT)
    return Scaling(Transparent(scaling.object,ft),scaling.scalar)
end
#TODO: Transparent for rgbftColors

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
using aoc2019.sif: image, resolve, checksum

input = joinpath(@__DIR__, "input")
bits = [parse(Int, char) for char in readline(input)]
img = image(bits, 25, 6)

p1 = checksum(img)
@assert p1 == 2806

p2 = resolve(img)

print("-----------------------------------------------------------------------\n")
print("space image format -- part one\n    corruption check : $p1\n")
print("space image format -- part two\n    password         :\n")
for y in 1:6
    print("        ")
    for x in 1:25
        glyph = p2[x, y] == 1 ? "#" : " "
        print("$glyph")
    end
    println()
end
print("-----------------------------------------------------------------------\n")


text:
function asap_ldcs(num_links, A)
	multiple_copies = true
	arch = TopLevel{A,2}("asap_ldcs")

	# Extra parameters
	num_fifos = 2

	# address for long distance circuit switched link processsors
	ldcs_proc_addrs = [(5,8),(15,8),(5,22),(15,22),(5,10),(5,20),(15,10),(15,20)]

    ####################
    # Normal Processor #
    ####################
    processor = build_processor_tile(num_links)
    for r in 1:24, c in 2:28
		in((r,c), ldcs_proc_addrs) && continue
        if multiple_copies
            add_child(arch, deepcopy(processor), CartesianIndex(r,c))
        else
            add_child(arch, processor, CartesianIndex(r,c))
        end
    end
    for r in 1:20, c in 29
		in((r,c), ldcs_proc_addrs) && continue
        if multiple_copies
            add_child(arch, deepcopy(processor), CartesianIndex(r,c))
        else
            add_child(arch, processor, CartesianIndex(r,c))
        end
    end

	#################################################
    # Long Distance Circuit Switched Link Processor #
    #################################################
    processor = build_ldcs_processor_tile(num_links)
    for ldcs_proc_addr in ldcs_proc_addrs
        if multiple_copies
            add_child(arch, deepcopy(processor), CartesianIndex(ldcs_proc_addr[1],ldcs_proc_addr[2]))
        else
            add_child(arch, processor, CartesianIndex(ldcs_proc_addr[1],ldcs_proc_addr[2]))
        end
    end

    ####################
    # Memory Processor #
    ####################
    memory_processor = build_memory_processor_tile(num_links)
    for r = 25, c = 2:28
        if multiple_copies
            add_child(arch, deepcopy(memory_processor), CartesianIndex(r,c))
        else
            add_child(arch, memory_processor, CartesianIndex(r,c))
        end
    end

    #################
    # 2 Port Memory #
    #################
    memory_2port = build_memory_2port()
    for r = 26, c in (2, 4, 6, 8, 11, 13, 15, 17, 20, 22, 24, 26)
        add_child(arch, memory_2port, CartesianIndex(r,c))
    end

    #################
    # 1 Port Memory #
    #################
    memory_1port = build_memory_1port()
    for r = 26, c in (10, 19)
        add_child(arch, memory_1port, CartesianIndex(r,c))
    end

    #################
    # Input Handler #
    #################
    input_handler = build_input_handler(num_links)
    for r ∈ (1, 13), c = 1
        add_child(arch, input_handler, CartesianIndex(r,c))
    end
    for r ∈ (12, 18), c = 30
        add_child(arch, input_handler, CartesianIndex(r,c))
    end

    ##################
    # Output Handler #
    ##################
    output_handler = build_output_handler(num_links)
    for (r,c) ∈ zip((12,18,1,14), (1, 1, 30, 30))
        add_child(arch, output_handler, CartesianIndex(r,c))
    end

    #######################
    # Global Interconnect #
    #######################
    connect_asap_ldcs_processors(arch, num_links)
    connect_io(arch, num_links)
    connect_memories(arch)
    return arch
end

function build_ldcs_processor_tile(num_links,
                              name = "ldcs_processor_tile",
                              include_memory = false;
                              directions = ("east", "north", "south", "west"),
                             )
    num_fifos = 2
    # Create a new component for the processor tile
    # No need to set the primtiive class or metadata because we won't
    # be needing it.
    comp = Component(name)
    # Add the circuit switched ports
    for dir in directions
        for (suffix,class)  in zip(("_in", "_out"), (Input, Output))
            port_name = join((dir, suffix))
            metadata = make_port_metadata(dir, class, num_links)
            add_port(comp, port_name, class, num_links, metadata = metadata)
        end
    end
    for (suffix,class)  in zip(("_in", "_out"), (Input, Output))
        port_name = join(("ldcs", suffix))
        metadata = make_port_metadata("ldcs", class, num_links)
        add_port(comp, port_name, class, num_links, metadata = metadata)
    end
    # Instantiate the processor primitive
    add_child(comp, build_processor(num_links, include_memory), "processor")

    # Instantiate the directional routing muxes
    routing_mux = build_mux((length(directions)+1),1)
    for dir in directions
        name = "$(dir)_mux"
        add_child(comp, routing_mux, name, num_links)
    end
	# Instantiate the ldcs routing muxes
    routing_mux = build_mux((length(directions)+1),1)
    add_child(comp, routing_mux, "ldcs_mux", num_links)
    # Instantiate the muxes routing data to the fifos
    num_fifo_entries = (length(directions)+1) * num_links + 1
    add_child(comp, build_mux(num_fifo_entries,1), "fifo_mux", num_fifos)
    # Add memory ports - only memory processor tiles will have the necessary
    # "memory_processor" attribute in the core to allow memory application
    # to be mapped to them.
    # Interconnect - Don't attach metadata and let the routing routine fill in
    # defaults to intra-tile routing.
    if include_memory
        metadata = make_port_metadata()
        add_port(comp, "memory_in", "input", metadata = metadata)
        add_port(comp, "memory_out", "output", metadata = metadata)
        add_link(comp, "processor.memory_out", "memory_out")
        add_link(comp, "memory_in", "processor.memory_in")
    end

    # Connect outputs of muxes to the tile outputs
    for dir in directions, i = 0:num_links-1
        mux_port = "$(dir)_mux[$i].out[0]"
        tile_port = "$(dir)_out[$i]"
        add_link(comp, mux_port, tile_port)
    end

    for i = 0:num_links-1
        mux_port = "ldcs_mux[$i].out[0]"
        tile_port = "ldcs_out[$i]"
        add_link(comp, mux_port, tile_port)
    end

    # Circuit switch output links
    for dir in directions, i = 0:num_links-1
        # Make the name for the processor.
        proc_port = "processor.$dir[$i]"
        mux_port = "$(dir)_mux[$i].in[0]"
        add_link(comp, proc_port, mux_port)
    end

    # Connect input fifos.
    for i = 0:num_fifos-1
        fifo_port = "fifo_mux[$i].out[0]"
        proc_port = "processor.fifo[$i]"
        add_link(comp, fifo_port, proc_port)
    end

    directions = ("east", "north", "south", "west", "ldcs")

    index_tracker = Dict((d,i) => 1 for d in directions, i in 0:num_links)
    # Add entries for the fifo
    for i in 0:num_fifos
        index_tracker[("fifo",i)] = 0
    end
    for dir in directions
        for i in 0:num_links-1
            # Create a source port for the tile input.
            source_port = "$(dir)_in[$i]"
            # Begin adding sink ports.
            # Need one for each mux on this layer and num_fifos for each of
            # the input fifos.
            sink_ports = String[]
            # Go through all fifos that directions that are not the current
            # directions.
            for d in Iterators.filter(x -> x != dir, directions)
                # Get the next free port for this mux.
                key = (d,i)
                index = index_tracker[key]
                mux_port = "$(d)_mux[$i].in[$index]"
                # Add this port to the list of sink ports
                push!(sink_ports, mux_port)
                # increment the index tracker
                index_tracker[key] += 1
            end
            # Add the fifo mux entries.
            for j = 0:num_fifos-1
                key = ("fifo", j)
                index = index_tracker[key]
                fifo_port = "fifo_mux[$j].in[$index]"
                push!(sink_ports, fifo_port)
                index_tracker[key] += 1
            end
            # Add the connection to the component.
            add_link(comp, source_port, sink_ports)
        end
    end

    check(comp)
    return comp
end

function connect_asap_ldcs_processors(tl, num_links)

    vals = ["processor", "input_handler", "output_handler"]
    fn = x -> search_metadata!(x, "attributes", vals, oneofin)
    src_rule = fn
    dst_rule = fn

    # Create offset rules.
    offsets = [CartesianIndex(-1,0),
               CartesianIndex(1,0),
               CartesianIndex(0,1),
               CartesianIndex(0,-1),
			   CartesianIndex(-10,0),
			   CartesianIndex(10,0),
			   CartesianIndex(0,-10),
			   CartesianIndex(0,10),]
    #=
    Create two tuples for the source ports and destination ports. In general,
    if the source link is going out of the north port, the destionation will
    be coming in the south port.
    =#
    src_dirs = ("north", "south", "east", "west", "ldcs", "ldcs", "ldcs", "ldcs")
    dst_dirs = ("south", "north", "west", "east", "ldcs", "ldcs", "ldcs", "ldcs")

    src_ports = [["$(src)_out[$i]" for i in 0:num_links-1] for src in src_dirs]
    dst_ports = [["$(dst)_in[$i]" for i in 0:num_links-1] for dst in dst_dirs]

    offset_rules = []
    for (o,s,d) in zip(offsets, src_ports, dst_ports)
        rules = [(o,i,j) for (i,j) in zip(s,d)]
        append!(offset_rules, rules)
    end
    # Build metadata dictionary for capacity and cost
    metadata = Dict(
        "cost"      => 1.0,
        "capacity"  => 1,
        "network"   => ["circuit_switched"]
    )

    connection_rule(tl, offset_rules, src_rule, dst_rule, metadata = metadata)
end


text:
# Sorting Direct Method, implementation following McCollum et al,
# "The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior"
# Comp. Bio. and Chem., 30, pg. 39-49 (2006).

mutable struct SortingDirectJumpAggregation{T,S,F1,F2,RNG,DEPGR} <: AbstractSSAJumpAggregator
    next_jump::Int
    next_jump_time::T
    end_time::T
    cur_rates::Vector{T}
    sum_rate::T
    ma_jumps::S
    rates::F1
    affects!::F2
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    dep_gr::DEPGR
    jump_search_order::Vector{Int}
    jump_search_idx::Int
  end

function SortingDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T,
                                      maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool},
                                      rng::RNG; num_specs, dep_graph=nothing, kwargs...) where {T,S,F1,F2,RNG}

    # a dependency graph is needed and must be provided if there are constant rate jumps
    if dep_graph === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use ConstantRateJumps with the SortingDirect algorithm a dependency graph must be supplied.")
        else
            dg = make_dependency_graph(num_specs, maj)
        end
    else
        dg = dep_graph
    end

    # make sure each jump depends on itself
    add_self_dependencies!(dg)

    # map jump idx to idx in cur_rates
    jtoidx = collect(1:length(crs))
    SortingDirectJumpAggregation{T,S,F1,F2,RNG,typeof(dg)}(nj, njt, et, crs, sr,
                                maj, rs, affs!, sps, rng, dg, jtoidx, zero(Int))
end

########### The following routines should be templates for all SSAs ###########

# condition for jump to occur
@inline function (p::SortingDirectJumpAggregation)(u, t, integrator)
    p.next_jump_time == t
end

# executing jump at the next jump time
function (p::SortingDirectJumpAggregation)(integrator)
    execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t)
    register_next_jump_time!(integrator, p, integrator.t)
    nothing
end

# setting up a new simulation
function (p::SortingDirectJumpAggregation)(dj, u, t, integrator) # initialize
    initialize!(p, integrator, u, integrator.p, t)
    register_next_jump_time!(integrator, p, t)
    nothing
end

############################# Required Functions ##############################

# creating the JumpAggregation structure (function wrapper-based constant jumps)
function aggregate(aggregator::SortingDirect, u, p, t, end_time, constant_jumps,
                   ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)

    build_jump_aggregation(SortingDirectJumpAggregation, u, p, t, end_time, ma_jumps,
                           rates, affects!, save_positions, rng; num_specs=length(u), kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::SortingDirectJumpAggregation, integrator, u, params, t)
    fill_rates_and_sum!(p, u, params, t)
    generate_jumps!(p, integrator, u, params, t)
    nothing
end

# execute one jump, changing the system state
function execute_jumps!(p::SortingDirectJumpAggregation, integrator, u, params, t)
    # execute jump
    u = update_state!(p, integrator, u)

    # update search order
    jso   = p.jump_search_order
    jsidx = p.jump_search_idx
    if jsidx != 1
        @inbounds tmp          = jso[jsidx]
        @inbounds jso[jsidx]   = jso[jsidx-1]
        @inbounds jso[jsidx-1] = tmp
    end

    # update current jump rates
    update_dependent_rates!(p, u, params, t)
    nothing
end

# calculate the next jump / jump time
function generate_jumps!(p::SortingDirectJumpAggregation, integrator, u, params, t)
    @fastmath p.next_jump_time = t + calc_next_jump!(p, u, params, t)
    nothing
end

######################## SSA specific helper routines #########################

# searches down the rate list for the next reaction
@fastmath function calc_next_jump!(p, u, params, t)

    # next jump type
    cur_rates = p.cur_rates
    numjumps  = length(cur_rates)
    jso       = p.jump_search_order
    rn        = p.sum_rate * rand(p.rng)
    @inbounds for idx = 1:numjumps
        rn -= cur_rates[jso[idx]]
        if rn < zero(rn)
            p.jump_search_idx = idx
            break
        end
    end
    @inbounds p.next_jump = jso[p.jump_search_idx]

    # return time to next jump
    randexp(p.rng) / p.sum_rate
end


text:
using SimpleGraphs

"""
    graph_square(G::SimpleGraph)
Create a new graph in which two vertices are adacent if there is a two-step path
between them in `G`.
"""
function graph_square(G::SimpleGraph{T}) where {T}
    A = adjacency(G)
    AA = A^2
    return SimpleGraph(AA)
end


text:
"""
    may(mutate!, args...)

Call `mutate!(args...)` if possible; fallback to the out-of-place version if not.
"""
Base.@propagate_inbounds may(mutate, args...) =
    if possible(mutate, args...)
        mutate(args...)
    else
        pure(mutate)(args...)
    end

_maybb(mutate) = something(_asbb(mutate), mutate)
_asbb(::Any) = nothing

const MaybeMutableContainer = Union{
    AbstractArray,
    AbstractDict,
    AbstractSet,
}

"""
    ismutable(x)

`true` if `x` is mutable as container.
"""
ismutable(x) = ismutable(typeof(x))
ismutable(::Type) = false
ismutable(::Type{<:ImmutableContainer}) = false
ismutable(::Type{<:MaybeMutableContainer}) = true
ismutable(::Type{<:AbstractString}) = false

"""
    ismutablestruct(x)

`true` if `x` is mutable as a struct/object.
"""
ismutablestruct(x) = ismutablestruct(typeof(x))
Base.@pure ismutablestruct(T::DataType) = T.mutable
ismutablestruct(::Type{<:NamedTuple}) = false

# trymutate(::typeof(push!)) = push!!
# trymutate(::typeof(append!)) = append!!

struct Undefined end
ismutable(::Undefined) = false


text:
# test the doc-build

using CalculusWithJulia
CalculusWithJulia.WeaveSupport.weave_all(; build_list=(:html,), force=false)


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

@doc raw"""VolumeMount describes a mounting of a Volume within a container.

    IoK8sApiCoreV1VolumeMount(;
        mountPath=nothing,
        mountPropagation=nothing,
        name=nothing,
        readOnly=nothing,
        subPath=nothing,
        subPathExpr=nothing,
    )

    - mountPath::String : Path within the container at which the volume should be mounted.  Must not contain ':'.
    - mountPropagation::String : mountPropagation determines how mounts are propagated from the host to container and the other way around. When not set, MountPropagationNone is used. This field is beta in 1.10.
    - name::String : This must match the Name of a Volume.
    - readOnly::Bool : Mounted read-only if true, read-write otherwise (false or unspecified). Defaults to false.
    - subPath::String : Path within the volume from which the container's volume should be mounted. Defaults to \"\" (volume's root).
    - subPathExpr::String : Expanded path within the volume from which the container's volume should be mounted. Behaves similarly to SubPath but environment variable references $(VAR_NAME) are expanded using the container's environment. Defaults to \"\" (volume's root). SubPathExpr and SubPath are mutually exclusive.
"""
mutable struct IoK8sApiCoreV1VolumeMount <: SwaggerModel
    mountPath::Any # spec type: Union{ Nothing, String } # spec name: mountPath
    mountPropagation::Any # spec type: Union{ Nothing, String } # spec name: mountPropagation
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    subPath::Any # spec type: Union{ Nothing, String } # spec name: subPath
    subPathExpr::Any # spec type: Union{ Nothing, String } # spec name: subPathExpr

    function IoK8sApiCoreV1VolumeMount(;mountPath=nothing, mountPropagation=nothing, name=nothing, readOnly=nothing, subPath=nothing, subPathExpr=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("mountPath"), mountPath)
        setfield!(o, Symbol("mountPath"), mountPath)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("mountPropagation"), mountPropagation)
        setfield!(o, Symbol("mountPropagation"), mountPropagation)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("subPath"), subPath)
        setfield!(o, Symbol("subPath"), subPath)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("subPathExpr"), subPathExpr)
        setfield!(o, Symbol("subPathExpr"), subPathExpr)
        o
    end
end # type IoK8sApiCoreV1VolumeMount

const _property_map_IoK8sApiCoreV1VolumeMount = Dict{Symbol,Symbol}(Symbol("mountPath")=>Symbol("mountPath"), Symbol("mountPropagation")=>Symbol("mountPropagation"), Symbol("name")=>Symbol("name"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("subPath")=>Symbol("subPath"), Symbol("subPathExpr")=>Symbol("subPathExpr"))
const _property_types_IoK8sApiCoreV1VolumeMount = Dict{Symbol,String}(Symbol("mountPath")=>"String", Symbol("mountPropagation")=>"String", Symbol("name")=>"String", Symbol("readOnly")=>"Bool", Symbol("subPath")=>"String", Symbol("subPathExpr")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1VolumeMount }) = collect(keys(_property_map_IoK8sApiCoreV1VolumeMount))
Swagger.property_type(::Type{ IoK8sApiCoreV1VolumeMount }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1VolumeMount[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1VolumeMount }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1VolumeMount[property_name]

function check_required(o::IoK8sApiCoreV1VolumeMount)
    (getproperty(o, Symbol("mountPath")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1VolumeMount }, name::Symbol, val)
end


text:
abstract type AbstractNode end

"""
    Node

This data type holds the basic Node structure. The type T is used to specify the type of the data
stored in the node.

* If `nchild` is `0` the Node is a leaf node.
* If `root` is `False` the Node is a child of another node.
* `inc_length` specifies the length of the incomming branch.
* `binary` specifies the path from the root to the Node. `1` and `0` represent left and right turns respectively.
"""
mutable struct GeneralNode{S<: AbstractString, R<:Real, A<:AbstractArray{<:Real},
                    C<:AbstractArray{<:Real}, I<:Integer, T<: AbstractString, B<:Bool} <: AbstractNode
    name::S
    data::A
    mother::Union{GeneralNode{S,R,A,C,I,T,B}, Missing}
    children::Vector{GeneralNode{S,R,A,C,I,T,B}}
    scaler::C
    nchild::I
    root::B
    inc_length::R
    binary::T
    num::I
    height::R
    IntExtMap::Vector{I}
    blv::Vector{R}
    initialized::B
    stats::Dict{String, Float64}
end # struct Node

const Node = GeneralNode{String, Float64, Array{Float64, 2}, Array{Float64, 2}, Int64, String, Bool}
const Node_cu = GeneralNode{String, Float64, CuArray{Float64}, CuArray{Float64}, Int64, String, Bool}

"""
    function Node()::Node

This function will initialize an empty node.
"""
function Node()::Node
        Node("no_name", ones(3,3), missing,Node[] ,ones(1,3),0,true,1.0,"0",1,1.0,Int64[],Float64[],false,  Dict{String, Float64}())
end

function Node(name::String; data::Array{A,2}=ones(2,3))::Node where A<:Real
        Node(name, data ,missing, Node[], ones(3,2), 0, true, 1.0, "0", 1, 1.0, Int64[], Float64[], false,  Dict{String, Float64}())
end

#################### Base functionality ####################

Base.:(==)(x::T, y::T) where T<:GeneralNode = x.num == y.num
Base.size(x::T) where T<:GeneralNode = size(post_order(x))
Base.length(x::T) where T<:GeneralNode = x.nchild

function Base.summary(io::IO, d::N) where N <: GeneralNode
    summary(io, d.name)
end

function Base.show(io::IO, d::N) where N <: GeneralNode
    print(io, "Tree with root:\n")
    show(io, d.name)
    if d.initialized
        print(io, "\nLength:\n")
        show(io, "text/plain", tree_length(d))
        print(io, "\nHeight:\n")
        show(io, "text/plain", tree_height(d))
    else
        print(io, "\nLength:\n")
        show(io, "text/plain", 0)
        print(io, "\nHeight:\n")
        show(io, "text/plain", 0)
        print(io, "\nNumber of leave nodes:\n")
        show(io, "text/plain",0)
    end
end

function showall(io::IO, d::N) where N <: GeneralNode
  show(io, d)
  print(io, "\nNode:\n")
  show(io, "text/plain", d.name)
  print(io, "\n#children:\n")
  show(io, d.nchild)
  print(io, "\nbinary:\n")
  show(io, d.binary)
end


text:
using Knet, ArgParse
import Base: start, next, done

function ncelm(args=ARGS)
    info("NCE language model")
    isa(args, AbstractString) && (args=split(args))
    opts = nce_parse_commandline(args)
    println(opts)
    opts["seed"] > 0 && setseed(opts["seed"])
    opts["dropout"] > 0 && error("TODO: implement dropout")
    ftype = opts["float64"]?Float64:Float32
    dict = Dict{Any,Int32}()
    global data = Any[]
    for f in opts["datafiles"]
        push!(data, LMData(f;  dict=dict,
                           batch=opts["batch_size"],
                           seqlen=opts["seq_length"],
                           ftype=ftype,
                           dense=opts["dense"]))
    end

    vocab_size = length(dict)
    psample = fill(ftype(1/vocab_size), vocab_size) # TODO: compute unigram instead

    global net = NCE(nce_rnn; layers = opts["layers"], rnn_size = opts["rnn_size"], vocab = vocab_size, psample=psample, nsample=opts["nsample"])
    lr = opts["lr"]
    setopt!(net, lr=lr, init = Uniform(-opts["init_weight"], opts["init_weight"]))

    perp = zeros(length(data))
    l=zeros(2); m=zeros(2)
    for ep=1:opts["max_max_epoch"]
        ep > opts["max_epoch"] && (lr /= opts["decay"]; setopt!(net, lr=lr))
        train(net, data[1], nothing; psample=psample, nsample=opts["nsample"], gclip=opts["max_grad_norm"], keepstate=true, losscnt=fill!(l,0), maxnorm=fill!(m,0))
        opts["gcheck"]>0 && gradcheck(net,data[1],softloss; gcheck=opts["gcheck"])
        perp[1] = exp(l[1]/l[2])
        for idata = 2:length(data)
            ldev = test(net, data[idata], softloss; keepstate=true)
            perp[idata] = exp(ldev)
        end
        @show (ep, perp..., m..., lr)
    end
    return (perp..., m...)
end

@knet function nce_rnn(word; layers=1, rnn_size=0, rnn_type=lstm, o...)
    wvec = wdot(word; o..., out=rnn_size)
    yrnn = repeat(wvec; o..., frepeat=rnn_type, nrepeat=layers, out=rnn_size)
end

type LMData; data; dict; batchsize; seqlength; ftype; dense; x; y; end

function LMData(fname::AbstractString; batch=20, seqlen=20, dict=Dict{Any,Int32}(), ftype=Float32, dense=false)
    data = Int32[]
    open(fname) do f
        for l in eachline(f)
            for w in split(l)
                push!(data, get!(dict, w, 1+length(dict)))
            end
            push!(data, get!(dict, "<eos>", 1+length(dict))) # end-of-sentence
        end
    end
    info("Read $fname: $(length(data)) words, $(length(dict)) vocab.")
    LMData(data, dict, batch, seqlen, ftype, dense, nothing, nothing)
end

function next(d::LMData,state)                              	# d.data is the whole corpus represented as a sequence of Int32's
    (nword, eos) = state                                        # eos is true if last output was nothing
    if !eos && (nword % (d.batchsize * d.seqlength) == 0)
        return (nothing, (nword, true))                         # output nothing indicating end of sequence
    end
    segsize = div(length(d.data), d.batchsize)                  # we split it into d.batchsize roughly equal sized segments
    offset = div(nword, d.batchsize)                            # this is how many words have been served so far from each segment
    if issparse(d.x)
        for b = 1:d.batchsize
            idata = (b-1)*segsize + offset
            d.x.rowval[b] = d.data[idata+1]
            d.y.rowval[b] = d.data[idata+2]
        end
    else
        fill!(d.x, 0)
        fill!(d.y, 0)
        for b = 1:d.batchsize
            idata = (b-1)*segsize + offset
            d.x[d.data[idata+1], b] = 1
            d.y[d.data[idata+2], b] = 1
        end
    end
    ((d.x, d.y), (nword + d.batchsize, false))	# each call to next will deliver a single word from each of the d.batchsize segments
end

# The state indicates number of words served, and whether the last output was end-of-sequence nothing
function start(d::LMData)
    ndict = length(d.dict)
    if d.x == nothing || size(d.x) != (ndict, d.batchsize)
        if d.dense
            d.x = zeros(d.ftype, ndict, d.batchsize)
            d.y = zeros(d.ftype, ndict, d.batchsize)
        else
            d.x = sponehot(d.ftype, ndict, d.batchsize)
            d.y = sponehot(d.ftype, ndict, d.batchsize)
        end
    end
    return (0, true)
end

function done(d::LMData,state)
    # nword is the number of sequences served
    # stop if there is not enough data for another full sequence
    (nword, eos) = state
    eos && nword + d.batchsize * d.seqlength > length(d.data)
end

function nce_parse_commandline(args)
    s = ArgParseSettings()
    @add_arg_table s begin
        "--batch_size"
        help = "minibatch size"
        arg_type = Int
        default = 128
        "--seq_length"
        help = "length of the input sequence"
        arg_type = Int
        default = 20
        "--layers"
        help = "number of lstm layers"
        arg_type = Int
        default = 1
        "--decay"
        help = "divide learning rate by this every epoch after max_epoch"
        arg_type = Float64
        default = 2.0
        "--rnn_size"
        help = "size of the lstm"
        arg_type = Int
        default = 200
        "--nsample"
        help = "number of noise samples per minibatch"
        arg_type = Int
        default = 100
        "--dropout"
        help = "dropout units with this probability"
        arg_type = Float64
        default =  0.0
        "--init_weight"
        help = "initialize weights uniformly in [-init_weight,init_weight]"
        arg_type = Float64
        default =  0.1
        "--lr"
        help = "learning rate"
        arg_type = Float64
        default = 1.0
        "--vocab_size"
        help = "size of the vocabulary"
        arg_type = Int
        default = 10000
        "--max_epoch"
        help = "number of epochs at initial lr"
        arg_type = Int
        default = 10000 # 4
        "--max_max_epoch"
        help = "number of epochs to train"
        arg_type = Int
        default = 1 # 13
        "--max_grad_norm"
        help = "gradient clip"
        arg_type = Float64
        default = 5.0
        "--gcheck"
        help = "gradient check this many units per layer"
        arg_type = Int
        default = 0
        "--seed"
        help = "random seed, use 0 to turn off"
        arg_type = Int
        default = 42
        "--float64"
        help = "Use Float64 (Float32 default)"
        action = :store_true
        "--dense"
        help = "Use dense matrix ops (sparse ops default)"
        action = :store_true
        "datafiles"
        help = "corpus files: first one will be used for training"
        nargs = '+'
        required = true
    end
    opts = parse_args(args,s)
    return opts
end

!isinteractive() && !isdefined(:load_only) && ncelm(ARGS)

# version 099b284
# [dy_052@hpc3004 examples]$ julia ncelm.jl ptb.valid.txt ptb.test.txt
# WARNING: requiring "Knet" in module "Main" did not define a corresponding module.
# INFO: NCE language model
# Dict{AbstractString,Any}("rnn_size"=>200,"lr"=>1.0,"decay"=>2.0,"dense"=>false,"batch_size"=>20,"float64"=>false,"dropout"=>0.0,"max_grad_norm"=>5.0,"max_epoch"=>10000,"init_weight"=>0.1,"gcheck"=>0,"layers"=>1,"vocab_size"=>10000,"seq_length"=>20,"nsample"=>100,"seed"=>42,"max_max_epoch"=>1,"datafiles"=>Any["ptb.valid.txt","ptb.test.txt"])
# INFO: Read ptb.valid.txt: 73760 words, 6022 vocab.
# INFO: Read ptb.test.txt: 82430 words, 7596 vocab.
#  19.979308 seconds (27.21 M allocations: 1.133 GB, 3.39% gc time)
#   6.901726 seconds (8.21 M allocations: 342.865 MB, 2.84% gc time)
# (ep,perp...,m...,lr) = (1,1.0383780365714608,816.3109741210938,999.4002685546875,38.147953033447266,1.0)