text:
export batched_scal, batched_scal!, batched_gemm!, batched_gemm

import LinearAlgebra: BLAS
import LinearAlgebra.BLAS: libblas, liblapack, @blasfunc, BlasInt

# level 1
"""
    batched_scal!(s, X)

Overwrite `X[:, :, i]` with `a[i] * X[:, :, i]`, where `i` is the batch dimension.
"""
function batched_scal! end

"""
    batched_scal(s, X)

Return `X` scaled by `a` for all the batch.
"""
function batched_scal end

batched_scal(s::AbstractVector{T}, X::AbstractArray{T, 3}) where T = batched_scal!(s, copy(X))

function batched_scal!(A::AbstractVector{T}, B::AbstractArray{T, 3}) where T
    @iterate_batch T A, B (1, 2) begin
        BLAS.scal!(stride(B, 3), A[batch_k], ptrB, 1)
    end
    return B
end

for (fname, elty, lib) in ((:dsyr_,:Float64,libblas),
    (:ssyr_,:Float32,libblas),
    (:zsyr_,:ComplexF64,liblapack),
    (:csyr_,:ComplexF32,liblapack))
    @eval begin
        function batched_syr!(uplo::AbstractChar, α::$elty, x::AbstractArray{$elty, 2}, A::AbstractArray{$elty, 3})
            @assert !has_offset_axes(A, x)
            n = checksquare(A)
            if length(x) != n
            throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))"))
            end

            @iterate_batch $(elty) (x, A) (1, 2) begin
                ccall((@blasfunc($fname), $lib), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                uplo, n, α, ptrx,
                stride(x, 1), ptrA, max(1,stride(A, 2)))
            end
            A
        end
    end
end

for (fname, elty, relty) in ((:zher_,:ComplexF64, :Float64),
                             (:cher_,:ComplexF32, :Float32))
    @eval begin
        function batched_her!(uplo::AbstractChar, α::$relty, x::AbstractMatrix{$elty}, A::AbstractArray{$elty, 3})
            @assert !has_offset_axes(A, x)
            n = checksquare(A)
            if length(x) != n
                throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))"))
            end
            @iterate_batch $(elty) (x, A) (1, 2) begin
                ccall((@blasfunc($fname), libblas), Cvoid,
                    (Ref{UInt8}, Ref{BlasInt}, Ref{$relty}, Ptr{$elty},
                    Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                    uplo, n, α, x,
                    stride(x, 1), A, max(1,stride(A,2)))
            end
            A
        end
    end
end

# TODO: use gemm_batch when mkl is available
"""
    batched_gemm!(transA, transB, alpha, A, B, beta, C)

Batched gemm!.
"""
function batched_gemm! end

"""
    batched_gemm(transA, transB, alpha, A, B)

Batched gemm.
"""
function batched_gemm end

"""
    batched_gemm(alpha, A, B)

"""
batched_gemm(alpha::T, A::AbstractArray{T, 3}, B::AbstractArray{T, 3}) where T = batched_gemm('N', 'N', alpha, A, B)

for (gemm, elty) in
        ((:dgemm_,:Float64),
         (:sgemm_,:Float32),
         (:zgemm_,:ComplexF64),
         (:cgemm_,:ComplexF32))
    @eval begin

        function batched_gemm!(
            transA::AbstractChar, transB::AbstractChar,
            alpha::($elty), A::AbstractArray{($elty), 3}, B::AbstractArray{($elty), 3},
            beta::($elty), C::AbstractArray{($elty), 3})

            @assert !BLAS.has_offset_axes(A, B, C)
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end

            @iterate_batch $(elty) A, B, C (2, 2, 2) begin
            ccall((BLAS.@blasfunc($gemm), BLAS.libblas), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BLAS.BlasInt}, Ref{BLAS.BlasInt},
                     Ref{BLAS.BlasInt}, Ref{$(elty)}, Ptr{$(elty)}, Ref{BLAS.BlasInt},
                     Ptr{$(elty)}, Ref{BLAS.BlasInt}, Ref{$(elty)}, Ptr{$(elty)},
                     Ref{BLAS.BlasInt}),
                     transA, transB, m, n,
                     ka, alpha, ptrA, max(1,stride(A,2)),
                     ptrB, max(1,stride(B,2)), beta, ptrC,
                     max(1,stride(C,2)))
            end
            return C
        end

        function batched_gemm(transA::AbstractChar, transB::AbstractChar, alpha::($elty), A::AbstractArray{$elty, 3}, B::AbstractArray{$elty, 3})
            batched_gemm!(transA, transB, alpha, A, B, zero($elty), similar(B, $elty, (size(A, transA == 'N' ? 1 : 2), size(B, transB == 'N' ? 2 : 1), size(B, 3))))
        end
        function batched_gemm(transA::AbstractChar, transB::AbstractChar, A::AbstractArray{$elty, 3}, B::AbstractArray{$elty, 3})
            batched_gemm(transA, transB, one($elty), A, B)
        end
    end
end


text:
abstract type AbstractLoss end

"Logisitic loss."
struct Logit <: AbstractLoss end

function loss(::Logit, x::T, y::T) where T
    return -log1pexp(x * y)
end

struct LogLoss <: AbstractLoss end

function loss(::LogLoss, x::T, y::T) where T
    return -logloss(x * y)
end

ext: jl


text:
export quick_sort!, hoare_qsort!

function quick_sort!(arr::Vector, h::Int64, t::Int64)
    len = t-h+1
    if len < 2
        return arr
    end

    n, i = h, h
    pivot = arr[t]
    for n = h:t
        if arr[n] <= pivot
            arr[i], arr[n] = arr[n], arr[i]
            i += 1
        end
    end

    if i > t
        i = t
    end

    quick_sort!(arr, h, i-1)
    quick_sort!(arr, i, t)

    return arr
end

function hoare_qsort!(arr::Vector, p::Int64, r::Int64)
    if p >= r
        return arr
    end
    x = arr[p]
    i = p - 1
    j = r + 1
    while true
        while true
            j -= 1
            arr[j] > x || break
        end
        while true
            i += 1
            arr[i] < x || break
        end

        if i < j
            arr[i], arr[j] = arr[j], arr[i]
        else
            break
        end
    end

    hoare_qsort!(arr, p, j);
    hoare_qsort!(arr, j+1, r);

    return arr
end


text:
"A fixed-length vector whose data is stored in a Blob."
struct BlobVector{T} <: AbstractArray{T, 1}
    data::Blob{T}
    length::Int64
end

Base.@propagate_inbounds function get_address(blob::BlobVector{T}, i::Int)::Blob{T} where T
    @boundscheck begin
        (0 < i <= blob.length) || throw(BoundsError(blob, i))
    end
    blob.data + (i-1)*self_size(T)
end

# blob interface

function Base.size(blob::Blob{BlobVector})
    (blob.length[],)
end

function Base.IndexStyle(_::Type{Blob{BlobVector{T}}}) where T
    Base.IndexLinear()
end

Base.@propagate_inbounds function Base.getindex(blob::Blob{BlobVector{T}}, i::Int)::Blob{T} where T
    get_address(blob[], i)
end

# array interface

function Base.size(blob::BlobVector)
    (blob.length,)
end

function Base.IndexStyle(_::Type{BlobVector{T}}) where T
    Base.IndexLinear()
end

Base.@propagate_inbounds function Base.getindex(blob::BlobVector{T}, i::Int)::T where T
    get_address(blob, i)[]
end

Base.@propagate_inbounds function Base.setindex!(blob::BlobVector{T}, v, i::Int)::T where T
    get_address(blob, i)[] = v
end

# copying, with correct handling of overlapping regions
function Base.copy!(
    dest::BlobVector{T}, doff::Int, src::BlobVector{T}, soff::Int, n::Int
) where T
    @boundscheck begin
        if doff < 1 || doff + n - 1 > length(dest)
            throw(BoundsError(dest, doff:doff+n-1))
        elseif soff < 1 || soff + n - 1 > length(src)
            throw(BoundsError(src, soff:soff+n-1))
        end
    end
    # Use memmove for speedy copying. Note: this correctly handles overlapping regions.
    blob_size = Blobs.self_size(T)
    ccall(
        :memmove,
        Ptr{Cvoid},
        (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
        Base.pointer(dest.data) + (doff - 1) * blob_size,
        Base.pointer(src.data) + (soff - 1) * blob_size,
        n * blob_size,
    )
end

# iterate interface

@inline function Base.iterate(blob::BlobVector, i=1)
    (i % UInt) - 1 < length(blob) ? (@inbounds blob[i], i + 1) : nothing
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 1.0.5
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# ### MRI reconstruction demo
# Single-coil under-sampled data reconstructed two ways:
# * edge-preserving regularization (akin to corner-rounding anisotropic TV) using gradient-based methods
# * synthesis sparsity via 1-norm of wavelet detail coefficients using fast proximal methods
#
# Although this is a 2D demo, the k-space sampling pattern used here would be appropriate only for a 3D scan with phase encoding in the two dimensions shown here.  So this demo is akin to a single slice of a 3D scan with 2D undersampling.
#
# This demo also reproduces Figs. 1 & 3 in the survey paper
# "Optimization methods for MR image reconstruction"
# in Jan 2020 IEEE Signal Processing Magazine
# [http://doi.org/10.1109/MSP.2019.2943645]
#
# 2019-03-06 Jeff Fessler, University of Michigan

# load all packages needed for this demo
# https://github.com/JeffFessler/MIRT.jl:
using MIRT: ellipse_im, prompt, embed, Afft, diffl_map, Aodwt
using MIRTjim: jim
using MIRT: ncg, ogm_ls, pogm_restart
using LinearAlgebra
using Plots; default(markerstrokecolor=:auto)
#using LinearMapsAA
using FFTW
using Random: seed!

# ### Create (synthetic) data

# shepp-logan phantom
nx,ny = 192,256
Xtrue = ellipse_im(ny, oversample=2)[Int((ny-nx)/2+1):Int(ny-(ny-nx)/2),:]
#Xtrue = ellipse_im(192, ny=256, oversample=2) # too small
nx,ny = size(Xtrue)
jim(Xtrue, "true image")
#prompt()

# +
#savefig("xtrue.pdf")

# +
# generate noisy, under-sampled k-space data
M,N = nx,ny
seed!(0); sampfrac = 0.3; samp = rand(M,N) .< sampfrac; sig = 1
mod2 = (N) -> mod.((0:N-1) .+ Int(N/2), N) .- Int(N/2)
samp .|= (abs.(mod2(M)) .< Int(M/8)) * (abs.(mod2(N)) .< Int(N/8))' # fully sampled center
@show sum(samp) / (M*N) # under-sampling factor

ytrue = fft(Xtrue)[samp]
y = ytrue + sig * √(2) * randn(ComplexF32, size(ytrue)) # complex noise!
ysnr = 20 * log10(norm(ytrue) / norm(y-ytrue))
@show ysnr;
# -

# display k-space sampling and zero-filled data
logger = (x;min=-6) -> log10.(max.(abs.(x) / maximum(abs.(x)), (10.)^min))
p1 = jim(samp, fft0=true, title="k-space sampling")
#savefig(p1, "samp.pdf")
jim(:abswarn, false) # suppress warnings about showing magnitude
p2 = jim(logger(embed(ytrue,samp)), fft0=true, title="k-space |data|")
plot(p1,p2)

# ### Prepare to reconstruct
# Creating a system matrix (encoding matrix) and an initial image
# The system matrix is a `LinearMapAA` object, akin to a `fatrix` in Matlab MIRT.

# system model from MIRT
F = Afft(samp) # operator!

# initial image based on zero-filled reconstruction
nrmse = (x) -> norm(x - Xtrue) / norm(Xtrue)
X0 = 1/(M*N) * (F' * y)
@show nrmse(X0)
jim(X0, "|X0|: initial image")
#prompt()

# ### Case 1: Edge-preserving regularization
#
# The cost function here is
# $
# 1/2 \|A x - y\|^2 + \beta \; \psi(T x)
# $
# where $\psi$ is the Fair potential function
# and $T$ is first-order finite differences in 2D,
# again implemented as a `LinearMapAA` object.

# Cost function for edge-preserving regularization
A = F
delta = 0.1
reg = 0.01 * M * N
# finite differences sparsifying transform for anisotropic TV:
T = diffl_map((M,N), 1:2 ; T=ComplexF32) # LinearMapAO
pot = (z,del) -> del^2 * (abs(z)/del - log(1 + abs(z)/del)) # Fair potential function
cost = (x) -> 1/2 * norm(A * x - y)^2 + reg * sum(pot.(T*x, delta))
dpot = (z,del) -> z / (1 + abs(z)/del) # potential derivative

# Nonlinear CG recon using edge-preserving regularization
fun = (x,iter) -> (cost(x), nrmse(x), time())
niter = 50
if !@isdefined(Xcg)
    Xcg, out_cg = ncg([A, T],
        [u -> u - y, v -> reg * dpot.(v,delta)], [t -> 1, t -> reg],
        X0, niter=niter, fun=fun)
end
@show nrmse(Xcg)

jim(Xcg, "CG: edge-preserving regularization")
#prompt()

# +
#savefig("xcg.pdf")
# -

# Nesterov FGM recon using edge-preserving regularization
if !@isdefined(Xfgm)
    fun_fgm = (iter,xk,yk,rs) -> (cost(yk), nrmse(yk), rs)
    f_grad = (x) -> A'*(A*x - y) + reg * (T' * dpot.(T * x, delta))
    L_f = N*M + 8*reg
    Xfgm, out_fgm = pogm_restart(X0, cost, f_grad, L_f;
        mom=:fpgm, niter=niter, fun=fun_fgm)
    @show nrmse(Xfgm)
end
jim(Xfgm, "FGM recon with edge-preserving regularization")
#prompt()

# OGM line-search recon using edge-preserving regularization
if !@isdefined(Xogm)
    Xogm, out_ogm = ogm_ls([A, T],
        [u -> u - y, v -> reg * dpot.(v,delta)], [t -> 1, t -> reg],
        X0, niter=niter, ninner=20, fun=fun)
end
@show nrmse(Xogm)

jim(Xogm, "OGM recon with edge-preserving regularization")
#prompt()

p1 = jim(Xtrue, "true")
p2 = jim(X0, "X0: initial")
p3 = jim(Xcg, "CG recon")
p5 = jim(X0 - Xtrue, "X0 error")
p6 = jim(Xcg - Xtrue, "Xcg error")
plot(p1, p2, p3, p1, p5, p6)
#plot(p1, p2, p3, layout=(1,3))
#plot(p5, p6)
#prompt()

# POGM - suboptimal in smooth case
if false
    fun_pogm = (iter,xk,yk,rs) -> (cost(xk), nrmse(xk), rs)
    f_grad = (x) -> A'*(A*x - y) + reg * (T' * dpot.(T * x, delta))
    L_f = N*M + 8*reg
    Xpogm, out_pogm = pogm_restart(X0, cost, f_grad, L_f;
        mom=:pogm, niter=niter, fun=fun_pogm)
    @show nrmse(Xpogm)
end

# plot cost vs iteration
cost_cg = [out_cg[k][1] for k=1:niter+1]
cost_fgm = [out_fgm[k][1] for k=1:niter+1]
cost_ogm = [out_ogm[k][1] for k=1:niter+1]
cost_min = min(minimum(cost_cg), minimum(cost_ogm))
plot(xlabel="iteration k", ylabel="Relative Cost", ytick=[0,8e4])
scatter!(0:niter, cost_fgm .- cost_min, markershape=:square, label="Cost FGM")
scatter!(0:niter, cost_ogm .- cost_min, markershape=:utriangle, label="Cost OGM")
scatter!(0:niter, cost_cg  .- cost_min, label="Cost CG")
#prompt()

# +
#savefig("cost_ogm_cg.pdf")
#savefig("cost_fgm_ogm_cg.pdf")
# -

# plot nrmse vs iteration
nrmse_cg = [out_cg[k][2] for k=1:niter+1]
nrmse_fgm = [out_fgm[k][2] for k=1:niter+1]
nrmse_ogm = [out_ogm[k][2] for k=1:niter+1]
#nrmse_pogm = [out_pogm[k][2] for k=1:niter+1]
plot(xlabel="iteration k", ylabel="NRMSE")
scatter!(0:niter, nrmse_fgm, markershape=:square, label="NRMSE FGM")
scatter!(0:niter, nrmse_ogm, markershape=:utriangle, label="NRMSE OGM")
scatter!(0:niter, nrmse_cg, label="NRMSE CG")
#scatter!(0:niter, nrmse_pogm, label="NRMSE POGM")
#prompt()

#savefig("nrmse_ogm_cg.pdf")

#nrmse(x_pogm), nrmse(X0)

#?pogm_restart

# ### Case 2: Wavelet sparsity in synthesis form
#
# The cost function here is
# $
# 1/2 \|A x - y\|^2 + \beta \; \|W x\|_1
# $
# where
# and $W$ is an orthogonal discrete (Haar) wavelet transform,
# again implemented as a `LinearMapAA` object.
# Because $W$ is unitary,
# we make the change of variables
# $z = W x$
# and solve for $z$
# and then let $x = W' z$
# at the end.

# Orthogonal discrete wavelet transform (LinearMapAO)
W, scales, _ = Aodwt((M,N) ; T = eltype(T))
jim(real(W * Xtrue) .* (scales .> 0), "wavelet detail coefficients")
#prompt()

# Cost function for edge-preserving regularization
Az = F * W'
Fnullz = (z) -> 0 # not needed
f_gradz = (z) -> Az' * (Az * z - y)
f_Lz = M*N # single coil cartesian DFT
regz = 0.03 * M * N # oracle from Xtrue wavelet coefficients!
pot = (z,del) -> del^2 * (abs(z)/del - log(1 + abs(z)/del)) # Fair potential function
costz = (z) -> 1/2 * norm(Az * z - y)^2 + reg * norm(z,1) # 1-norm regularizer
soft = (z,c) -> sign(z) * max(abs(z) - c, 0) # soft thresholding
g_prox = (z,c) -> soft.(z, regz * c) # proximal operator
z0 = W * X0;

# #### Run ISTA=PGM and FISTA=FPGM and POGM, the latter two with adaptive restart
# See [http://doi.org/10.1007/s10957-018-1287-4]

# +
# Run proximal methods
function fun_ista(iter, xk_z, yk, is_restart)
    if mod(iter, 50) == 0
#        @show iter
    end
    xh = W' * xk_z
return (costz(xk_z), nrmse(xh), is_restart) # , psnr(xh)) # time()
end

function fun_fista(iter, xk, yk_z, is_restart)
    xh = W' * yk_z
    return (costz(yk_z), nrmse(xh), is_restart) # , psnr(xh)) # time()
end

z_ista, out_ista = pogm_restart(z0, Fnullz, f_gradz, f_Lz; mom=:pgm, niter=niter,
    restart=:none, restart_cutoff=0., g_prox=g_prox, fun=fun_ista)
Xista = W'*z_ista
@show nrmse(Xista)

z_fista, out_fista = pogm_restart(z0, Fnullz, f_gradz, f_Lz; mom=:fpgm, niter=niter,
    restart=:gr, restart_cutoff=0., g_prox=g_prox, fun=fun_fista)
Xfista = W'*z_fista
@show nrmse(Xfista)

z_pogm, out_pogm = pogm_restart(z0, Fnullz, f_gradz, f_Lz; mom=:pogm, niter=niter,
    restart=:gr, restart_cutoff=0., g_prox=g_prox, fun=fun_fista)
Xpogm = W'*z_pogm
@show nrmse(Xpogm)

jim(Xfista, "FISTA/FPGM")
jim(Xpogm, "POGM with ODWT")
#prompt()

# +
#savefig("xpogm_odwt.pdf")
# -

# plot cost vs iteration
cost_ista = [out_ista[k][1] for k=1:niter+1]
cost_fista = [out_fista[k][1] for k=1:niter+1]
cost_pogm = [out_pogm[k][1] for k=1:niter+1]
cost_min = min(minimum(cost_ista), minimum(cost_pogm))
plot(xlabel="iteration k", ylabel="Relative cost")
scatter!(0:niter, cost_ista  .- cost_min, label="Cost ISTA")
scatter!(0:niter, cost_fista .- cost_min, markershape=:square, label="Cost FISTA")
scatter!(0:niter, cost_pogm  .- cost_min, markershape=:utriangle, label="Cost POGM")
#prompt()

# +
#savefig("cost_pogm_odwt.pdf")
# -

# plot nrmse vs iteration
nrmse_ista = [out_ista[k][2] for k=1:niter+1]
nrmse_fista = [out_fista[k][2] for k=1:niter+1]
nrmse_pogm = [out_pogm[k][2] for k=1:niter+1]
plot(xlabel="iteration k", ylabel="NRMSE", ylim=[0.01,0.08])
scatter!(0:niter, nrmse_ista, label="NRMSE ISTA")
scatter!(0:niter, nrmse_fista, markershape=:square, label="NRMSE FISTA")
scatter!(0:niter, nrmse_pogm, markershape=:utriangle, label="NRMSE POGM")
#prompt()

# +
#savefig("nrmse_pogm_odwt.pdf")
# -

p1 = jim(Xtrue, "true")
p2 = jim(X0, "X0: initial")
p3 = jim(Xpogm, "POGM recon")
p5 = jim(X0 - Xtrue, "X0 error")
p6 = jim(Xpogm - Xtrue, "Xpogm error")
plot(p1, p2, p3, p1, p5, p6)
#plot(p1, p2, p3, layout=(1,3))
#plot(p5, p6)
#prompt()

max_stars_repo_name: chunjiw/mirt-demo

max_issues_repo_name: chunjiw/mirt-demo

max_forks_repo_name: chunjiw/mirt-demo


text:
export Operator
module Operator
using KernelAbstractions.Extras: @unroll

using CUDA
using GemmKernels
using GemmKernels.Tiling

# -------------------------------------
# Default definition for padded layouts
# -------------------------------------

for f in (:fragtype_a, :fragtype_b, :fragtype_accum, :load_a, :load_b, :load_c, :store_d)
    @eval @inline $f(op, ::Type{Layout.Padded{L, P}}, args...) where {L, P} = $f(op, L, args...)
end

# ----
#  SIMT Op
# ----
struct SIMTOp end

@inline shape(::Type{SIMTOp}) where {M, N, K} = (M = 8, N = 4, K = 1)

# convert_index_func: function used to transpose the index in case of a row-major layout
# 2 types of optimizations:
# 1. coalescing of memory ops
# 2. vectorization of memory ops (like vload in layout.jl?)
for (layout_type, convert_index_func) in [
                                        (Layout.AlignedColMajor, identity),
                                        (Layout.AlignedRowMajor, x -> reverse(Tuple(x)))
                                       ]

    @eval begin
        # fragtype: type of per thread input/result of mma operations
        @inline fragtype_a(::Type{SIMTOp}, ::Type{$layout_type{T}}) where {T} = T

        @inline fragtype_b(::Type{SIMTOp}, ::Type{$layout_type{T}}) where {T} = T
        @inline fragtype_accum(::Type{SIMTOp}, ::Type{$layout_type{T}}) where {T} = T

        # define load / stores based on layout types
        @inline function load_a(::Type{SIMTOp}, ::Type{$layout_type{T}}, workspace, tile::Tile) where {T}
            laneId = (threadIdx().x - 1) % 32 + 1
            row = ((laneId - 1) % 8)

            y, x = $convert_index_func((tile.base.M  + tile.offset.M + 1 + row, tile.base.K + tile.offset.K + 1))
            @inbounds return workspace[y, x]
        end

        @inline function load_b(::Type{SIMTOp}, ::Type{$layout_type{T}}, workspace, tile::Tile) where {T}
            laneId = (threadIdx().x - 1) % 32 + 1
            col = ((laneId - 1) ÷ 8)

            y, x = $convert_index_func((tile.base.K + tile.offset.K + 1, tile.base.N + tile.offset.N + 1 + col))
            @inbounds return workspace[y, x]
        end

        @inline function load_c(op::Type{SIMTOp}, layout::Type{$layout_type{T}}, workspace, tile::Tile) where {T}
            laneId = (threadIdx().x - 1) % 32 + 1
            row = ((laneId - 1) % 8)
            col = ((laneId - 1) ÷ 8)

            y, x = $convert_index_func((tile.base.M + tile.offset.M + 1 + row, tile.base.N + tile.offset.N + 1 + col))
            @inbounds return workspace[y, x]
        end

        @inline function store_d(::Type{SIMTOp}, ::Type{$layout_type{T}}, workspace, frag, tile::Tile) where {T}
            laneId = (threadIdx().x - 1) % 32 + 1

            row = ((laneId - 1) % 8)
            col = ((laneId - 1) ÷ 8)

            y, x = $convert_index_func((tile.base.M + tile.offset.M + 1 + row, tile.base.N + tile.offset.N + 1 + col))
            @inbounds workspace[y, x] = frag
        end
    end
end

@inline function mma(op::Type{SIMTOp}, a_frag, b_frag, c_frag)
    @inbounds return a_frag * b_frag + c_frag
end

# ----

struct WMMAOp{M, N, K, T} end

@inline shape(::Type{WMMAOp{M, N, K, T}}) where {M, N, K, T} = (M = M, N = N, K = K)

# convert_index_func: function used to transpose the index in case of a row-major layout
for (layout_type, wmma_layout_type, convert_index_func) in [
                                        (Layout.AlignedColMajor, WMMA.ColMajor, identity),
                                        (Layout.AlignedRowMajor, WMMA.RowMajor, x -> reverse(Tuple(x)))
                                       ]
    @eval begin
        @inline fragtype_a(::Type{WMMAOp{16, 16, 16, T}}, ::Type{$layout_type{Float16}}) where {T} = WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixA}
        @inline fragtype_b(::Type{WMMAOp{16, 16, 16, T}}, ::Type{$layout_type{Float16}}) where {T} = WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixB}
        @inline fragtype_accum(::Type{WMMAOp{16, 16, 16, T}}, ::Type{$layout_type{T}}) where {T} = WMMA.Fragment{16, 16, 16, 8, T, WMMA.Unspecified, WMMA.Accumulator}

        @inline function load_a(::Type{WMMAOp{M, N, K, T}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K, T}
            conf = WMMA.Config{M, N, K, T}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(Float16)
            return WMMA.load_a(ptr, size(workspace, 1), $wmma_layout_type, conf)
        end

        @inline function load_b(::Type{WMMAOp{M, N, K, T}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K, T}
            conf = WMMA.Config{M, N, K, T}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(Float16)
            return WMMA.load_b(ptr, size(workspace, 1), $wmma_layout_type, conf)
        end

        @inline function load_c(::Type{WMMAOp{M, N, K, T}}, ::Type{$layout_type{T}}, workspace, tile::Tile) where {M, N, K, T}
            conf = WMMA.Config{M, N, K, T}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(T)
            return WMMA.load_c(ptr, size(workspace, 1), $wmma_layout_type, conf)
        end

        @inline function store_d(::Type{WMMAOp{M, N, K, T}}, ::Type{$layout_type{T}}, workspace, frag, tile::Tile) where {M, N, K, T}
            conf = WMMA.Config{M, N, K, T}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(T)
            WMMA.store_d(ptr, frag, size(workspace, 1), $wmma_layout_type, conf)
        end
    end
end

function mma(::Type{WMMAOp{M, N, K, T}}, a_frag, b_frag, c_frag) where {M, N, K, T}
    conf = WMMA.Config{M, N, K, T}
    return WMMA.mma(a_frag, b_frag, c_frag, conf)
end

# -----------
# WMMAComplex
# -----------

struct WMMAComplexOp{M, N, K} end

@inline shape(::Type{WMMAComplexOp{M, N, K}}) where {M, N, K} = (M = M, N = N, K = K)

# convert_index_func: function used to transpose the index in case of a row-major layout
for (layout_type, wmma_layout_type, convert_index_func) in [
                                        (Layout.SplitColMajor, WMMA.ColMajor, identity),
                                        (Layout.SplitRowMajor, WMMA.RowMajor, x -> reverse(Tuple(x))),
                                       ]
    @eval begin
        @inline fragtype_a(::Type{WMMAComplexOp{16, 16, 16}}, ::Type{$layout_type{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixA}}
        @inline fragtype_b(::Type{WMMAComplexOp{16, 16, 16}}, ::Type{$layout_type{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixB}}
        @inline fragtype_accum(::Type{WMMAComplexOp{16, 16, 16}}, ::Type{$layout_type{Float32}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 8, Float32, WMMA.Unspecified, WMMA.Accumulator}}

        @inline function load_a(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{16, 16, 16, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            return (WMMA.load_a(pointer(workspace, ind), size(workspace)[1], $wmma_layout_type, conf),
                    WMMA.load_a(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], $wmma_layout_type, conf))
        end

        @inline function load_b(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{16, 16, 16, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            return (WMMA.load_b(pointer(workspace, ind), size(workspace)[1], $wmma_layout_type, conf),
                    WMMA.load_b(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], $wmma_layout_type, conf))
        end

        @inline function load_c(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float32}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            return (WMMA.load_c(pointer(workspace, ind), size(workspace)[1], $wmma_layout_type, conf),
                    WMMA.load_c(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], $wmma_layout_type, conf))
        end

        @inline function store_d(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float32}}, workspace, frag, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            WMMA.store_d(pointer(workspace, ind), frag[1], size(workspace)[1], $wmma_layout_type, conf)
            WMMA.store_d(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), frag[2], size(workspace)[1], $wmma_layout_type, conf)
        end
    end
end

using LLVM

@inline function mma(::Type{WMMAComplexOp{M, N, K}}, a_frag, b_frag, c_frag) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}

    c_re = c_frag[1]
    c_im = c_frag[2]

    c_re = WMMA.mma(a_frag[1],  b_frag[1], c_re, conf)
    c_re = WMMA.mma(.-(a_frag[2]), b_frag[2], c_re, conf)

    c_im = WMMA.mma(a_frag[1], b_frag[2], c_im, conf)
    c_im = WMMA.mma(a_frag[2], b_frag[1], c_im, conf)

    return (c_re, c_im)
end

# --------
# WMMADual
# --------

struct WMMADualOp{M, N, K} end

@inline shape(::Type{WMMADualOp{M, N, K}}) where {M, N, K} = (M = M, N = N, K = K)

@inline fragtype_a(::Type{WMMADualOp{16, 16, 16}}, ::Type{Layout.SplitColMajor{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, WMMA.ColMajor, WMMA.MatrixA}}
@inline fragtype_b(::Type{WMMADualOp{16, 16, 16}}, ::Type{Layout.SplitColMajor{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, WMMA.ColMajor, WMMA.MatrixB}}
@inline fragtype_accum(::Type{WMMADualOp{16, 16, 16}}, ::Type{Layout.SplitColMajor{Float32}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 8, Float32, WMMA.Unspecified, WMMA.Accumulator}}

@inline function load_a(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float16}}, workspace, tile::Tile) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    return (WMMA.load_a(pointer(workspace, ind), size(workspace)[1], WMMA.ColMajor, conf),
            WMMA.load_a(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], WMMA.ColMajor, conf))
end

@inline function load_b(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float16}}, workspace, tile::Tile) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    return (WMMA.load_b(pointer(workspace, ind), size(workspace)[1], WMMA.ColMajor, conf),
            WMMA.load_b(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], WMMA.ColMajor, conf))
end

@inline function load_c(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float32}}, workspace, tile::Tile) where {M, N, K}
    conf = WMMA.Config{M, N, K, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    return (WMMA.load_c(pointer(workspace, ind), size(workspace)[1], WMMA.ColMajor, conf),
            WMMA.load_c(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], WMMA.ColMajor, conf))
end

@inline function store_d(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float32}}, workspace, frag, tile::Tile) where {M, N, K}
    conf = WMMA.Config{M, N, K, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    WMMA.store_d(pointer(workspace, ind), frag[1], size(workspace)[1], WMMA.ColMajor, conf)
    WMMA.store_d(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), frag[2], size(workspace)[1], WMMA.ColMajor, conf)
end

@inline function mma(::Type{WMMADualOp{M, N, K}}, a_frag, b_frag, c_frag) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}

    c_re = c_frag[1]
    c_du = c_frag[2]

    c_re = WMMA.mma(a_frag[1],  b_frag[1], c_re, conf)

    c_du = WMMA.mma(a_frag[1], b_frag[2], c_du, conf)
    c_du = WMMA.mma(a_frag[2], b_frag[1], c_du, conf)

    return (c_re, c_du)
end

end


text:
module PoincareInvariant2ndTest

    using PoincareInvariants
    using GeometricIntegrators
    using GeometricIntegrators.Utils
    using SymPy

    const nx = 100
    const ny = 100
    const Δt = 10.
    const B₀ = 1.
    const r₀ = 0.5
    const z₀ = 0.0
    const z₁ = 0.1
    const u₀ = 5E-1
    const u₁ = 5E-2
    const kz = 2
    const ku = 4

    function B(t, q)
        B₀ * (1 + q[1]^2 + q[2]^2)
    end

    function gcϑ(t, q, ϑ)
        ϑ[1] =-q[2] * (1 + B(t,q)) / 4
        ϑ[2] =+q[1] * (1 + B(t,q)) / 4
        ϑ[3] = q[4]
        ϑ[4] = 0
        nothing
    end

    function gcω(t, q, Ω)
        Ω[1,1] = 0
        Ω[1,2] =-B(t,q)
        Ω[1,3] = 0
        Ω[1,4] = 0

        Ω[2,1] =+B(t,q)
        Ω[2,2] = 0
        Ω[2,3] = 0
        Ω[2,4] = 0

        Ω[3,1] = 0
        Ω[3,2] = 0
        Ω[3,3] = 0
        Ω[3,4] =+q[4]

        Ω[4,1] = 0
        Ω[4,2] = 0
        Ω[4,3] =-q[4]
        Ω[4,4] = 0

        nothing
    end

    function gc_surface_q(s,t)
        x  = r₀*(s-0.5)
        y  = r₀*(t-0.5)
        z  = z₀ + z₁ * cos(2π*kz*s) * cos(2π*kz*t)
        u  = u₀ + u₁ * sin(2π*ku*s) * sin(2π*ku*t)

        [x, y, z, u]
    end

    function gc_surface_p(s,t)
        q = gc_surface_q(s,t)
        p = zero(q)
        gcϑ(zero(eltype(q)), q, p)
        p
    end

    function gc_dummy(t, a, b, c)
        nothing
    end

    function gc_dummy_pode(q₀, p₀)
        PODE(gc_dummy, gc_dummy, q₀, p₀)
    end

    function gc_dummy_iode(q₀)
        p₀ = zero(q₀)

        if ndims(q₀) == 1
            gcϑ(zero(eltype(q₀)), q₀, p₀)
        else
            tq = zeros(eltype(q₀), size(q₀,1))
            tp = zeros(eltype(p₀), size(p₀,1))

            for i in 1:size(q₀,2)
                simd_copy_xy_first!(tq, q₀, i)
                gcϑ(zero(eltype(q₀)), tq, tp)
                simd_copy_yx_first!(tp, p₀, i)
            end
        end

        IODE(gc_dummy, gc_dummy, gc_dummy, gc_dummy, q₀, p₀)
    end

    function compute_canonical_invariant_numerical()
        pinv = PoincareInvariant2ndCanonical(gc_dummy_pode, gc_surface_q, gc_surface_p, Δt, 4, nx, ny, 0)
        sol  = Solution(pinv.equ, Δt, 0)

        I = evaluate_poincare_invariant(pinv, sol)

        return I[1]
    end

    function compute_invariant_numerical()
        pinv = PoincareInvariant2nd(gc_dummy_iode, gc_surface_q, gcω, Δt, 4, nx, ny, 0)
        sol  = Solution(pinv.equ, Δt, 0)

        I, J, K, L = evaluate_poincare_invariant(pinv, sol)

        return I[1], J[1]
    end

    function compute_invariant_analytical()
        s, t = Sym("s, t")
        x, y, z, u = gc_surface_q(s,t)

        S = simplify(B₀ * (1 + x^2 + y^2) * ( diff(x, s) * diff(y, t) + diff(x, t) * diff(y, s) )
                   + u * ( diff(z, s) * diff(u, t) + diff(z, t) * diff(u, s) ) )
        I = N(SymPy.integrate(S, (s, 0, 1), (t, 0, 1)), 20)
    end

    export compute_canonical_invariant_numerical, compute_invariant_numerical, compute_invariant_analytical

end

using .PoincareInvariant2ndTest

### compute analytical invariant ###

Ia = compute_invariant_analytical()

### compute and check canonical invariant ###

Jn = compute_canonical_invariant_numerical()

@assert Jn ≈ Ia atol=2eps()

### compute and check noncanonical invariant ###

In, Jn = compute_invariant_numerical()

@assert In ≈ Ia atol=2eps()
@assert Jn ≈ Ia atol=2eps()


text:
@testset "variablelength" begin

    testvalues = [
        ([UInt8(0b10000001), UInt8(0b00000001)], 0b10000001),
        ([UInt8(0b00001111)], 0b00001111),
        ([UInt8(0b10011001), UInt8(0b00100101)], 0b0000110010100101),
        ([UInt8(0b10000001), UInt8(0b01111111)], 0b11111111)
    ]

    @testset "it should correctly read a variable length number" begin
        for (input, output) in testvalues
            @test MIDI.readvariablelength(IOBuffer(input)) == output
        end
    end

    @testset "it should correctly write a variable length number" begin
        for (output, input) in testvalues
            buf = IOBuffer()
            MIDI.writevariablelength(buf, Int(input))
            @test take!(buf) == output
        end
    end
end


text:
function ∂  end
function ∂ₛ end

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

function (sim::Simulation)(::typeof(∂),dir::Direction,gridtype::GridType)
    D(dir,sim.grid,sim.e⁻ⁱᵏᴸ,gridtype)
end

function (sim::Simulation)(::typeof(∂ₛ),dir::Direction,gridtype::GridType)
       S =  get_S_comp(sim,  dir)
       spdiagm((1 ./ S)[:])  *  sim(∂,dir,gridtype)
end

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

  function D(dir::Direction{NDir},g::Grid{NDim},e⁻ⁱᵏᴸ::NTuple,gridtype::GridType) where {NDim,NDir}
      if NDir > NDim return spzeros(ncells(g),ncells(g)) end
      dw = spacing(g,dir);
      nw = size(g,dir);
      ∂w = D(nw,e⁻ⁱᵏᴸ[NDir],gridtype) / dw;
      return kron(∂w, dir, g)
    end

 ################################################################################
 function  D(N::Int, e⁻ⁱᵏᴸ::Number,gridtype::GridType{:Primal})
    T = promote_type(Float, eltype(e⁻ⁱᵏᴸ))
    Df =  spdiagm(0=>  -ones(T,N),1=> ones(T,N-1))
    Df[N,1] =  1.0 * e⁻ⁱᵏᴸ
    return  Df
 end

 function  D(N::Int, e⁻ⁱᵏᴸ::Number,gridtype::GridType{:Dual})
    T = promote_type(Float, eltype(e⁻ⁱᵏᴸ))
    Db =  spdiagm(0=> ones(T,N),  -1=> -ones(T,N-1))
    Db[1,N] = -1.0 * (e⁻ⁱᵏᴸ)^-1
    return  Db
 end
 #############################################################################


text:
using Test
using Ju
using Random:seed!
using StatsBase:mean

@testset "Ju" begin

include("base.jl")
include("approximators.jl")
include("spaces.jl")
include("buffers.jl")
include("action_selectors.jl")
include("policies.jl")
include("environments.jl")

end


text:
include("main.jl")

BUSES, LINES, GENERATORS, rPTDF = load_case_data()
ERROR_VARIANCES_TRUE, ERROR_DIST_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

# objective, bus_results, line_results = Solve_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_SAMPLE)
# ddobjective, bus_results, line_results = Solve_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_WC)

buses_lc = change_load_same_pf(1)

objective, bus_results, line_results = Solve_DD_DR_CCLinDist(buses_lc, LINES, GENERATORS, ERROR_VARIANCES_SAMPLE)
ddobjective, ddbus_results, ddline_results = Solve_DD_DR_CCLinDist(buses_lc, LINES, GENERATORS, ERROR_VARIANCES_WC)

info("Results of standard cc. Objective: $objective")
println(bus_results)
println(line_results)

info("Results of dd cc. Objective: $ddobjective")
println(ddbus_results)
println(ddline_results)

sample_size = 1000 # max = 10000
# Create global set of true observations to make tests comparable
OBSERVATIONS = rand(ERROR_DIST_TRUE, sample_size)
test_results(bus_results, line_results, objective; samples = sample_size, message="Standard CC")
test_results(ddbus_results, ddline_results, ddobjective; samples = sample_size, message="DD CC")


text:
struct Certificate
    id::String
    name::String
    notafter::DateTime
    sha1fingerprint::String
    created_at::DateTime
    dnsnames::Array{String, 1}
    state::String
    certificatetype::String

    function Certificate(data::Dict{String})
        data["not_after"] = DateTime(data["not_after"][1:end-1])
        data["created_at"] = DateTime(data["created_at"][1:end-1])

        new(
            data["id"],
            data["name"],
            data["not_after"],
            data["sha1_fingerprint"],
            data["created_at"],
            data["dns_names"],
            data["state"],
            data["type"]
        )
    end
end

function show(io::IO, c::Certificate)
    print(io, "Certificate ($(c.name))")
end

# List all certificates
function getallcertificates!(client::AbstractClient)
    uri = joinpath(ENDPOINT, "certificates?per_page=$MAXOBJECTS")
    getalldata!(client, uri, Certificate)
end

# Retrieve an existing certificate
function getcertificate!(client::AbstractClient, certificateid::String)
    uri = joinpath(ENDPOINT, "certificates", "$(certificateid)")
    Certificate(getdata!(client, uri))
end

function getcertificate!(client::AbstractClient, certificate::Certificate)
    getcertificate!(client, certificate.id)
end

function createcertificate!(client::AbstractClient; name::String, private_key::String,
                            leaf_certificate::String, kwargs...)
    postbody = Dict{String, Any}("name" => name, "private_key" => private_key,
                                 "leaf_certificate" => leaf_certificate)
    merge!(postbody, Dict{String, Any}([String(k[1]) => k[2] for k in kwargs]))

    uri = joinpath(ENDPOINT, "certificates")
    Certificate(postdata!(client, uri, postbody))
end

function deletecertificate!(client::AbstractClient, certificateid::String)
    uri = joinpath(ENDPOINT, "certificates", certificateid)
    deletedata!(client, uri)
end

function deletecertificate!(client::AbstractClient, certificate::Certificate)
    deletecertificate!(client, certificate.id)
end


text:
function GAMMAH(N::Integer, over::Real)
    #INTEGER N,I,J
    #DOUBLE PRECISION OVER
    I = N
    J = 2 * I - 1
    retval = 1.0
    while (J>=1) && retval < over
         retval = retval * J/2.0
         I = I-1
         J = 2*I-1
    end
    if J > 1
        retval = 0.0
    end
    return retval
end

function PSI(M::Integer, K::Integer)
    FACTR1=0.0
    FACTR2=0.0
    N=2*K+M
    J=1
    if K>=1
        while J <= K
            FACTR1 = FACTR1+1.0/J
            J += 1
        end
    end
    I = 1
    while I <= N
        FACTR2 = FACTR2 + 1.0/(2.0*I-1.0)
        I=I+1
    end
    if N != 0
        FACTR2 = 0.0
    end

    return FACTR1-2.0*FACTR2

end

function FRAC(Z::Real, M::Integer, NMAX::Integer, ε::Real, TINYSQ::Real)
    N = NMAX
    MM = 0
    DZ2 = Z + Z
    DN0 = N + M
    DN1 = N + 1.0
    DN2 = DN0 + 0.5
    DN3 = DN0 - 0.5
    DN4 = 1.0 - M - M
    B = 2.0*DN1*Z/DN2
    A = 1.0
    FC = TINYSQ
    C0 = FC
    D0 = 0.0
    while true
        D0 = B + A*D0
        if abs(D0)  < TINYSQ
            D0 = TINYSQ
        end
        C0 = B + A/C0
        if abs(C0) < TINYSQ
            C0 = TINYSQ
        end
        D0 = 1.0/D0
        DELTA = C0*D0
        FC = FC*DELTA
        MM = MM + 1
        A = -(1.0 + DN4/(DN3 + MM))
        B = DZ2*(DN1 + MM)/(DN2 + MM)
        if MM < 1000000
            if abs(DELTA-1.0) > ε
                continue
            else
                return FC
            end
        elseif MM == 1000000
            println("CF CONVERGENCE FAILS")
            return FC
        end
    end
    return FC
end

function FACTOR(M::Integer, K::Integer)
    retval = 1.0
    if K >= 1
        X1 = M + 0.5
        N = 2*K - 1
        I = K
        J = N
        while J >= 0
            retval = retval*(J+X1)/I/I
            J = J-1
            I = I -1
            if I == 0
                I = 1
            end
        end
    end
    return retval
end

function EXPAN(Z::Real, MODE::Integer, IPRE::Integer, OVER::Real, QARGU::Real, M::Integer)
    PRECI = [1.e-13 1.e-9]
    πSQ=sqrt(π)
    DB = 2.0*log(2.0)
    FL = M/2.0
    CC = abs(trunc(FL)-FL)
    if (CC < 0.4)
        AR = 1.0
    else
        AR = -1.0
    end
    DZ = 1.0
    for I = 1:M
        DZ = DZ/QARGU
    end
    if MODE == 0
        GAMMA = GAMMAH(M, OVER)*AR*πSQ
    else
        GAMMA = AR
    end
    DFAC = 2.0/π*DZ*GAMMA/πSQ
    DF1 = log(2.0*Z)
    A0 = 1.0/sqrt(2.0*Z)
    Z2I = 1.0/(Z*Z)
    DELTA = 1.0
    SUM = 0.0
    K = 0
    DA2 = FACTOR(M, 0)
    DA1 = DB + PSI(M, 0)
    while true
        DELTA = (DF1+DA1)*DA2*A0
        SUM = SUM+DELTA
        DCC = 0.5+M+2.0*K
        DCCP = DCC + 1.0
        DKK = K + 1.0
        DA2 = DA2*DCCP*DCC/(DKK*DKK)
        DA1 = DA1+1.0/DKK-1.0/DCCP-1.0/DCC
        K = K + 1
        A0 = A0*0.25*Z2I
        if abs(DELTA/SUM) > PRECI[IPRE]
            continue
        else
            break
        end
    end
    return SUM*DFAC
end

function FRACPS(QZ::Real, M::Integer, N::Integer, ε::Real, TINYSQ::Real)
    MM = 0
    DN2 = N*N
    DZ2 = QZ+QZ
    DM = M-0.5
    B = DZ2*M
    A = DN2-DM*DM
    FC = TINYSQ
    C0 = FC
    D0 = 0.0
    while true
        D0 = B + A*D0
        if abs(D0) < TINYSQ
            D0 = TINYSQ
        end
        C0 = B + A/C0
        if abs(C0) < TINYSQ
            C0 = TINYSQ
        end
        D0 = 1.0/D0
        DELTA = C0*D0
        FC = FC*DELTA
        MM = MM + 1
        A = DN2-(MM+DM)*(MM+DM)
        B = DZ2*(MM+M)
        if MM < 1000000
            if  abs(DELTA-1.0) > ε
                continue
            else
                return FC
            end
        elseif MM == 1000000
            println("CF CONVERGENCE FAILS")
            return FC
        end
    end
    return FC
end

function FACTCO(N::Integer, PL::Real, M::Integer)
    retval = 1.0/PL
    X1 =  M + 0.50
    X2 = -M + 0.50
    J = N
    while J >= 0
        retval = retval*(J+X1)/(J+X2)
        J -= 1
    end
    return retval
end

function ELLIP2(AK)
    Q=(1.0-AK)*(1.0+AK)
    return DRF(Q)-(AK)^2*DRD(Q)/3.0
end

function ELLIP1(AK)
    return DRF((1.0-AK)*(1.0+AK))
end

function DRF(Y::Real)
    ERRTOL = 1.e-8
    C1 = 1.0/24.0
    C2 = 3.0/44.0
    C3 = 1.0/14.0

    DRF = 0.0
    XN = 0.
    YN = Y
    ZN = 1.0

    MU = 0.0
    XNDEV = 0.0
    YNDEV = 0.0
    ZNDEV = 0.0

    while true
        MU = (XN + YN + ZN)/3.0
        XNDEV = 2.0 - (MU+XN)/MU
        YNDEV = 2.0 - (MU+YN)/MU
        ZNDEV = 2.0 - (MU+ZN)/MU
        EPSLON = max(abs(XNDEV), abs(YNDEV), abs(ZNDEV))
        if EPSLON < ERRTOL
            break
        end
        XNROOT = sqrt(XN)
        YNROOT = sqrt(YN)
        ZNROOT = sqrt(ZN)
        LAMDA = XNROOT*(YNROOT + ZNROOT) + YNROOT*ZNROOT
        XN = (XN + LAMDA)*0.250
        YN = (YN + LAMDA)*0.250
        ZN = (ZN + LAMDA)*0.250
    end
    E2 = XNDEV*YNDEV - ZNDEV*ZNDEV
    E3 = XNDEV*YNDEV*ZNDEV
    S  = 1.0 + (C1*E2-0.10-C2*E3)*E2 + C3*E3
    return S/sqrt(MU)
end

function DRD(Y::Real)
    ERRTOL = 1.e-8
    C1 = 3.0/14.0
    C2 = 1.0/6.0
    C3 = 9.0/22.0
    C4 = 3.0/26.0

    DRD = 0.0

    XN = 0.0
    YN = Y
    ZN = 1.0
    SIGMA = 0.0
    POWER4 = 1.0

    MU = 0.0
    XNDEV = 0.0
    YNDEV = 0.0
    ZNDEV = 0.0

    while true
        MU = (XN + YN + 3.0*ZN)*0.20
        XNDEV = (MU-XN)/MU
        YNDEV = (MU-YN)/MU
        ZNDEV = (MU-ZN)/MU
        EPSLON = max(abs(XNDEV), abs(YNDEV), abs(ZNDEV))
        if EPSLON < ERRTOL
            break
        end
        XNROOT = sqrt(XN)
        YNROOT = sqrt(YN)
        ZNROOT = sqrt(ZN)
        LAMDA = XNROOT*(YNROOT + ZNROOT) + YNROOT*ZNROOT
        SIGMA = SIGMA + POWER4/(ZNROOT*(ZN+LAMDA))
        POWER4 = POWER4*0.250
        XN = (XN+LAMDA)*0.250
        YN = (YN+LAMDA)*0.250
        ZN = (ZN+LAMDA)*0.250
    end

    EA = XNDEV*YNDEV
    EB = ZNDEV*ZNDEV
    EC = EA - EB
    ED = EA - 6.0*EB
    EF = ED + EC + EC
    S1 = ED*(-C1+0.250*C3*ED-1.50*C4*ZNDEV*EF)
    S2 = ZNDEV*(C2*EF+ZNDEV*(-C3*EC+ZNDEV*C4*EA))
    return 3.0*SIGMA + POWER4*(1.0 + S1 + S2)/(MU*sqrt(MU))
 end


text:
##
# ProductFun represents f(x,y) by Fun(.coefficients[k](x),.space[2])(y)
# where all coefficients are in the same space
##

export ProductFun

immutable ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,T}<:BivariateFun{T}
    coefficients::Vector{Fun{S,T}}     # coefficients are in x
    space::SS
end

ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,T<:Number,P}(cfs::Vector{Fun{S,T}},sp::AbstractProductSpace{Tuple{S,V},P,2}) =
    ProductFun{S,V,typeof(sp),T}(cfs,sp)
ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,
         W<:UnivariateSpace,T<:Number,P}(cfs::Vector{Fun{S,T}},sp::AbstractProductSpace{Tuple{W,V},P,2}) =
   ProductFun{W,V,typeof(sp),T}(Fun{W,T}[Fun(cfs[k],columnspace(sp,k)) for k=1:length(cfs)],sp)

Base.size(f::ProductFun,k::Integer) =
    k==1?mapreduce(ncoefficients,max,f.coefficients):length(f.coefficients)
Base.size(f::ProductFun) = (size(f,1),size(f,2))

## Construction in an AbstractProductSpace via a Matrix of coefficients

function ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,T<:Number,P}(cfs::Matrix{T},sp::AbstractProductSpace{Tuple{S,V},P,2};
                                                                       tol::Real=100eps(T),chopping::Bool=false)
    if chopping
        ncfs,kend=norm(cfs,Inf),size(cfs,2)
        if kend > 1 while isempty(chop(cfs[:,kend],ncfs*tol)) kend-=1 end end
        ret=Fun{S,T}[Fun(chop(cfs[:,k],ncfs*tol),columnspace(sp,k)) for k=1:max(kend,1)]
        ProductFun{S,V,typeof(sp),T}(ret,sp)
    else
        ret=Fun{S,T}[Fun(cfs[:,k],columnspace(sp,k)) for k=1:size(cfs,2)]
        ProductFun{S,V,typeof(sp),T}(ret,sp)
    end
end

## Construction in a ProductSpace via a Vector of Funs

function ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,T<:Number}(M::Vector{Fun{S,T}},dy::V)
    funs=Fun{S,T}[Mk for Mk in M]
    ProductFun{S,V,ProductSpace{S,V},T}(funs,ProductSpace(S[space(fun) for fun in funs],dy))
end

## Adaptive construction

function ProductFun{S<:UnivariateSpace,V<:UnivariateSpace}(f::Function,sp::AbstractProductSpace{Tuple{S,V}};tol=100eps())
    for n = 50:100:5000
        X = coefficients(ProductFun(f,sp,n,n;tol=tol))
        if size(X,1)<n && size(X,2)<n
            return ProductFun(X,sp;tol=tol)
        end
    end
    warn("Maximum grid size of ("*string(5000)*","*string(5000)*") reached")
    ProductFun(f,sp,5000,5000;tol=tol,chopping=true)
end

## ProductFun values to coefficients

function ProductFun(f::Function,S::AbstractProductSpace,M::Integer,N::Integer;tol=100eps())
    xy = checkpoints(S)
    T = promote_type(eltype(f(first(xy)...)),eltype(S))
    ptsx,ptsy=points(S,M,N)
    vals=T[f(ptsx[k,j],ptsy[k,j]) for k=1:size(ptsx,1), j=1:size(ptsx,2)]
    ProductFun(transform!(S,vals),S;tol=tol,chopping=true)
end
ProductFun(f::Function,S::TensorSpace) = ProductFun(LowRankFun(f,S))

ProductFun(f,dx::Space,dy::Space)=ProductFun(f,TensorSpace(dx,dy))

## Domains promoted to Spaces

ProductFun(f::Function,D::BivariateDomain,M::Integer,N::Integer)=ProductFun(f,Space(D),M,N)
ProductFun(f,d::Domain)=ProductFun(f,Space(d))
ProductFun(f,dx::UnivariateDomain,dy::UnivariateDomain)=ProductFun(f,Space(dx),Space(dy))
ProductFun(f::Function) = ProductFun(f,Interval(),Interval())

## Conversion from other 2D Funs

ProductFun(f::LowRankFun)=ProductFun(coefficients(f),space(f,1),space(f,2))
ProductFun{S<:AbstractProductSpace}(f::Fun{S})=ProductFun(coefficientmatrix(f),space(f))

## Conversion to other ProductSpaces with the same coefficients

ProductFun(f::ProductFun,sp::TensorSpace)=space(f)==sp?f:ProductFun(coefficients(f,sp),sp)
ProductFun{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS},sp::ProductDomain)=ProductFun(f,Space(sp))

function ProductFun(f::ProductFun,sp::AbstractProductSpace)
    u=Array(Fun{typeof(columnspace(sp,1)),eltype(f)},length(f.coefficients))

    for k=1:length(f.coefficients)
        u[k]=Fun(f.coefficients[k],columnspace(sp,k))
    end

    ProductFun(u,sp)
end

## For specifying spaces by anonymous function

ProductFun(f::Function,SF::Function,T::Space,M::Integer,N::Integer)=ProductFun(f,typeof(SF(1))[SF(k) for k=1:N],T,M)

## Conversion of a constant to a ProductFun

ProductFun(c::Number,sp::BivariateSpace)=ProductFun([Fun(c,columnspace(sp,1))],sp)
ProductFun(f::Fun,sp::BivariateSpace)=ProductFun([Fun(f,columnspace(sp,1))],sp)

## Utilities

function funlist2coefficients{S,T}(f::Vector{Fun{S,T}})
    A=zeros(T,mapreduce(ncoefficients,max,f),length(f))
    for k=1:length(f)
        A[1:ncoefficients(f[k]),k]=f[k].coefficients
    end
    A
end

function pad{S,V,SS,T}(f::ProductFun{S,V,SS,T},n::Integer,m::Integer)
    ret=Array(Fun{S,T},m)
    cm=min(length(f.coefficients),m)
    for k=1:cm
        ret[k]=pad(f.coefficients[k],n)
    end

    for k=cm+1:m
        ret[k]=zero(T,columnspace(f,k))
    end
    ProductFun{S,V,SS,T}(ret,f.space)
end

function pad!{S,V,SS,T}(f::ProductFun{S,V,SS,T},::Colon,m::Integer)
    cm=length(f.coefficients)
    resize!(f.coefficients,m)

    for k=cm+1:m
        f.coefficients[k]=zero(T,columnspace(f,k))
    end
    f
end

coefficients(f::ProductFun)=funlist2coefficients(f.coefficients)

function coefficients(f::ProductFun,ox::Space,oy::Space)
    T=eltype(f)
    m=size(f,1)
    B=Array(T,m,length(f.coefficients))
    # convert in x direction
    #TODO: adaptively grow in x?
    for k=1:length(f.coefficients)
        B[:,k]=pad!(coefficients(f.coefficients[k],ox),m)
    end

    # convert in y direction
    for k=1:size(B,1)
        ccfs=coefficients(vec(B[k,:]),space(f,2),oy)
        if length(ccfs)>size(B,2)
            B=pad(B,size(B,1),length(ccfs))
        end
        B[k,1:length(ccfs)]=ccfs
        for j=length(ccfs)+1:size(B,2)
            B[k,j]=zero(T)
        end
    end

    B
end

@compat (f::ProductFun)(x,y) = evaluate(f,x,y)

coefficients(f::ProductFun,ox::TensorSpace) = coefficients(f,ox[1],ox[2])

values{S,V,SS,T}(f::ProductFun{S,V,SS,T}) = itransform!(space(f),coefficients(f))

vecpoints{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS},k) = points(f.space[k],size(f,k))

space(f::ProductFun) = f.space
space(f::ProductFun,k) = space(space(f),k)
columnspace(f::ProductFun,k) = columnspace(space(f),k)

domain(f::ProductFun) = domain(f.space)
#domain(f::ProductFun,k)=domain(f.space,k)

canonicalevaluate{S,V,SS,T}(f::ProductFun{S,V,SS,T},x::Number,::Colon) =
    Fun(T[fc(x) for fc in f.coefficients],space(f,2))
canonicalevaluate(f::ProductFun,x::Number,y::Number) = canonicalevaluate(f,x,:)(y)
canonicalevaluate{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS},x::Colon,y::Number) =
    evaluate(f.',y,:)  # doesn't make sense For general product fon without specifying space

canonicalevaluate(f::ProductFun,xx::AbstractVector,yy::AbstractVector) =
    hcat([evaluate(f,x,:)(yy) for x in xx]...).'

evaluate(f::ProductFun,x,y) = canonicalevaluate(f,tocanonical(f,x,y)...)
evaluate(f::ProductFun,x) = evaluate(f,x...)

*{F<:ProductFun}(c::Number,f::F) = F(c*f.coefficients,f.space)
*(f::ProductFun,c::Number) = c*f

function chop{S}(f::ProductFun{S},es...)
    kend=size(f,2)
    if kend > 1 while isempty(chop(f.coefficients[kend].coefficients,es...)) kend-=1 end end
    ret=Fun{S,eltype(f)}[Fun(chop(f.coefficients[k].coefficients,es...),space(f.coefficients[k])) for k=1:max(kend,1)]

    typeof(f)(ret,f.space)
end

##TODO: following assumes f is never changed....maybe should be deepcopy?
function +{F<:ProductFun}(f::F,c::Number)
    cfs=copy(f.coefficients)
    cfs[1]+=c
    F(cfs,f.space)
end
+(c::Number,f::ProductFun) = f+c
-(f::ProductFun,c::Number) = f+(-c)
-(c::Number,f::ProductFun) = c+(-f)

function +(f::ProductFun,g::ProductFun)
    if f.space == g.space
        if size(f,2) >= size(g,2)
            @assert f.space==g.space
            cfs = copy(f.coefficients)
            for k=1:size(g,2)
                cfs[k]+=g.coefficients[k]
            end

            ProductFun(cfs,f.space)
        else
            g+f
        end
    else
        s=conversion_type(f.space,g.space)
        ProductFun(f,s)+ProductFun(g,s)
    end
end

-(f::ProductFun) = (-1)*f
-(f::ProductFun,g::ProductFun) = f+(-g)

*(B::Fun,f::ProductFun) = ProductFun(map(c->B*c,f.coefficients),space(f))
*(f::ProductFun,B::Fun) = (B*f.').'

LowRankFun{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS}) = LowRankFun(f.coefficients,space(f,2))
LowRankFun(f::Fun) = LowRankFun(ProductFun(f))

function differentiate{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS},j::Integer)
    if j==1
        df=map(differentiate,f.coefficients)
        ProductFun(df,space(first(df)),space(f,2))
    else
        differentiate(f.',1).'
    end
end

# If the transpose of the space exists, then the transpose of the ProductFun exists
Base.transpose{S,V,SS,T}(f::ProductFun{S,V,SS,T}) =
    ProductFun(transpose(coefficients(f)),transpose(space(f)))

for op in (:(Base.sin),:(Base.cos))
    @eval ($op)(f::ProductFun) =
        Fun(transform!(space(f),$op(values(pad(f,size(f,1)+20,size(f,2))))),space(f))
end

.^(f::ProductFun,k::Integer) =
    Fun(transform!(space(f),values(pad(f,size(f,1)+20,size(f,2))).^k),space(f))

for op = (:(Base.real),:(Base.imag),:(Base.conj))
    @eval ($op){S,V<:Space{RealBasis},SS<:TensorSpace}(f::ProductFun{S,V,SS}) =
        ProductFun(map($op,f.coefficients),space(f))
end

#For complex bases
Base.real{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS}) =
    real(ProductFun(real(u.coefficients),space(u)).').'-imag(ProductFun(imag(u.coefficients),space(u)).').'
Base.imag{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS}) =
    real(ProductFun(imag(u.coefficients),space(u)).').'+imag(ProductFun(real(u.coefficients),space(u)).').'

## Call LowRankFun version
# TODO: should cumsum and integrate return TensorFun or lowrankfun?
for op in (:(Base.sum),:(Base.cumsum),:integrate)
    @eval $op{S,V,SS<:TensorSpace}(f::ProductFun{S,V,SS},n...) = $op(LowRankFun(f),n...)
end

## ProductFun transform

# function transform{ST<:Space,N<:Number}(::Type{N},S::Vector{ST},T::Space,V::Matrix)
#     @assert length(S)==size(V,2)
#     # We assume all S spaces have same domain/points
#     C=Array(N,size(V)...)
#     for k=1:size(V,1)
#         C[k,:]=transform(T,vec(V[k,:]))
#     end
#     for k=1:size(C,2)
#         C[:,k]=transform(S[k],C[:,k])
#     end
#     C
# end
# transform{ST<:Space,N<:Real}(S::Vector{ST},T::Space{Float64},V::Matrix{N})=transform(Float64,S,T,V)
# transform{ST<:Space}(S::Vector{ST},T::Space,V::Matrix)=transform(Complex{Float64},S,T,V)

for op in (:tocanonical,:fromcanonical)
    @eval $op(f::ProductFun,x...) = $op(space(f),x...)
end


text:
# setup paths -
_PATH_TO_ROOT = pwd()
_PATH_TO_SRC = joinpath(_PATH_TO_ROOT, "src")
_PATH_TO_DATA = joinpath(_PATH_TO_ROOT, "data")
_PATH_TO_CONFIG = joinpath(_PATH_TO_ROOT, "config")

# setup my package -
import Pkg
Pkg.activate(_PATH_TO_ROOT)
Pkg.update()
Pkg.instantiate()

# load external packages -
using TOML
using DataFrames
using CSV
using Plots
using DelimitedFiles
using Optim

# load my codes -
my_code_array = [
    "Data.jl", "Model.jl", "Optimize.jl"
]
for my_code_file in my_code_array
    path_to_code_file = joinpath(_PATH_TO_SRC, my_code_file)
    include(path_to_code_file)
end


text:
using Distributed
addprocs(3);
@everywhere include("./src/MCPhylo.jl")
@everywhere using .MCPhylo

mt_st, df_st = make_tree_with_data("notebook/data-st-64-110.paps.nex"); # load your own nexus file
mt_ie, df_ie = make_tree_with_data("notebook/data-ie-42-208.paps.nex"); # load your own nexus file
mt_aa, df_aa = make_tree_with_data("notebook/data-aa-58-200.paps.nex"); # load your own nexus file
mt_an, df_an = make_tree_with_data("notebook/data-an-45-210.paps.nex"); # load your own nexus file
mt_pn, df_pn = make_tree_with_data("notebook/data-pn-67-183.paps.nex"); # load your own nexus file
mt_st2 = deepcopy(mt_st)
mt_ie2 = deepcopy(mt_ie)
mt_an2 = deepcopy(mt_an)
mt_pn2 = deepcopy(mt_pn)
mt_aa2 = deepcopy(mt_aa)
mt_st3 = deepcopy(mt_st)
mt_ie3 = deepcopy(mt_ie)
mt_aa3 = deepcopy(mt_aa)
mt_an3 = deepcopy(mt_an)
mt_pn3 = deepcopy(mt_pn)

randomize!(mt_st2)
randomize!(mt_ie2)
randomize!(mt_aa2)
randomize!(mt_an2)
randomize!(mt_pn2)

randomize!(mt_st3)
randomize!(mt_ie3)
randomize!(mt_aa3)
randomize!(mt_an3)
randomize!(mt_pn3)

my_data = Dict{Symbol, Any}(
  :mtree => [mt_st, mt_ie, mt_aa, mt_an, mt_pn],
  :df => [df_st, df_ie, df_aa, df_an, df_pn],
  :nnodes => [size(df_st)[3],size(df_ie)[3],size(df_an)[3]],
  :nbase => [size(df_st)[1],size(df_ie)[1],size(df_an)[1]],
  :nsites => [size(df_st)[2],size(df_ie)[2],size(df_an)[2]],
);

# model setup
model =  Model(
    df_ie = Stochastic(3, (tree_ie, pi_ie, rates) -> PhyloDist(tree_ie, pi_ie, [1.0], rates[1:4], Restriction), false, false),
    df_st = Stochastic(3, (tree_st, pi_st, rates) -> PhyloDist(tree_st, pi_st, [1.0], rates[5:8], Restriction), false, false),
    df_aa = Stochastic(3, (tree_aa, pi_aa, rates) -> PhyloDist(tree_aa, pi_aa, [1.0], rates[9:12], Restriction), false, false),
    df_an = Stochastic(3, (tree_an, pi_an, rates) -> PhyloDist(tree_an, pi_an, [1.0], rates[13:16], Restriction), false, false),
    df_pn = Stochastic(3, (tree_pn, pi_pn, rates) -> PhyloDist(tree_pn, pi_pn, [1.0], rates[17:20], Restriction), false, false),
    pi_ie = Stochastic(1, ()-> Dirichlet(2,1)),
    pi_st = Stochastic(1, ()-> Dirichlet(2,1)),
    pi_aa = Stochastic(1, ()-> Dirichlet(2,1)),
    pi_an = Stochastic(1, ()-> Dirichlet(2,1)),
    pi_pn = Stochastic(1, ()-> Dirichlet(2,1)),
    tree_ie = Stochastic(Node(), (a,b,c,d) -> CompoundDirichlet(a[1],b[1],c[1],d[1]), true),
    tree_st = Stochastic(Node(), (a,b,c,d) -> CompoundDirichlet(a[2],b[2],c[2],d[2]), true),
    tree_aa = Stochastic(Node(), (a,b,c,d) -> CompoundDirichlet(a[3],b[3],c[3],d[3]), true),
    tree_an = Stochastic(Node(), (a,b,c,d) -> CompoundDirichlet(a[4],b[4],c[4],d[4]), true),
    tree_pn = Stochastic(Node(), (a,b,c,d) -> CompoundDirichlet(a[5],b[5],c[5],d[5]), true),
    a = Stochastic(1, ()->Gamma()),
    b = Stochastic(1, ()->Gamma()),
    c = Stochastic(1, ()->Beta(2,2)),
    d = Stochastic(1, ()->Gamma()),
    rates = Logical(1, (αs) -> vcat(discrete_gamma_rates.(αs, αs, 4)...),false),
    αs = Stochastic(1, () -> Gamma(), true),
     )

# intial model values
inits = [ Dict{Symbol, Union{Any, Real}}(
    :tree_ie => mt_ie,
    :tree_st => mt_st,
    :tree_pn => mt_pn,
    :tree_aa => mt_aa,
    :tree_an => mt_an,
    :pi_ie=> rand(Dirichlet(2, 1)),
    :pi_st=> rand(Dirichlet(2, 1)),
    :pi_an=> rand(Dirichlet(2, 1)),
    :pi_aa=> rand(Dirichlet(2, 1)),
    :pi_pn=> rand(Dirichlet(2, 1)),
    :df_ie => my_data[:df][2],
    :df_st => my_data[:df][1],
    :df_aa => my_data[:df][3],
    :df_an => my_data[:df][4],
    :df_pn => my_data[:df][5],
    :αs => rand(5),
    :a => rand(5),
    :b => rand(5),
    :c => rand(5),
    :d => rand(5),
    ),
    Dict{Symbol, Union{Any, Real}}(
    :tree_ie => mt_ie2,
    :tree_st => mt_st2,
    :tree_pn => mt_pn2,
    :tree_aa => mt_aa2,
    :tree_an => mt_an2,
    :pi_ie=> rand(Dirichlet(2, 1)),
    :pi_st=> rand(Dirichlet(2, 1)),
    :pi_an=> rand(Dirichlet(2, 1)),
    :pi_aa=> rand(Dirichlet(2, 1)),
    :pi_pn=> rand(Dirichlet(2, 1)),
    :df_ie => my_data[:df][2],
    :df_st => my_data[:df][1],
    :df_aa => my_data[:df][3],
    :df_an => my_data[:df][4],
    :df_pn => my_data[:df][5],
    :αs => rand(5),
    :a => rand(5),
    :b => rand(5),
    :c => rand(5),
    :d => rand(5),
    ),
    Dict{Symbol, Union{Any, Real}}(
    :tree_ie => mt_ie3,
    :tree_st => mt_st3,
    :tree_pn => mt_pn3,
    :tree_aa => mt_aa3,
    :tree_an => mt_an3,
    :pi_ie=> rand(Dirichlet(2, 1)),
    :pi_st=> rand(Dirichlet(2, 1)),
    :pi_an=> rand(Dirichlet(2, 1)),
    :pi_aa=> rand(Dirichlet(2, 1)),
    :pi_pn=> rand(Dirichlet(2, 1)),
    :df_ie => my_data[:df][2],
    :df_st => my_data[:df][1],
    :df_aa => my_data[:df][3],
    :df_an => my_data[:df][4],
    :df_pn => my_data[:df][5],
    :αs => rand(5),
    :a => rand(5),
    :b => rand(5),
    :c => rand(5),
    :d => rand(5),
    ),
    ]

scheme = [PNUTS(:tree_ie, target=0.8, targetNNI=7),
          PNUTS(:tree_st, target=0.8, targetNNI=7),
          PNUTS(:tree_aa, target=0.8, targetNNI=7),
          PNUTS(:tree_an, target=0.8, targetNNI=7),
          PNUTS(:tree_pn, target=0.8, targetNNI=7),
          SliceSimplex(:pi_ie),
          SliceSimplex(:pi_an),
          SliceSimplex(:pi_st),
          SliceSimplex(:pi_aa),
          SliceSimplex(:pi_pn),
          Slice([:αs, :a, :b, :c, :d], 1.0, Multivariate),
          ]

setsamplers!(model, scheme);

# do the mcmc simmulation. if trees=true the trees are stored and can later be
# flushed ot a file output.
sim = mcmc(model, my_data, inits, 25, burnin=5,thin=1, chains=1, trees=true)
sim = mcmc(sim, 10, trees=true)
to_file(sim, "testMult_")


text:
# Statistics
# ==========

mutable struct WorkerStats
    # the number of calls emitted
    ncalls::Int

    # the number of calls finished
    nfinished::Int

    function WorkerStats()
        return new(0, 0)
    end
end

mutable struct Stats
    # statistics state {:uninitialized, :initialized, :finalized}
    state::Symbol

    # the number of workers
    nworkers::Int

    # stats of each worker
    workerstats::Dict{Int,WorkerStats}

    # initialized time
    inittime::DateTime

    # initialized time (nano second)
    inittime_ns::UInt64

    # finalized time (nano second)
    fintime_ns::UInt64

    # total call time (nano second)
    calltime_ns::UInt64

    # total worker select time (nano second)
    selecttime_ns::UInt64

    # total remotecall time (nano second)
    remotecalltime_ns::UInt64

    # total waitargs time (nano second)
    waitargstime_ns::UInt64

    # teardown time (nano second)
    teardowntime_ns::UInt64

    function Stats()
        return new(
            :uninitialized, 0, Dict{Int,WorkerStats}(),
            DateTime(), 0, 0, 0, 0, 0, 0, 0)
    end
end

function Base.show(io::IO, stats::Stats)
    println(io, "Statistics:")
    println(io, "  state=$(stats.state)")
    if stats.state == :uninitialized
        return
    end
    println(io, "  nworkers=$(stats.nworkers)")
    println(io, "  worker stats:")
    for worker in sort!(collect(keys(stats.workerstats)))
        workerstats = stats.workerstats[worker]
        println(io, "    worker=$(worker) ncalls=$(workerstats.ncalls) nfinished=$(workerstats.nfinished)")
    end
    println(io, "  initialization: ", stats.inittime)
    if stats.state == :finalized
        total = stats.fintime_ns - stats.inittime_ns
        println(io, "  total time:   $(readable_time(total)) (100.00%)")
        println(io, "    call:       $(readable_time(stats.calltime_ns)) ($(readable_ratio(stats.calltime_ns, total)))")
        println(io, "    select:     $(readable_time(stats.selecttime_ns)) ($(readable_ratio(stats.selecttime_ns, total)))")
        println(io, "    remotecall: $(readable_time(stats.remotecalltime_ns)) ($(readable_ratio(stats.remotecalltime_ns, total)))")
        println(io, "    waitargs:   $(readable_time(stats.waitargstime_ns)) ($(readable_ratio(stats.waitargstime_ns, total)))")
        println(io, "    teardown:   $(readable_time(stats.teardowntime_ns)) ($(readable_ratio(stats.teardowntime_ns, total)))")
    end
end

function readable_time(ns)
    scale = 9
    for unit in ["s", "ms", "μs"]
        if ns ≥ 10^scale
            return @sprintf("%.2f%s", ns / 10^scale, unit)
        end
        scale -= 3
    end
    return @sprintf("%.2fns", ns)
end

function readable_ratio(num, den)
    return @sprintf("%.2f%%", num / den * 100)
end

function initstats!(stats, workers)
    stats.inittime = now()
    stats.inittime_ns = time_ns()
    stats.nworkers = length(workers)
    for worker in workers
        stats.workerstats[worker] = WorkerStats()
    end
    stats.state = :initialized
    return stats
end

function finstats!(stats)
    stats.fintime_ns = time_ns()
    stats.state = :finalized
    return stats
end

macro tic()
    quote
        __start__ = time_ns()
    end |> esc
end

macro toc(metric)
    quote
        $(metric) += time_ns() - __start__
    end |> esc
end


text:
include("updateA!.jl")
include("normcol_equal.jl")

# return inv(τ*A*A' + λ*B*B' + γ*I)
function getinv(τ::Float64, λ::Float64, γ::Float64, A::Matrix{Float64}, B::Matrix{Float64})
    C = A*A'
    D = B*B'
    Dim = size(C, 1)

    for j=1:Dim, i=1:Dim
        if i==j
            @inbounds C[i,j] = τ*C[i,j] + λ*D[i,j] + γ
        else
            @inbounds C[i,j] = τ*C[i,j] + λ*D[i,j]
        end
    end

    inv(C)
end

function initialization(Data::Matrix{Float64}, Label::Matrix{Int64}, DictSize::Int64, τ::Float64, λ::Float64, γ::Float64, DEMO=false)
    # In this initialization function, we do the following things:
    # 1. Random initialization of dictionary pair D and P for each class
    # 2. Compute the class-specific inverse matrix used in Eq. (10)
    # 3. Compute matrix class-specific code matrix A by Eq. (8)
    #    with the random initialized D and P
    #
    # The randn seeds are set to make sure the results in our paper are
    # reproducible. The randn seed setting can be removed, our algorithm is
    # not sensitive to the initialization of D and P. In most cases, different
    # initialization will lead to the same recognition accuracy on a wide range
    # of testing databases.
    ClassNum   = maximum(Label)
    Dim        = size(Data, 1)

    DataMat    = Array(Any, ClassNum)
    DictMat    = Array(Any, ClassNum)
    P          = Array(Any, ClassNum)
    DataInvMat = Array(Any, ClassNum)
    CoefMat    = Array(Any, ClassNum)

    for i=1:ClassNum
        @inbounds TempData   = Data[:, find(Label .== i)]
        @inbounds DataMat[i] = TempData

        DEMO && srand(i)
        TempRand = randn(Dim, DictSize)
        TempRand = normcol_equal(TempRand)
        @inbounds DictMat[i] = TempRand

        DEMO && srand(2i)
        TempRand = randn(Dim, DictSize)
        TempRand = normcol_equal(TempRand)
        @inbounds P[i]       = TempRand'

        @inbounds TempDataC = Data[:, find(Label .!=i)]
        @inbounds DataInvMat[i] = getinv(τ, λ, γ, TempData, TempDataC)
    end

    updateA!(CoefMat, DictMat, DataMat, P, τ, DictSize)

    DataMat, DictMat, P, DataInvMat, CoefMat, ClassNum, Dim
end


text:
# tests String value choice points, i.e. using choose(StringDataType, regex) construct

@generator SCWildcardGen begin # prefix SC (for String value Choice points) to avoid type name clashes
    start() = choose(String, ".")
end

@generator SCQuantifiersGen begin
    start() = choose(String, "a?b+c*d{4}e{5,6}f{7,}g{8,8}")
end

@generator SCAlternationGen begin
    start() = choose(String, "foo|bar|baz")
end

@generator SCBracketsGen begin
    start() = choose(String, "a[uvw][x-z0-3]b")
end

@generator SCParenthesesGen begin
    start() = choose(String, "a(bc|de)")
end

@generator SCClassesGen begin
    start() = choose(String, "\\s\\S\\d\\D\\w\\W")
end

@generator SCEscapesGen begin
    start() = choose(String, "\\.\\[\\]\\|\\?\\+\\*\\\\")
end

@testset "string value choice points" begin

	@testset "choose(String) using regex containing wilcard" begin

		gn = SCWildcardGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^.$", td)
		    @mtest_values_vary td
		end

	end

	@testset "choose(String) using regex containing quantifiers" begin

		gn = SCQuantifiersGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^a?b+c*d{4}e{5,6}f{7,}g{8,8}$", td)
		    @mtest_values_include [0,1] count(x->x=='a',td)
		    @mtest_values_include [1,2,3] count(x->x=='b',td)
		    @mtest_values_include [0,1,2] count(x->x=='c',td)
		    @mtest_values_include [4,] count(x->x=='d',td)
		    @mtest_values_include [5,6] count(x->x=='e',td)
		    @mtest_values_include [7,8,9] count(x->x=='f',td)
		    @mtest_values_are [8] count(x->x=='g',td)
		end

	end

	@testset "choose(String) using regex containing alternation" begin

		gn = SCAlternationGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^foo|bar|baz$", td)
		    @mtest_values_are ["foo","bar","baz"] td
		end

	end

	@testset "choose(String) using regex containing bracket" begin

		gn = SCBracketsGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^a[uvw][x-z0-3]b$", td)
		    @mtest_values_are ['u','v','w'] td[2]
		    @mtest_values_are ['x','y','z','0','1','2','3'] td[3]
		end

	end

	@testset "choose(String) using regex containing parentheses" begin

		gn = SCParenthesesGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^a(bc+|de+)$", td)
		    @mtest_values_are ['b','d'] td[2]
		    @mtest_values_are ['c','e'] td[3]
		end

	end

	@testset "choose(String) using regex containing classes" begin

		gn = SCClassesGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^\s\S\d\D\w\W$", td)
		    @mtest_values_vary td
		end

	end

	@testset "choose(String) using regex that escapes metacharacters" begin

		gn = SCEscapesGen()

		@mtestset reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn)
		    @test typeof(td) == String
		    @test ismatch(r"^\.\[\]\|\?\+\*\\$", td)
		end

	end

end


text:
function Circle(arg0::Circle)
    return Circle((Circle,), arg0)
end

function Circle(arg0::S2Point, arg1::S2Point, arg2::jdouble)
    return Circle((S2Point, S2Point, jdouble), arg0, arg1, arg2)
end

function Circle(arg0::Vector3D, arg1::jdouble)
    return Circle((Vector3D, jdouble), arg0, arg1)
end

function copy_self(obj::Circle)
    return jcall(obj, "copySelf", Circle, ())
end

function empty_hyperplane(obj::Circle)
    return jcall(obj, "emptyHyperplane", SubCircle, ())
end

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

function get_arc(obj::Circle, arg0::S2Point, arg1::S2Point)
    return jcall(obj, "getArc", Arc, (S2Point, S2Point), arg0, arg1)
end

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

function get_inside_arc(obj::Circle, arg0::Circle)
    return jcall(obj, "getInsideArc", Arc, (Circle,), arg0)
end

function get_offset(obj::Circle, arg0::Point)
    return jcall(obj, "getOffset", jdouble, (Point,), arg0)
end

function get_offset(obj::Circle, arg0::Vector3D)
    return jcall(obj, "getOffset", jdouble, (Vector3D,), arg0)
end

function get_phase(obj::Circle, arg0::Vector3D)
    return jcall(obj, "getPhase", jdouble, (Vector3D,), arg0)
end

function get_point_at(obj::Circle, arg0::jdouble)
    return jcall(obj, "getPointAt", Vector3D, (jdouble,), arg0)
end

function get_pole(obj::Circle)
    return jcall(obj, "getPole", Vector3D, ())
end

function get_reverse(obj::Circle)
    return jcall(obj, "getReverse", Circle, ())
end

function get_tolerance(obj::Circle)
    return jcall(obj, "getTolerance", jdouble, ())
end

function get_transform(::Type{Circle}, arg0::Rotation)
    return jcall(Circle, "getTransform", Transform, (Rotation,), arg0)
end

function get_x_axis(obj::Circle)
    return jcall(obj, "getXAxis", Vector3D, ())
end

function get_y_axis(obj::Circle)
    return jcall(obj, "getYAxis", Vector3D, ())
end

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

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

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

function project(obj::Circle, arg0::Point)
    return jcall(obj, "project", Point, (Point,), arg0)
end

function reset(obj::Circle, arg0::Vector3D)
    return jcall(obj, "reset", void, (Vector3D,), arg0)
end

function revert_self(obj::Circle)
    return jcall(obj, "revertSelf", void, ())
end

function same_orientation_as(obj::Circle, arg0::Hyperplane)
    return jcall(obj, "sameOrientationAs", jboolean, (Hyperplane,), arg0)
end

function to_space(obj::Circle, arg0::Point)
    return jcall(obj, "toSpace", S2Point, (Point,), arg0)
end

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

function to_sub_space(obj::Circle, arg0::Point)
    return jcall(obj, "toSubSpace", S1Point, (Point,), arg0)
end

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

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

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

function whole_hyperplane(obj::Circle)
    return jcall(obj, "wholeHyperplane", SubCircle, ())
end

function whole_space(obj::Circle)
    return jcall(obj, "wholeSpace", SphericalPolygonsSet, ())
end


text:
#############################
### 2D Coordinate systems ###
#############################
"""
`Polar{T,A}(r::T, θ::A)` - 2D polar coordinates
"""
struct Polar{T,A}
    r::T
    θ::A

    Polar{T, A}(r, θ) where {T, A} = new(r, θ)
end

function Polar(r, θ)
    r2, θ2 = promote(r, θ)

    return Polar{typeof(r2), typeof(θ2)}(r2, θ2)
end

Base.show(io::IO, x::Polar) = print(io, "Polar(r=$(x.r), θ=$(x.θ) rad)")
Base.isapprox(p1::Polar, p2::Polar; kwargs...) = isapprox(p1.r, p2.r; kwargs...) && isapprox(p1.θ, p2.θ; kwargs...)

"`PolarFromCartesian()` - transformation from `AbstractVector` of length 2 to `Polar` type"
struct PolarFromCartesian <: Transformation; end
"`CartesianFromPolar()` - transformation from `Polar` type to `SVector{2}` type"
struct CartesianFromPolar <: Transformation; end

Base.show(io::IO, trans::PolarFromCartesian) = print(io, "PolarFromCartesian()")
Base.show(io::IO, trans::CartesianFromPolar) = print(io, "CartesianFromPolar()")

function (::PolarFromCartesian)(x::AbstractVector)
    length(x) == 2 || error("Polar transform takes a 2D coordinate")

    Polar(hypot(x[1], x[2]), atan(x[2], x[1]))
end

function transform_deriv(::PolarFromCartesian, x::AbstractVector)
    length(x) == 2 || error("Polar transform takes a 2D coordinate")

    r = hypot(x[1], x[2])
    f = x[2] / x[1]
    c = one(eltype(x))/(x[1]*(one(eltype(x)) + f*f))
    @SMatrix [ x[1]/r    x[2]/r ;
              -f*c       c      ]
end
transform_deriv_params(::PolarFromCartesian, x::AbstractVector) = error("PolarFromCartesian has no parameters")

function (::CartesianFromPolar)(x::Polar)
    s,c = sincos(x.θ)
    SVector(x.r * c, x.r * s)
end
function transform_deriv(::CartesianFromPolar, x::Polar)
    sθ, cθ = sincos(x.θ)
    @SMatrix [cθ  -x.r*sθ ;
              sθ   x.r*cθ ]
end
transform_deriv_params(::CartesianFromPolar, x::Polar) = error("CartesianFromPolar has no parameters")

Base.inv(::PolarFromCartesian) = CartesianFromPolar()
Base.inv(::CartesianFromPolar) = PolarFromCartesian()

compose(::PolarFromCartesian, ::CartesianFromPolar) = IdentityTransformation()
compose(::CartesianFromPolar, ::PolarFromCartesian) = IdentityTransformation()

# For convenience
Base.convert(::Type{Polar}, v::AbstractVector) = PolarFromCartesian()(v)
@inline Base.convert(::Type{V}, p::Polar) where {V <: AbstractVector} = convert(V, CartesianFromPolar()(p))
@inline Base.convert(::Type{V}, p::Polar) where {V <: StaticVector} = convert(V, CartesianFromPolar()(p))

#############################
### 3D Coordinate Systems ###
#############################
"""
Spherical(r, θ, ϕ) - 3D spherical coordinates

There are many Spherical coordinate conventions and this library uses a somewhat exotic one.
Given a vector `v` with Cartesian coordinates `xyz`, let `v_xy = [x,y,0]` be the orthogonal projection of `v` on the `xy` plane.

* `r` is the radius. It is given by `norm(v, 2)`.
* `θ` is the azimuth. It is the angle from the x-axis to `v_xy`
* `ϕ` is the latitude. It is the angle from `v_xy` to `v`.

```jldoctest
julia> using CoordinateTransformations

julia> v = randn(3);

julia> sph = SphericalFromCartesian()(v);

julia> r = sph.r; θ=sph.θ; ϕ=sph.ϕ;

julia> v ≈ [r * cos(θ) * cos(ϕ), r * sin(θ) * cos(ϕ), r*sin(ϕ)]
true
"""
struct Spherical{T,A}
    r::T
    θ::A
    ϕ::A

    Spherical{T, A}(r, θ, ϕ) where {T, A} = new(r, θ, ϕ)
end

function Spherical(r, θ, ϕ)
    r2, θ2, ϕ2 = promote(r, θ, ϕ)

    return Spherical{typeof(r2), typeof(θ2)}(r2, θ2, ϕ2)
end

Base.show(io::IO, x::Spherical) = print(io, "Spherical(r=$(x.r), θ=$(x.θ) rad, ϕ=$(x.ϕ) rad)")
Base.isapprox(p1::Spherical, p2::Spherical; kwargs...) = isapprox(p1.r, p2.r; kwargs...) && isapprox(p1.θ, p2.θ; kwargs...) && isapprox(p1.ϕ, p2.ϕ; kwargs...)

"""
Cylindrical(r, θ, z) - 3D cylindrical coordinates
"""
struct Cylindrical{T,A}
    r::T
    θ::A
    z::T

    Cylindrical{T, A}(r, θ, z) where {T, A} = new(r, θ, z)
end

function Cylindrical(r, θ, z)
    r2, θ2, z2 = promote(r, θ, z)

    return Cylindrical{typeof(r2), typeof(θ2)}(r2, θ2, z2)
end

Base.show(io::IO, x::Cylindrical) = print(io, "Cylindrical(r=$(x.r), θ=$(x.θ) rad, z=$(x.z))")
Base.isapprox(p1::Cylindrical, p2::Cylindrical; kwargs...) = isapprox(p1.r, p2.r; kwargs...) && isapprox(p1.θ, p2.θ; kwargs...) && isapprox(p1.z, p2.z; kwargs...)

"`SphericalFromCartesian()` - transformation from 3D point to `Spherical` type"
struct SphericalFromCartesian <: Transformation; end
"`CartesianFromSpherical()` - transformation from `Spherical` type to `SVector{3}` type"
struct CartesianFromSpherical <: Transformation; end
"`CylindricalFromCartesian()` - transformation from 3D point to `Cylindrical` type"
struct CylindricalFromCartesian <: Transformation; end
"`CartesianFromCylindrical()` - transformation from `Cylindrical` type to `SVector{3}` type"
struct CartesianFromCylindrical <: Transformation; end
"`CylindricalFromSpherical()` - transformation from `Spherical` type to `Cylindrical` type"
struct CylindricalFromSpherical <: Transformation; end
"`SphericalFromCylindrical()` - transformation from `Cylindrical` type to `Spherical` type"
struct SphericalFromCylindrical <: Transformation; end

Base.show(io::IO, trans::SphericalFromCartesian) = print(io, "SphericalFromCartesian()")
Base.show(io::IO, trans::CartesianFromSpherical) = print(io, "CartesianFromSpherical()")
Base.show(io::IO, trans::CylindricalFromCartesian) = print(io, "CylindricalFromCartesian()")
Base.show(io::IO, trans::CartesianFromCylindrical) = print(io, "CartesianFromCylindrical()")
Base.show(io::IO, trans::CylindricalFromSpherical) = print(io, "CylindricalFromSpherical()")
Base.show(io::IO, trans::SphericalFromCylindrical) = print(io, "SphericalFromCylindrical()")

# Cartesian <-> Spherical
function (::SphericalFromCartesian)(x::AbstractVector)
    length(x) == 3 || error("Spherical transform takes a 3D coordinate")

    Spherical(hypot(x[1], x[2], x[3]), atan(x[2], x[1]), atan(x[3], hypot(x[1], x[2])))
end
function transform_deriv(::SphericalFromCartesian, x::AbstractVector)
    length(x) == 3 || error("Spherical transform takes a 3D coordinate")
    T = eltype(x)

    r = hypot(x[1], x[2], x[3])
    rxy = hypot(x[1], x[2])
    fxy = x[2] / x[1]
    cxy = one(T)/(x[1]*(one(T) + fxy*fxy))
    f = -x[3]/(rxy*r*r)

    @SMatrix [ x[1]/r   x[2]/r  x[3]/r;
          -fxy*cxy  cxy     zero(T);
           f*x[1]   f*x[2]  rxy/(r*r) ]
end
transform_deriv_params(::SphericalFromCartesian, x::AbstractVector) = error("SphericalFromCartesian has no parameters")

function (::CartesianFromSpherical)(x::Spherical)
    sθ, cθ = sincos(x.θ)
    sϕ, cϕ = sincos(x.ϕ)
    SVector(x.r * cθ * cϕ, x.r * sθ * cϕ, x.r * sϕ)
end
function transform_deriv(::CartesianFromSpherical, x::Spherical{T}) where T
    sθ, cθ = sincos(x.θ)
    sϕ, cϕ = sincos(x.ϕ)
    @SMatrix [cθ*cϕ -x.r*sθ*cϕ -x.r*cθ*sϕ ;
              sθ*cϕ  x.r*cθ*cϕ -x.r*sθ*sϕ ;
              sϕ     zero(T)    x.r * cϕ  ]
end
transform_deriv_params(::CartesianFromSpherical, x::Spherical) = error("CartesianFromSpherical has no parameters")

# Cartesian <-> Cylindrical
function (::CylindricalFromCartesian)(x::AbstractVector)
    length(x) == 3 || error("Cylindrical transform takes a 3D coordinate")

    Cylindrical(hypot(x[1], x[2]), atan(x[2], x[1]), x[3])
end

function transform_deriv(::CylindricalFromCartesian, x::AbstractVector)
    length(x) == 3 || error("Cylindrical transform takes a 3D coordinate")
    T = eltype(x)

    r = hypot(x[1], x[2])
    f = x[2] / x[1]
    c = one(T)/(x[1]*(one(T) + f*f))
    @SMatrix [ x[1]/r   x[2]/r   zero(T) ;
              -f*c      c        zero(T) ;
               zero(T)  zero(T)  one(T)  ]
end
transform_deriv_params(::CylindricalFromCartesian, x::AbstractVector) = error("CylindricalFromCartesian has no parameters")

function (::CartesianFromCylindrical)(x::Cylindrical)
    sθ, cθ = sincos(x.θ)
    SVector(x.r * cθ, x.r * sθ, x.z)
end
function transform_deriv(::CartesianFromCylindrical, x::Cylindrical{T}) where {T}
    sθ, cθ = sincos(x.θ)
    @SMatrix [cθ      -x.r*sθ  zero(T) ;
              sθ       x.r*cθ  zero(T) ;
              zero(T)  zero(T) one(T)  ]
end
transform_deriv_params(::CartesianFromPolar, x::Cylindrical) = error("CartesianFromCylindrical has no parameters")

function (::CylindricalFromSpherical)(x::Spherical)
    sϕ, cϕ = sincos(x.ϕ)
    Cylindrical(x.r*cϕ,x.θ,x.r*sϕ)
end
function transform_deriv(::CylindricalFromSpherical, x::Spherical)
    M1 = transform_deriv(CylindricalFromCartesian(), CartesianFromSpherical()(x))
    M2 = transform_deriv(CartesianFromSpherical(), x)
    return M1*M2
end
transform_deriv_params(::CylindricalFromSpherical, x::Spherical) = error("CylindricalFromSpherical has no parameters")

function (::SphericalFromCylindrical)(x::Cylindrical)
    Spherical(hypot(x.r,x.z),x.θ,atan(x.z,x.r))
end
function transform_deriv(::SphericalFromCylindrical, x::Cylindrical)
    M1 = transform_deriv(SphericalFromCartesian(), CartesianFromCylindrical()(x))
    M2 = transform_deriv(CartesianFromCylindrical(), x)
    return M1*M2
end
transform_deriv_params(::SphericalFromCylindrical, x::Cylindrical) = error("SphericalFromCylindrical has no parameters")

Base.inv(::SphericalFromCartesian)   = CartesianFromSpherical()
Base.inv(::CartesianFromSpherical)   = SphericalFromCartesian()
Base.inv(::CylindricalFromCartesian) = CartesianFromCylindrical()
Base.inv(::CartesianFromCylindrical) = CylindricalFromCartesian()
Base.inv(::CylindricalFromSpherical) = SphericalFromCylindrical()
Base.inv(::SphericalFromCylindrical) = CylindricalFromSpherical()

# Inverse composition
compose(::SphericalFromCartesian,   ::CartesianFromSpherical)   = IdentityTransformation()
compose(::CartesianFromSpherical,   ::SphericalFromCartesian)   = IdentityTransformation()
compose(::CylindricalFromCartesian, ::CartesianFromCylindrical) = IdentityTransformation()
compose(::CartesianFromCylindrical, ::CylindricalFromCartesian) = IdentityTransformation()
compose(::CylindricalFromSpherical, ::SphericalFromCylindrical) = IdentityTransformation()
compose(::SphericalFromCylindrical, ::CylindricalFromSpherical) = IdentityTransformation()

# Cyclic compositions
compose(::SphericalFromCartesian,   ::CartesianFromCylindrical) = SphericalFromCylindrical()
compose(::CartesianFromSpherical,   ::SphericalFromCylindrical) = CartesianFromCylindrical()
compose(::CylindricalFromCartesian, ::CartesianFromSpherical)   = CylindricalFromSpherical()
compose(::CartesianFromCylindrical, ::CylindricalFromSpherical) = CartesianFromSpherical()
compose(::CylindricalFromSpherical, ::SphericalFromCartesian)   = CylindricalFromCartesian()
compose(::SphericalFromCylindrical, ::CylindricalFromCartesian) = SphericalFromCartesian()

# For convenience
Base.convert(::Type{Spherical}, v::AbstractVector) = SphericalFromCartesian()(v)
Base.convert(::Type{Cylindrical}, v::AbstractVector) = CylindricalFromCartesian()(v)

Base.convert(::Type{V}, s::Spherical) where {V <: AbstractVector} = convert(V, CartesianFromSpherical()(s))
Base.convert(::Type{V}, c::Cylindrical) where {V <: AbstractVector} = convert(V, CartesianFromCylindrical()(c))
Base.convert(::Type{V}, s::Spherical) where {V <: StaticVector} = convert(V, CartesianFromSpherical()(s))
Base.convert(::Type{V}, c::Cylindrical) where {V <: StaticVector} = convert(V, CartesianFromCylindrical()(c))

Base.convert(::Type{Spherical}, c::Cylindrical) = SphericalFromCylindrical()(c)
Base.convert(::Type{Cylindrical}, s::Spherical) = CylindricalFromSpherical()(s)


text:
cd(dirname(@__FILE__))

using LinearAlgebra
using Random
using Statistics
using Test

include("../src/SimpleBirthDeathProcess.jl")
using .SimpleBirthDeathProcess

Random.seed!(43314697)

function runtests()
  tests = [
    "common/functions",
    "common/hypergeometric",
    "common/meixner_ortho_poly",
    "simple/transition_probabilities",
    "simple/simulation",
    "simple/gradient",
    "simple/hessian"
  ]

  for t in tests
    f = string(t, ".jl")
    println("Going through tests in '", f, "'...")
    include(f)
    println("PASSED!")
  end

  nothing
end

runtests()


text:
using DataStructures

"""
Represents the tiles in a single player's hand.
"""
struct Hand
    tiles::Vector{SuitedTile}
    hidden::Vector{SuitedTile}
    exposed::Vector{SuitedTile}

    @doc """
        Hand()
        Hand(tiles)

    Create an empty hand, or one that starts off with the `tiles` given (all hidden).

    ## Example

    ```jldoctest
    julia> using MahjongTiles: character, wind, bamboo, Hand

    julia> Hand([character(1), character(2), wind(:north), bamboo(5)])
    🀃🀇🀈🀔
    ```
    """
    function Hand(initial_tiles::AbstractVector{<: Tile} = [])
        hidden_tiles = sort!(initial_tiles)
        return new(hidden_tiles, hidden_tiles, SuitedTile[])
    end
end

@doc """
    push_hidden!(hand, tiles...)
    Base.push!(hand, tiles...)

Add a tile or `tiles` to the hidden part of the hand.
"""
function push_hidden!(hand::Hand, tile...)
    push!(hand.hidden, tile...)
    push!(hand.tiles, tile...)
    sort!(hand.hidden, alg=InsertionSort())
    sort!(hand.tiles, alg=InsertionSort())

    return hand
end
Base.push!(hand::Hand, tile...) = push_hidden!(hand, tile)

@doc """
    push_exposed!(hand, tiles...)

Add a tile or `tiles` to the open/exposed part of the hand.
"""
function push_exposed!(hand::Hand, tile...)
    push!(hand.exposed, tile...)
    push!(hand.tiles, tile...)
    sort!(hand.exposed, alg=InsertionSort())
    sort!(hand.tiles, alg=InsertionSort())

    return hand
end

@doc """
    expose!(hand, tiles...)

Expose/make open the tiles from the hand.
"""
function expose!(hand::Hand, tiles...)

end

Base.first(hand::Hand) = first(hand.tiles)
Base.first(hand::Hand, num::Integer) = first(hand.tiles, num)

Base.iterate(hand::Hand) = iterate(hand.tiles)
Base.iterate(hand::Hand, state) = iterate(hand.tiles, state)

function Base.show(io::IO, mime::MIME"text/plain", hand::Hand)
    for tile in hand.tiles
        show(io, mime, tile)
    end
end

@doc """
    is_pure(hand)

Determine whether the hand is pure (清).
A hand is considered pure if there are no winds or dragons.

## Example

```jldoctest
julia> using MahjongTiles: wind, character, dragon, bamboo, circle, Hand, is_pure

julia> is_pure(Hand([wind(:east), character(3)]))
false

julia> is_pure(Hand([character(3), dragon(:green)]))
false

julia> is_pure(Hand([bamboo(8), character(3), circle(1)]))
true
```
"""
function is_pure(hand::Hand)
    return all(!is_suit(tile, Wind) && !is_suit(tile, Dragon) for tile in hand)
end

@doc """
    is_single_suit(hand)

Determine whether the hand consists of a single suit (一色),
which can include non-suits like winds or flowers.
To determine if the hand is only exactly a single suit,
excluding winds, dragons, and flowers/seasons,
combine this function with `is_pure`.

## Example

```jldoctest
julia> using MahjongTiles: wind, flower, character, circle, bamboo, Hand, is_single_suit

julia> is_single_suit(Hand([character(1), character(4), character(8), flower(3)]))
true

julia> is_single_suit(Hand([character(5), circle(5), bamboo(5)]))
false

julia> is_single_suit(Hand([circle(1), circle(7), wind(:north), wind(:east)]))
true
```
"""
function is_single_suit(hand::Hand)
    only_suit = suit(first(hand))
    return all(suit(tile) == only_suit || suit(tile) in (Wind, Dragon, Flower, Season) for tile in hand)
end

@doc """
    has_dragon_triple(hand)

Determine whether the hand has a triple of any one of the three dragon tiles.

## Example

```jldoctest
julia> using MahjongTiles: dragon, wind, bamboo, Hand, has_dragon_triple

julia> has_dragon_triple(Hand([dragon(:red), dragon(:red), dragon(:red), bamboo(2), wind(:east)]))
true

julia> has_dragon_triple(Hand([dragon(:red), dragon(:red), dragon(:green), dragon(:white)]))
false
```
"""
function has_dragon_triple(hand::Hand)
    counts = counter(hand.tiles)
    return counts[dragon(:red)] >= 3 || counts[dragon(:green)] >= 3 || counts[dragon(:white)] >= 3
end


text:
"""
  create_grid_stochastic(a::Array{Float64,1}, b::Array{Float64,1}, nums::Int64)
Function to create a grid. a is a vector of lower
bounds, b a vector of upper bounds and numPoints is the number of points
to be generated. The output is an Array{Float64,2}, where each row is a new point
and each column is a dimension of this points.
"""
function create_grid_stochastic(a::Array{Float64,1}, b::Array{Float64,1}, numPoints::Int64; gridType::Symbol = :uniform, alpha::Float64 = 0.025)

    #Each row is a new point and each column is a dimension of this points.
    #---------------------------------------------------------------------
    pointsFound = zeros(numPoints, length(a))

	  #Safety checks
	  #-------------
	  if numPoints < 1
	    Base.error("The input numPoints should be >= 1. numPoints = $(numPoints).")
	  end

	  if length(a) != length(b)
	    Base.error("length(a) != length(b)")
	  end

    if sum(a .< b) != length(b)
      Base.error("a .< b should hold")
    end

    # Loop until getting the desired number of points
    #-------------------------------------------------
    nbPointsFound = 0
    #Drawing from a normal
    #----------------------
    if gridType == :normal

      # Let's create a variance-covariance matrix such that the point
      # has 95% chances of being within the upper and upper bounds
      #---------------------------------------------------------------
      arrayStd = generate_std(a, b, numPoints, alpha)
      varCov = diagm(arrayStd.^2)

      d = MvNormal((a .+ b)./2, varCov)

      while nbPointsFound < numPoints
        draws = rand(d, numPoints)

        #Check whether point is within upper and lower bound
        #----------------------------------------------------
        for i=1:size(draws,2)
          # Take the sum to check all the inequalities at once
          if sum(draws[:,i] .>= a) == length(a) && sum(draws[:,i] .<= b) == length(a)
            nbPointsFound +=1
            # Check that the output is not "full"
            if nbPointsFound <= numPoints
              pointsFound[nbPointsFound,:] = draws[:,i]
            end

          end
        end

      end

    # Drawing from uniform
    #---------------------
    elseif gridType == :uniform

      # uniform along each dimension
      #-----------------------------
      # Loop over columns
      for j=1:size(pointsFound,2)
        # Loop over rows
        for i=1:size(pointsFound,1)

          # Remember: each row is a new point and each column is a dimension
          # of this points.
          #------------------------------------------------------------------
          d = Uniform(a[j], b[j])

          pointsFound[i,j] = rand(d, 1)[1]
          nbPointsFound +=1

        end
      end

    else
      Base.error("gridType has to be :normal or :uniform")
    end

    return pointsFound

end

"""
  generate_std(a::Array{Float64,1}, b::Array{Float64,1}, numPoints::Int64, alpha::Float64)
Function to generate standard errors such that there is (1-alpha)% chances for
tge point to fall within the upper and lower bound, when sampling from the normal distribution
"""
function generate_std(a::Array{Float64,1}, b::Array{Float64,1}, numPoints::Int64, alpha::Float64)

  # Output:
  #--------
  arrayStd = zeros(length(b))

  #Safety checks
  #-------------
  if numPoints < 1
    Base.error("The input numPoints should be >= 1. numPoints = $(numPoints).")
  end

  if length(a) != length(b)
    Base.error("length(a) != length(b)")
  end

  if alpha <= 0 && alpha > 1
    Base.error("alpha should be such that 0 < alpha < 1")
  end

  if sum(a .< b) != length(b)
    Base.error("a .< b should hold")
  end

  meanValue = (a .+ b)./2

  for i=1:length(b)

      arrayStd[i] = (sqrt(numPoints)*(b[i]-meanValue[i]))/quantile(Normal(), 1 - (alpha/2))

  end

  return arrayStd

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

include("testdefs.jl")
addprocs(4; topology="master_slave")
using Base.Test

@test_throws RemoteException remotecall_fetch(()->remotecall_fetch(myid, 3), 2)

function test_worker_counts()
    # check if the nprocs/nworkers/workers are the same on the remaining workers
    np=nprocs()
    nw=nworkers()
    ws=sort(workers())

    for p in workers()
        @test (true, true, true) == remotecall_fetch(p, np, nw, ws) do x,y,z
            (x==nprocs(), y==nworkers(), z==sort(workers()))
        end
    end
end

function remove_workers_and_test()
    while nworkers() > 0
        rmprocs(workers()[1]; waitfor=2.0)
        test_worker_counts()
        if nworkers() == nprocs()
            break
        end
    end
end

remove_workers_and_test()

# connect even pids to other even pids, odd to odd.
mutable struct TopoTestManager <: ClusterManager
    np::Integer
end

function Base.launch(manager::TopoTestManager, params::Dict, launched::Array, c::Condition)
    dir = params[:dir]
    exename = params[:exename]
    exeflags = params[:exeflags]

    for i in 1:manager.np
        io, pobj = open(pipeline(detach(
            setenv(`$(Base.julia_cmd(exename)) $exeflags --bind-to $(Base.LPROC.bind_addr) --worker $(Base.cluster_cookie())`, dir=dir)); stderr=STDERR), "r")
        wconfig = WorkerConfig()
        wconfig.process = pobj
        wconfig.io = io
        wconfig.ident = i
        wconfig.connect_idents = collect(i+2:2:manager.np)
        push!(launched, wconfig)
    end

    notify(c)
end

const map_pid_ident=Dict()
function Base.manage(manager::TopoTestManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :register
        map_pid_ident[id] = get(config.ident)
    elseif op == :interrupt
        kill(get(config.process), 2)
    end
end

addprocs(TopoTestManager(8); topology="custom")

while true
    if any(x->get(map_pid_ident, x, 0)==0, workers())
        yield()
    else
        break
    end
end

for p1 in workers()
    for p2 in workers()
        i1 = map_pid_ident[p1]
        i2 = map_pid_ident[p2]
        if (iseven(i1) && iseven(i2)) || (isodd(i1) && isodd(i2))
            @test p2 == remotecall_fetch(p->remotecall_fetch(myid, p), p1, p2)
        else
            @test_throws RemoteException remotecall_fetch(p->remotecall_fetch(myid, p), p1, p2)
        end
    end
end

remove_workers_and_test()

max_stars_repo_name: felixrehren/julia

max_issues_repo_name: felixrehren/julia

max_forks_repo_name: felixrehren/julia