text:
export @synapse_str
macro synapse_str(synapse)
  delayMax=maximum(net.delay)
  delayPeriod=ceil(Int,delayMax/net.dt)
  if delay!=0
    strInit=""
    str=""
    for line in split(synapse,"\n")
      var1=match(r"^.*(?=\+)",line).match
      var2=match(r"(?<==).*",line).match

      strInit*="net.$(var1)Temp=zeros(net.N,$delayPeriod)\n";
      str*="net.$(var1)+=net.$(var1)Temp[:,mod1(net.t/net.dt,$delayPeriod)]\n"
      str*="forwardIndex=mod1(round(Int,(net.t+net.delay[post,pre])/net.dt),$delayPeriod)\n"
      str*="net.$(var1)Temp[post,forwardIndex]+=net.$(var2)[post,pre]\n"
    end
    str="for pre = net.i\nfor post=1:N\n"*str*"end\nend\n"
  end
  if delay==0
    str=""
    strInit=""
    for line in split(synapse,"\n")
      var1=match(r"^.*(?=\+)",line).match
      var2=match(r"(?<==).*",line).match
      str*="net.$(var1)+=sum(net.$(var2)[:,net.i],2)\n"
    end
  end

  print(strInit)
  map(x->eval(parse(x)),split(strInit,"\n"))

  str="function synapse!(net::Network$timestr)\n"*str
  str*="end\n"
  print(str)
  eval(parse(str))
end


text:
using BooSTjl
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
function test()
    x::Matrix{Float64} = [randn(10) randn(10)]
    y::Vector{Float64} = randn(10)

    t1 = BooST(x,y, M = 5,display = false)
    t2 = predictBooST(t1,x)
    t3 = BooSTMore(x,y,t1,M = 5, display = false)
    t4 = estimate_derivatives(t1,x,1)
    return t1, t2, t3, t4

end

(t1,t2,t3,t4) = test()

@test typeof(t1) == Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64,Array{Float64,1}}
@test typeof(t2) == Vector{Float64}
@test typeof(t3) == Tuple{Dict{Float64,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Float64,Tuple{Float64,Float64,Int64,Array{Float64,1},Int64,Float64},Int64,Array{Float64,1}}
@test typeof(t4) == Vector{Float64}

ext: jl

max_stars_repo_name: gabrielrvsc/BooSTjl

max_issues_repo_name: gabrielrvsc/BooSTjl

max_forks_repo_name: gabrielrvsc/BooSTjl


text:
using DocumenterLaTeX, Markdown
import Documenter: Documents, Documenter, Writers, Utilities
import Documenter.Writers.LaTeXWriter: piperun, _print

const LaTeX_CC="xelatex"
const DOCKER_IMAGE = "tianjun2018/documenter-latex:latest"

function Documenter.Writers.LaTeXWriter.latexinline(io, math::Markdown.LaTeX)
    # Handle MathJax and TeX inconsistency since the first wants `\LaTeX` wrapped
    # in math delims, whereas actual TeX fails when that is done.
    math.formula == "\\LaTeX" ? _print(io, " ", math.formula, " ") : _print(io, " \\(", math.formula, "\\) ")
end

function Documenter.Writers.LaTeXWriter.compile_tex(doc::Documents.Document, settings::LaTeX, texfile::String)
    if settings.platform == "native"
        Sys.which("latexmk") === nothing && (@error "LaTeXWriter: latexmk command not found."; return false)
        @info "LaTeXWriter: using latexmk to compile tex."
        piperun(`latexmk -f -interaction=nonstopmode -view=none -$(LaTeX_CC) -shell-escape $texfile`)
        return true

        # NOTE: 中文排版依然有问题，我们暂时不管他们。输出的pdf大部分情况下依然可以使用。
        # try
        #     piperun(`latexmk -f -interaction=nonstopmode -view=none -$(LaTeX_CC) -shell-escape $texfile`)
        #     return true
        # catch err
        #     logs = cp(pwd(), mktempdir(); force=true)
        #     @error "LaTeXWriter: failed to compile tex with latexmk. " *
        #            "Logs and partial output can be found in $(Utilities.locrepr(logs))." exception = err
        #     return false
        # end
    elseif settings.platform == "docker"
        Sys.which("docker") === nothing && (@error "LaTeXWriter: docker command not found."; return false)
        @info "LaTeXWriter: using docker to compile tex."
        script = """
            mkdir /home/zeptodoctor/build
            cd /home/zeptodoctor/build
            cp -r /mnt/. .
            latexmk -f -interaction=nonstopmode -view=none -$(LaTeX_CC) -shell-escape $texfile
            """
        try
            piperun(`docker run -itd -u zeptodoctor --name latex-container -v $(pwd()):/mnt/ --rm $(DOCKER_IMAGE)`)
            piperun(`docker exec -u zeptodoctor latex-container bash -c $(script)`)
            piperun(`docker cp latex-container:/home/zeptodoctor/build/. .`)
            return true
        catch err
            logs = cp(pwd(), mktempdir(); force=true)
            @error "LaTeXWriter: failed to compile tex with docker. " *
                   "Logs and partial output can be found in $(Utilities.locrepr(logs))." exception = err
            return false
        finally
            try; piperun(`docker stop latex-container`); catch; end
        end
    end
end


text:
## device interface

# hybrid Tausworthe and Linear Congruent generator from
# https://developer.nvidia.com/gpugems/GPUGems3/gpugems3_ch37.html

function TausStep(z::Unsigned, S1::Integer, S2::Integer, S3::Integer, M::Unsigned)
    b = (((z << S1) ⊻ z) >> S2)
    return (((z & M) << S3) ⊻ b)
end

LCGStep(z::Unsigned, A::Unsigned, C::Unsigned) = A * z + C

make_rand_num(::Type{Float64}, tmp) = 2.3283064365387e-10 * Float64(tmp)
make_rand_num(::Type{Float32}, tmp) = 2.3283064f-10 * Float32(tmp)

function next_rand(::Type{FT}, state::NTuple{4, T}) where {FT, T <: Unsigned}
    state = (
        TausStep(state[1], Cint(13), Cint(19), Cint(12), T(4294967294)),
        TausStep(state[2], Cint(2), Cint(25), Cint(4), T(4294967288)),
        TausStep(state[3], Cint(3), Cint(11), Cint(17), T(4294967280)),
        LCGStep(state[4], T(1664525), T(1013904223))
    )
    tmp = (state[1] ⊻ state[2] ⊻ state[3] ⊻ state[4])
    return (
        state,
        make_rand_num(FT, tmp)
    )
end

function gpu_rand(::Type{T}, state, randstate::AbstractVector{NTuple{4, UInt32}}) where T
    threadid = GPUArrays.threadidx_x(state)
    stateful_rand = next_rand(T, randstate[threadid])
    randstate[threadid] = stateful_rand[1]
    return stateful_rand[2]
end

floattype(::Type{T}) where T <: Union{Int64, UInt64} = Float64
floattype(::Type{T}) where T <: Union{Int32, UInt32} = Float32

to_number_range(x::AbstractFloat, ::Type{T}) where T <: Unsigned = T(round(x * typemax(T)))

to_number_range(x::F, ::Type{T}) where {T <: Signed, F <: AbstractFloat} = Base.unsafe_trunc(T, round(((x - F(0.5)) * typemax(T)) * T(2)))

function gpu_rand(::Type{T}, state, randstate::AbstractVector{NTuple{4, UInt32}}) where T <: Integer
    f = gpu_rand(floattype(T), state, randstate)
    return to_number_range(f, T)
end

## host interface

struct RNG <: AbstractRNG
    state::GPUArray{NTuple{4,UInt32},1}
end

const GLOBAL_RNGS = Dict()
function global_rng(A::GPUArray)
    dev = GPUArrays.device(A)
    get!(GLOBAL_RNGS, dev) do
        N = GPUArrays.threads(dev)
        state = similar(A, NTuple{4, UInt32}, N)
        copyto!(state, [ntuple(i-> rand(UInt32), 4) for i=1:N])
        RNG(state)
    end
end

function Random.rand!(rng::RNG, A::GPUArray{T}) where T <: Number
    gpu_call(A, (rng.state, A,)) do state, randstates, a
        idx = linear_index(state)
        idx > length(a) && return
        @inbounds a[idx] = gpu_rand(T, state, randstates)
        return
    end
    A
end

Random.rand!(A::GPUArray) = rand!(global_rng(A), A)


text:
using ModelObjectsLH, ModelParams, Test

mdl = ModelParams;

function pv_collection_test()
    @testset "Main" begin
        n = 5;
        pvv = mdl.make_test_pvector_collection(n);

        @test !isempty(pvv);
        @test length(pvv) == n;

        for (objId, pvec) in pvv
            @test objId == get_object_id(pvec);
            @test has_pvector(pvv, objId);
            pv = find_pvector(pvv, objId);
            @test isequal(get_object_id(pv), objId);
        end
    end
end

function constructors_test()
    @testset "Constructors" begin
        n = 3;
        pv = [mdl.make_test_pvector(2+j; objId = Symbol("obj$j"))  for j = 1 : n];
        pvecV = PVectorCollection(pv);
        @test length(pvecV) == n;
    end
end

function make_dict_test()
    @testset "Make Dict" begin
        n = 6;
        pvv = mdl.make_test_pvector_collection(n);
        d = mdl.make_dict(pvv; isCalibrated = true);
        @test length(d) == n
        for objIdStr in keys(d)
            objId = make_object_id(objIdStr);
            pvec = find_pvector(pvv, objId);
            @test isequal(d[objIdStr],  make_dict(pvec; isCalibrated = true, valueType = :value))
        end
    end
end

function find_param_test()
    @testset "Find param" begin
        n = 5;
        pvv = mdl.make_test_pvector_collection(n);
        d = find_param(pvv, :notThere);
        @test isempty(d);
        pName = :p4;
        d = find_param(pvv, pName);
        for (objId, p) in d
            pvec = find_pvector(pvv, objId);
            @test param_exists(pvec, pName);
        end
    end
end

@testset "PVectorCollection" begin
    constructors_test();
    pv_collection_test();
    make_dict_test();
    find_param_test();
end

# -----------


text:
#=
Absolute:
- Julia version: 1.4.2
- Author: akhil
- Date: 2020-05-28
=#
"""
Evaluate Absolute Value of a number

# Examples
```jldoctest
julia> Maths.abs_value(5)
5
```
"""
function abs_value(x)
    return x>0 ? x : -x
end


text:
using Distributions
using ReverseDiffOverload
import Distributions.logpdf

test(dist, x, params...) = begin
    @assert isapprox(logpdf(dist, x, params...), logpdf(dist(params...), x))
    if dist <: DiscreteDistribution
        testdiff((params...) -> logpdf(dist, x, params...), params...)
    else
        testdiff((x, params...) -> logpdf(dist, x, params...), x, params...)
    end
end

#Normal
valid(::Type{Normal}, x, μ, σ) = begin
    @assert all(σ .> 0)
end
#Differentiable versions.
logpdf(::Type{Normal}, x, μ, σ) = begin
    #valid(Normal, x, μ, σ)
    σ2 = σ.*σ
    diff = x-μ
    -(log(2*pi*σ2) + diff.*diff ./ σ2)/2
end
test(Normal, 1, 2, 3)

#MvNormal
valid(::Type{MvNormal}, x::Vector, μ::Vector, Σ::Matrix) = begin
    @assert isposdef(Σ)
end
#Differentiable versions.
logpdf(::Type{MvNormal}, x::Vector, μ::Vector, Σ::Matrix) = begin
    #valid(MvNormal, x, μ, Σ)
    diff = x-μ
    -(length(x)*log(2*pi) + logdet(Σ) + (diff'/Σ)*x)/2
end
#test(MvNormal, [1., 2], [3., 4], eye(2) + ones(2,2))

#Gamma
valid(::Type{Gamma}, x, shape, scale) = begin
    @assert all(x .>= 0)
    @assert all(shape .>= 0)
    @assert all(scale .>= 0)
end
logpdf(::Type{Gamma}, x, shape, scale) = begin
    #valid(Gamma, x, shape, scale)
    -logΓ(shape) - shape.*log(scale) + (shape-1).*log(x) - x./scale
end
test(Gamma, 1, 2, 3)

#InverseGamma
valid(::Type{InverseGamma}, x, shape, scale) = begin
    @assert all(x .>= 0)
    @assert all(shape .>= 0)
    @assert all(scale .>= 0)
end
logpdf(::Type{InverseGamma}, x, shape, scale) = begin
    #valid(InverseGamma, x, shape, scale)
    -logΓ(shape) + shape.*log(scale) - (shape+1).*log(x) - scale./x
end
test(InverseGamma, 1, 2, 3)

#Bernoulli
valid(::Type{Bernoulli}, x, p) = begin
    @assert all(0 .<= p .<= 1)
    @assert all((x .== 1)|(x .== 0))
end
logpdf(::Type{Bernoulli}, x, p) = begin
    #valid(Bernoulli, x, p)
    x*log(p) + (1-x)*log(1-p)
end
test(Bernoulli, 1., 0.9)
test(Bernoulli, 0., 0.9)

#Categorical
valid(::Type{Categorical}, x, ps::Vector{Float64}) = begin
    @assert isapprox(sum(ps), 1)
    @assert all(0 .<= ps .<= 1)
end
logpdf(::Type{Categorical}, x, ps) = begin
    #valid(Categorical, x, ps)
    log(ps[x])
end
test(Categorical, 2, [0.3, 0.4, 0.3])

#Beta
valid(::Type{Beta}, x, α, β) = begin
    @assert all(α .> 0)
    @assert all(β .> 0)
    @assert all(0 .<= x .<= 1)
end
logB(α, β) = logΓ(α) + logΓ(β) - logΓ(α + β)
logpdf(::Type{Beta}, x, α, β) = begin
    #valid(Beta, x, α, β)
    (α-1).*log(x) + (β-1).*log(1-x) - logB(α, β)
end
test(Beta, 0.6, 1.3, 1.5)

#Dirichlet
valid(::Type{Dirichlet}, x, α::Vector{Float64}) = begin
    @assert all(0 .<= x)
    @assert all(0 .<= α)
end
logB(α) = sum(logΓ(α)) - logΓ(sum(α))
logpdf(::Type{Dirichlet}, x, α) = begin
    #valid(Dirichlet, x, α)
    sum((α-1).*log(x)) - logB(α)
end
test(Dirichlet, [1.2, 1.6, 1.3], [1.3, 1.6, 1.3])

#Poisson
valid(::Type{Poisson}, k, λ) = begin
    @assert all(k .>= 0)
    @assert all(λ .>= 0)
end
logpdf(::Type{Poisson}, k, λ) = begin
    #valid(Poisson, k, λ)
    k*log(λ) - log(factorial(k)) - λ
end
test(Poisson, 3, 5)

#Non-differentiable fallbacks from Distributions.jl
logpdf{D <: Distribution}(d::Type{D}, x, params...) =
    logpdf(d(params...), x)


text:
using Test, MathieuFunctions
using LinearAlgebra
using DelimitedFiles

readcsv(f) = DelimitedFiles.readdlm(f, ',')

function tapprox(a, b; atol=1e-15)
    normval = norm(a - b, Inf)
    @info "normval = $normval"
    isapprox(a, b; norm= x -> norm(x, Inf), atol = atol)
end

@testset "basic" begin
    @test charλ(0, 1, k=1:10) == [1.0, 1.0, 9.0, 9.0, 25.0, 25.0, 49.0, 49.0, 81.0, 81.0]
    @test maximum(charA(0,k=0:100) - [0:100;].^2) == 0
    @test norm(charB(0,k=1:100) - [1:100;].^2) == 0
end

filename = "MathieuCharacteristicA-1.csv"
@testset "$filename" begin
    test1 = readcsv(filename)
    r = reduce(hcat,Vector{Float64}[charA(q,k=0:10) for q in [-10:.01:10;]])
    @test tapprox(test1, r;  atol=7.5e-13)
end

filename = "MathieuCharacteristicA-2.csv"
@testset "$filename" begin
    test1 = readcsv(filename)
    r = reduce(hcat,Vector{Float64}[charA(q,k=0:3) for q in [30:.01:50;]])
    @test tapprox(test1, r;  atol=7.6e-13) # NOTE: was 7.5e-13
end

filename = "MathieuCharacteristicB-1.csv"
@testset "$filename" begin
    test1 = readcsv(filename)
    r = reduce(hcat,Vector{Float64}[charB(q,k=1:10) for q in [-10:.01:10;]])
    @test tapprox(test1, r;  atol=7.5e-13)
end

filename = "MathieuCharacteristicB-2.csv"
@testset "$filename" begin
    test1 = readcsv(filename)
    r = reduce(hcat,Vector{Float64}[charB(q,k=1:3) for q in [30:.01:50;]])
    @test tapprox(test1, r; atol=2.8e-11)
end

filename = "MathieuCharacteristicL-1.csv"
@testset "$filename" begin
    test1 = readcsv(filename)[1:100,:]
    test2 = Float64[charλ(q,ν,k=1:1)[1] for ν in [0:.01:0.99;], q in [-5:.01:5;]]
    @test tapprox(test1, test2, atol=7.5e-15)
    # TODO: test ν > 1 (currently failing)
end

filename = "MathieuCharacteristicL-2.csv"
@testset "$filename" begin
    test1 = readcsv(filename)[1:100,:]
    test2 = Float64[charλ(q,ν,k=1:1)[1] for ν in [0:.01:0.99;], q in [30:.01:50;]]
    @test tapprox(test1, test2, atol=4.5e-14)
    # TODO: test ν > 1 (currently failing)
end


text:
#------------------------------------------------------------------------------

# GroebnerMetainfo
# Stores parameters for groebner algorithm
struct GroebnerMetainfo{Rng}
    # if set, then use fglm algorithm for order conversion
    usefglm::Bool
    # output polynomials order
    targetord::Symbol
    # order for computation
    computeord::Symbol

    # correctness checks levels
    heuristiccheck::Bool
    randomizedcheck::Bool
    guaranteedcheck::Bool

    # linear algebra backend to be used
    linalg::Symbol

    rng::Rng
end

function set_metaparameters(ring, ordering, certify, forsolve, linalg, rng)
    usefglm = false
    targetord = :lex
    computeord = :lex

    if forsolve
        targetord = :lex
        usefglm = true
        if ordering in (:deglex, :degrevlex)
            computeord = ordering
        else
            computeord = :degrevlex
        end
        # TODO: lex fglm
        computeord = :lex
    else
        if ordering == :input
            ordering = ring.ord
        end
        targetord = ordering
        usefglm = false
        computeord = targetord
    end

    heuristiccheck = true
    # heuristiccheck = false
    randomizedcheck = true
    if certify
        guaranteedcheck = true
    else
        guaranteedcheck = false
    end

    @info "Computing in $computeord order, result is in $targetord order"
    @info "Using fglm: $usefglm"

    GroebnerMetainfo(usefglm, targetord, computeord,
                        heuristiccheck, randomizedcheck, guaranteedcheck,
                        linalg, rng)
end

function select_tablesize(ring, exps)
    nvars = ring.nvars
    sz = length(exps)

    tablesize = 2^10
    if nvars > 4
        tablesize = 2^14
    end
    if nvars > 7
        tablesize = 2^16
    end

    if sz < 3
        tablesize = div(tablesize, 2)
    end
    if sz < 2
        tablesize = div(tablesize, 2)
    end

    tablesize
end


text:
"""
    StreamConfig

Holds the conguration for livestream, defaults to `nothing`

#Fields
- `livestreamto::Symbol` Livestream platform `:local` or `:twitch`
- `protocol::String` The streaming protocol to be used. Defaults to UDP
- `address::String` The IP address for the `:local` stream(ignored in case of `:twitch`)
- `port::Int` The port for the `:local` stream(ignored in case of `:twitch`)
- `twitch_key::String` Twitch stream key for your account
"""
struct StreamConfig
    livestreamto::Symbol
    protocol::String
    address::String
    port::Int
    twitch_key::String
end


text:
# Makes a vector in the standard basis

export eVec

""" `vec = eVec(d, j)`

Returns a sparse vector with all 0's except at the position j.
"""

function eVec(dim::Int, j::Int)
	@assert dim >= 0 "Dimension must be positive."
	@assert j >= 0 "Index must be positive."
	@assert j <= dim "Index must be less than dimension."
	return sparsevec([j], [1], dim)
end


text:
module Verify

export create_nb_dict, get_nb

using ..Common

immutable Link{G}
    generators::Vector{G}
    even::Bool

    function Link(generators::Vector, even::Bool, k::Int)
        new(except_index(generators, k), even $ isodd(length(generators)-k))
    end
end
Base.isequal(l1::Link, l2::Link) = (l1.generators == l2.generators && l1.even == l2.even)
Base.hash(link::Link) = hash(link.generators) $ link.even

function add_links!{F<:Facet,L<:Link}(links::Dict{L,F}, facets)
    for facet in facets
        generators, even = get_canonical_winding(facet)
        for k in 1:nfacet(facet)
            link = L(generators, even, k)
            if haskey(links, link); error("Link already exists!"); end
            links[link] = facet
        end
    end
end

function create_nb_dict{F<:Facet}(::Type{F}, facets)
    links = Dict{Link{gtype(F)}, F}()
    add_links!(links, facets)
    links
end

function get_nb{F<:Facet,L<:Link}(links::Dict{L,F}, facet::F, g)
    generators, even = get_canonical_winding(facet)
    k = indexof(generators, g)
    link = L(generators, !even, k)
    links[link]
end

end # module


text:
type JamshidianSwaptionEngine{S <: ShortRateModel} <: PricingEngine
  model::S
  # ts::Y

  # function call{S}(::Type{JamshidianSwaptionEngine}, model::S)
  #   new{S, YieldTermStructure}(model)
  # end
  #
  # function call{S, Y}(::Type{JamshidianSwaptionEngine}, model::S, yts::Y)
  #   new{S, Y}(model, yts)
  # end
end

# methods #
function _calculate!(pe::JamshidianSwaptionEngine, swaption::Swaption)
  tsmodel = pe.model
  ref_date = reference_date(tsmodel.ts)
  dc = tsmodel.ts.dc

  amounts = copy(swaption.swap.args.fixedCoupons)
  amounts[end] += swaption.swap.nominal

  maturity = year_fraction(dc, ref_date, swaption.exercise.dates[1])

  fixedPayTimes = zeros(length(swaption.swap.args.fixedPayDates))
  valueTime = year_fraction(dc, ref_date, swaption.swap.args.fixedResetDates[1])
  for i = 1:length(fixedPayTimes)
    fixedPayTimes[i] = year_fraction(dc, ref_date, swaption.swap.args.fixedPayDates[i])
  end

  finder = RStarFinder(tsmodel, swaption.swap.nominal, maturity, valueTime, fixedPayTimes, amounts)

  minStrike = -10.0
  maxStrike = 10.0
  slv = BrentSolver(10000, true, true, minStrike, maxStrike)

  rStar = solve(slv, operator(finder), 1e-8, 0.05, minStrike, maxStrike)

  w = isa(swaption.swap.swapT, Payer) ? Put() : Call()

  _size = length(swaption.swap.args.fixedCoupons)

  val = 0.0
  _B = discount_bond(tsmodel, maturity, valueTime, rStar)

  @simd for i = 1:_size
    @inbounds fixedPayTime = year_fraction(dc, ref_date, swaption.swap.args.fixedPayDates[i])

    strike = discount_bond(tsmodel, maturity, fixedPayTime, rStar) / _B

    dboValue = discount_bond_option(tsmodel, w, strike, maturity, valueTime, fixedPayTime)

    @inbounds val += amounts[i] * dboValue
  end

  swaption.results.value = val

  return swaption
end


text:
using Revise
using AlphaZero
using ProgressLogging
using TensorBoardLogger
using Logging
using ParameterSchedulers: Scheduler, Cos
import Flux
import FileIO
import Random
import CUDA

import LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

include("../mu_game_wrapper.jl")
include("../network.jl")
# include("../alphazerolike.jl")
include("../trace.jl")
include("../play.jl")
include("../simulations.jl")
include("../training.jl")
include("../learning.jl")
include("../benchmark.jl")
include("../probe-games.jl")

device = Flux.cpu
# device = Flux.gpu
# CUDA.allowscalar(true)

gspec = Examples.games["connect-four"]
# gspec = ProbeGames.games["a1r1"]
# gspec = ProbeGames.games["a1o2r2"]
# gspec = ProbeGames.games["a1r1s2"] # representation collapses states into the same one
# gspec = ProbeGames.games["a2r2"]
# gspec = ProbeGames.games["twoplayer"]
# gspec = ProbeGames.games["simplertictactoe"]

n=32
self_play = (;
  sim=SimParams(
    num_games=100,
    num_workers=n,
    # TODO make batch sizes for each network
    batch_size=16,
    use_gpu=false,
    # todevice=Flux.cpu, #?
    reset_every=4, #not used, mcts resets everytime
    flip_probability=0.,
    alternate_colors=false),
  mcts = MctsParams(
    num_iters_per_turn=128, #1000 benchmark
    cpuct=1.25,
    temperature=ConstSchedule(1.0),
    dirichlet_noise_ϵ=0.25,
    dirichlet_noise_α=0.1),
  device=Flux.gpu)

# # for mcts rollout generation
# self_play = (;
#   sim=SimParams(
#     num_games=10_000,
#     num_workers=n,
#     batch_size=n,
#     use_gpu=false,
#     reset_every=4, #not used, mcts resets everytime
#     flip_probability=0.,
#     alternate_colors=false),
#   mcts = MctsParams(
#     num_iters_per_turn=500, #1000 benchmark
#     cpuct=2.5,
#     temperature=ConstSchedule(1.0),
#     dirichlet_noise_ϵ=0.25,
#     dirichlet_noise_α=0.1))

arena = (;
  sim=SimParams(
    num_games=50,
    num_workers=16,
    batch_size=16,
    use_gpu=false,
    reset_every=1,
    flip_probability=0., #0.5
    alternate_colors=true),
  mcts = MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.3),
    dirichlet_noise_ϵ=0.1),
  update_threshold=0.0,
  device=Flux.cpu)

learning_params = (;
  num_unroll_steps=10, #if =0, g is not learning, muzero-general=20
  td_steps=50, # with max length=9, always go till the end of the game, rootvalues don't count
  discount=0.997,
  #// value_loss_weight = 0.25, #TODO
  l2_regularization=1f-4, #Float32
  #// l2_regularization=0f0, #Float32
  loss_computation_batch_size=512,
  batches_per_checkpoint=10,
  num_checkpoints=1,
  opt=Scheduler(
    Cos(λ0=1e-3, λ1=1e-5, period=10000), # cosine annealing, google 2e4, generat doesn't use any
    Flux.ADAM()
  ),
  device=Flux.gpu
)

benchmark_sim = SimParams(
    num_games=400,
    num_workers=64,
    batch_size=32,
    use_gpu=false,
    reset_every=1,
    flip_probability=0.0, #0.5
    alternate_colors=true)

bench_mcts = MctsParams(
  num_iters_per_turn=500, #1000 benchmark
  cpuct=1.25,
  temperature=ConstSchedule(0.3),
  dirichlet_noise_ϵ=0.25,
  dirichlet_noise_α=0.1)

benchmark = (;
  vanilla_mcts=Benchmark.Duel(
    Mu(arena.mcts),
    Benchmark.MctsRollouts(bench_mcts),
    benchmark_sim),
  # minmax_d5=Benchmark.Duel(
  #   Mu(arena.mcts),
  #   Benchmark.MinMaxTS(depth=5, amplify_rewards=true, τ=1.),
  #   benchmark_sim)
)

  #memory = CircularBuffer{Trace{GI.state_type(gspec)}}(1024)

  # μparams = MuParams(self_play, learning_params, arena, 6, 3000)
  μparams = MuParams(self_play, learning_params, arena, 10_000, 10_000)
  # μparams = MuParams(self_play, learning_params, arena, 2, 3000)

  # αnns = (
  #   f=PredictionNetwork(
  #     gspec,
  #     PredictionHP(hiddenstate_shape=27, width=256, depth_common=6,
  #       use_batch_norm=true, batch_norm_momentum=1.)),
  #   g=AlphaDynamics(gspec), # set num_unroll_steps=0 to not use it during learning
  #   h=AlphaRepresentation())

  # αnns = (
  #   f=AlphaPrediction(
  #       PredictionNetwork(
  #         gspec,
  #         PredictionHP(hiddenstate_shape=27, width=256, depth_common=6,
  #           use_batch_norm=true, batch_norm_momentum=1.))),
  #   g=AlphaDynamics(gspec), # set num_unroll_steps=0 to not use it during learning
  #   h=AlphaRepresentation())
  Random.seed!(2137)
  # env = MuEnv(gspec, μparams, αnns, experience=FileIO.load("memory_minmax.jld2", "mem"))
  # env = MuEnv(gspec, μparams, αnns, experience=FileIO.load("memory_minmax_d5_t025.jld2", "mem"))

# muzero-general params
μNetworkHP = MuNetworkHP(
  gspec,
  PredictionHP(hiddenstate_shape=64, width=64, depth_common=0,
    depth_vectorhead=0, depth_scalarhead=0, use_batch_norm=false, batch_norm_momentum=1.),
  DynamicsHP(hiddenstate_shape=64, width=64, depth_common=-1,
   depth_vectorhead=0, depth_scalarhead=0,  use_batch_norm=false, batch_norm_momentum=1.),
  RepresentationHP(hiddenstate_shape=64, width=0, depth=-1))

  env = MuEnv(gspec, μparams, MuNetwork(μNetworkHP))
  # env = MuEnv(gspec, μparams, deepcopy(env.bestnns))
  # env = MuEnv(gspec, μparams, MuNetwork(μNetworkHP), experience=FileIO.load("results/c4_5x4/memory_mctsrollout500.jld2", "mem"))
  # env = MuEnv(gspec, μparams, μNetwork, experience=FileIO.load("results/c4_5x4/memory_mctsrollout500.jld2", "mem"))

##

path = "results/c4_5x4/" * format(now(),"yyyy-mm-ddTHHMM") * "_custominit33_testbatching/"
tblogger=TBLogger(path, min_level=Logging.Info)
  # TRAIN
  with_logger(tblogger) do
    @info "params" params=μparams
    # @info "alphalike params" αnns.f.net.hyper
    # @timed learning_step!(env)
    @info "network params" μNetworkHP
    # @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player"
    @info "Benchmark" run_duel(env, benchmark)...
    mutrain!(env, benchmark=benchmark, path=path)
  end

  # SUPERVISED LEARNING
  with_logger(tblogger) do
    @info "params" params=μparams
    @info "network params" μNetworkHP
    # @info "alphalike params" αnns.f.net.hyper
    # train!(env; benchmark=benchmark)
    @info "Memory Analysis" memory_analysis(env.memory)... "generated by MCTS player"
    @info "Benchmark" run_duel(env, benchmark)...
    function supervised_learning!(env, benchmark, num_iterations, path)
      while env.itc <= num_iterations
        _, time_learning = @timed learning_step!(env)
        @info "Training" stage="iteration finished" env.itc time_learning
        if env.itc % 50 == 0
          @info "Benchmark" run_duel(env, benchmark)...
          FileIO.save(path*"$(format(now(),"yyyy-mm-ddTHHMM"))_env_$(env.itc).jld2", "env", env)
        end
        env.itc += 1
      end
    end
    supervised_learning!(env, benchmark, 1500, path)
  end

  with_logger(tblogger) do
  self_play_step!(env)
end

  @enter self_play_step!(env)
  @profview self_play_step!(env)
  @timed self_play_step!(env)
  @timed learning_step!(env)
  @profview run_duel(env, benchmark)
  mutrain!(env)

  # TRAIN without benchmark
  with_logger(tblogger) do
    @info "params" params=μparams
    # @info "alphalike params" αnns.f.net.hyper
    # @timed learning_step!(env)
    @info "network params" μNetworkHP
    # @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player"
    mutrain!(env)
  end

  # TRAIN with pretrain
  with_logger(tblogger) do
    @info "params" params=μparams
    # @info "alphalike params" αnns.f.net.hyper
    @info "network params" μNetworkHP
    # @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player"
    @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player"
    _ = run_duel(env, benchmark)
    for _ in 1:50
      learning_step!(env)
    end
    _ = run_duel(env, benchmark)
    mutrain!(env, benchmark=benchmark)
  end

benchmark = [
  Benchmark.Duel(
    Mu(arena.mcts),
    Benchmark.MctsRollouts(self_play.mcts),
    benchmark_sim)]
benchmark = [
  Benchmark.Duel(
    Mu(arena.mcts),
    Benchmark.MctsRollouts(arena.mcts),
    benchmark_sim)]

run_duel(env,benchmark)
AlphaZero.Util.@printing_errors train!(env, benchmark=benchmark)

self_play_step!(env)
  @enter rewards, redundancy = pit_networks(gspec, env.curnns, env.bestnns, arena)
  rewards, redundancy = pit_networks(gspec, env.curnns, env.bestnns, arena)
  mean(rewards)
# FileIO.save("env_longminmaxlearning.jld2", "env", env)
# FileIO.save("env_muzerogeneral.jld2", "env", env)
# FileIO.save("env_muzerogeneral_minmaxpretrained.jld2", "env", env)
  train!(env)

#generate minmax samples
  # with_logger(tblogger) do
  env = MuEnv(gspec, μparams, αnns)

    @info "params" params=μparams
    AlphaZero.Util.@printing_errors @timed self_play_step!(env)
    @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player τ=0.25, depth=5"
    # FileIO.save("memory_minmax_d5_t025.jld2", "mem", env.memory)
    # FileIO.save("memory_minmax_randominits_d5_t025.jld2", "mem", env.memory)
    # FileIO.save("memory_minmax_d7.jld2", "mem", env.memory)
    # FileIO.save("memory_mcts_rollout_300.jld2", "mem", env.memory)
    # FileIO.save("memory_mctsrollout300_custominit1.jld2", "mem", env.memory)
    # FileIO.save("memory_mctsrollout64.jld2", "mem", env.memory)
    FileIO.save("results/c4_5x4/memory_mctsrollout500.jld2", "mem", env.memory)
  # end

count(last(t.rewards)==(1) for t in env.memory)
env = envsim
  hyper = env.params.learning_params
  traces = [sample_trace(env.memory) for _ in 1:2] #? change to iterator
  trace_pos_idxs = [sample_position(t) for t in traces]
  sample = make_target(gspec, traces[1], 2, hyper)
  samples = [make_target(gspec, t, i, hyper) for (t,i) in zip(traces, trace_pos_idxs)]
  # FileIO.save("sample.jld2", "sample_st6",sample)
  sample = FileIO.load("sample.jld2", "sample_st6")
samples = [sample]
  X       = Flux.batch(smpl.x       for smpl in samples)
  A_mask  = Flux.batch(smpl.a_mask  for smpl in samples)
  As      = Flux.batch(smpl.as      for smpl in samples)
  Ps      = Flux.batch(smpl.ps      for smpl in samples)
  Vs      = Flux.batch(smpl.vs      for smpl in samples)
  Rs      = Flux.batch(smpl.rs      for smpl in samples)
  f32(arr) = convert(AbstractArray{Float32}, arr)
data =  map(x -> f32(x), (; X, A_mask, As, Ps, Vs, Rs))
f32(As)
(data.X, data.As)

  # data = [sample_batch(gspec, env.memory, env.params.learning_params) for _ in 1:1]
  losses(env.curnns, env.params.learning_params, data)
  @enter losses(env.curnns, env.params.learning_params, data)

  @enter sample_batch(gspec, envres.memory, envres.params.learning_params)
  samples = (sample_batch(gspec, envres.memory, envres.params.learning_params)|>gpu for _ in 1:n)
  for spl in samples
    @info spl
  end

  t = @timed self_play_step!(env)
  @timed learning_step!(env)
  # @enter self_play_step!(env)
  @enter learning_step!(env)

  # FileIO.save("env.jld2", env)
  mem_test = FileIO.load("memory_minmax.jld2", "mem")

  # for (b1,b2) in zip(env.memory, mem_test)
  #   @assert b1 ≂ b2
  # end
  # mem_test["mem"].buffer  env.memory.buffer
  # mem_test["mem"].buffer[1]
  # env.memory.buffer[1]

  # learning_step
lp = env.params.learning_params
tr = MuTrainer(env.gspec, env.curnns|>lp.device|>Flux.trainmode!, env.memory, env.params.learning_params, env.params.learning_params.opt)
nbatches = lp.batches_per_checkpoint
env.curnns|>Flux.gpu
update_weights!(tr, nbatches)
@timed samples = [sample_batch(tr.gspec, tr.memory, tr.hyper) for _ in 1:10] |> Flux.gpu

counter=0
for d in data
  counter += 1
end
counter
d = samples[1]
losses(tr.nns, tr.hyper, d)
@info d.X[:,:,1,1]

#TODO check with MCTS, MinMax
#Revise.jl, tensorboard_logger

# # checking Losses
# zero(Ps)
# one.(Ps)
# Flux.Losses.mse(zero(Vs), one.(Vs), agg=x->mean(∇_scale .* x)/sum(∇_scale))
# Flux.Losses.mse(zero(Vs), one.(Vs), agg=x->mean(x .* ∇_scale))
# Flux.Losses.mse(zero(Vs), one.(Vs), agg=x->sum(∇_scale .* x) / sum(∇_scale))
# Flux.Losses.mse(zero(Vs), one.(Vs))
# sum(∇_scale)
# one.(Vs)
# Vs
# one.(Ps) .* ∇_scale

# Flux.Losses.crossentropy(one.(Ps).*ℯ, one.(Ps), agg=x->mean(x,∇_scale))
# log.(one.(Ps) .* ℯ)
# log.(zero(Ps)+eps.(Ps))
# one.(Ps) .*
# one.(Ps) .* log.(one.(Ps).*ℯ)
# mean(.-sum(one.(Ps) .* log.(one.(Ps).*ℯ), dims=1), weights(∇_scale))
# ∇_scale .* 2

"Lenient comparison operator for `struct`, both mutable and immutable (type with \\eqsim)."
@generated function ≂(x, y)
  if !isempty(fieldnames(x)) && x == y
    mapreduce(n -> :(x.$n ≂ y.$n), (a,b)->:($a && $b), fieldnames(x))
  else
    :(x == y)
  end
end

# nns = deepcopy(tr.nns)
# nns == tr.nns
# nns ≂ tr.nns
# env.curnns ≂ env.bestnns
# nns = deepcopy(env.bestnns)
# nns ≂ env.bestnns
# fieldnames(typeof(nns))
# typeof(nns) == typeof(env.bestnns)
# nns.f.common

(X, A_mask, As, Ps, Vs, Rs) = sample_batch(env.gspec, env.memory, env.params.learning_params)
@enter losses(env.curnns, learning_params, (X, A_mask, As, Ps, Vs, Rs))
As
A = As[1:1,:]
Vs

env.curnns = deepcopy(env.bestnns)
@assert env.curnns ≂ env.bestnns    # all their fields are the same
@assert !(env.curnns == env.bestnns) # they have not the same memory adress

tr = MuTrainer(env.gspec, env.curnns, env.memory, env.params.learning_params, Flux.ADAM())
@assert env.curnns ≂ tr.nns
@assert env.curnns == tr.nns # same memory adress (shallow copy)

@timed update_weights!(tr, tr.hyper.batches_per_checkpoint)
@assert env.curnns == tr.nns
@assert !(env.bestnns ≂ tr.nns) # after update! tr.nns and env.curnns changed
@assert !(env.curnns ≂ env.bestnns)

learning_step!(env)
train!(env)

# gnerator with rand
function foo(x)
  y = x+1
  z = x+1
  return (; x,y,z)
end
rands = (rand() for _ in 1:5)
foos = (foo(r) for r in rands)
X = [f.x for f in foos] # Flux.batch(f.x for f in foos)
Y = [f.y for f in foos]
Z = [f.z for f in foos]
X[1], Y[1], Z[1]


text:
module StochasticOptimalTransport

using LogExpFunctions: LogExpFunctions

using LinearAlgebra: LinearAlgebra
using Random: Random
using Statistics: Statistics

include("utils.jl")
include("discrete.jl")
include("semidiscrete.jl")

@doc raw"""
    wasserstein([rng, ]c, μ, ν[, ε; kwargs...])

Estimate the (entropic regularization of the) Wasserstein distance
```math
W_{ε}(μ, ν) = \min_{π ∈ Π(μ,ν)} \int c(x, y) \,π(\mathrm{d}(x,y)) +
ε \mathrm{KL}(π \,|\, μ ⊗ ν)
```
with respect to cost function `c` using stochastic optimization.

If both measures `μ` and `ν` are [`DiscreteMeasure`](@ref)s, then the
Wasserstein distance is approximated with stochastic averaged gradient
descent (SAG). The convergence rate of SAG is ``O(1 / k)``, where
``k`` is the number of iterations, and hence converges faster than stochastic
gradient descent (SGD) at the price of increased memory consumption.

If only one of the measures `μ` and `ν` is a [`DiscreteMeasure`](@ref), then the
Wasserstein distance is approximated with stochastic gradient descent with
averaging (SGA). In this case, it is required that samples of the non-discrete measure `λ`
can be obtained with `rand(rng, λ)`. The convergence rate of SGA is ``O(1/√k)``,
where ``k`` is the number of iterations.

If `ε` is `nothing` (the default), then the unregularized Wasserstein distance is
approximated. Otherwise, the entropic regularization with `ε > 0` is estimated.

The SAG algorithm uses a constant step size, whereas the SGA algorithm uses the step size
schedule
```math
    τᵢ = \frac{τ₁}{1 + \sqrt{(i - 1) / w}}
```
for the ``i``th iteration, where ``τ₁`` corresponds to the initial step size and ``w``
indicates the number of iterations serving as a warm-up phase.

# Keyword arguments
- `maxiters::Int = 10_000`: maximum number of gradient steps
- `stepsize = 1`: constant stepsize of SAG or initial step size ``τ₁`` of SGA
- `warmup_phase = 1`: warm-up phase ``w`` of SGA
- `atol = 0`: absolute tolerance
- `rtol = iszero(atol) ? typeof(float(atol))(1 // 10_000) : 0`: relative tolerance
- `montecarlo_samples = 10_000`: Number of Monte Carlo samples from `μ` for approximating
  an expectation with respect to `μ` in the semi-discrete optimal transport problem

# References

Genevay et al. (2016). Stochastic Optimization for Large-Scale Optimal Transport. Advances in Neural Information Processing Systems (NIPS 2016), 29:3440-3448.

Peyré, Gabriel, & Marco Cuturi (2019). Computational Optimal Transport. Foundations and Trends in Machine Learning, 11(5-6):355-607.
"""
wasserstein(args...; kwargs...) = wasserstein(Random.GLOBAL_RNG, args...; kwargs...)
function wasserstein(
    rng::Random.AbstractRNG, c, μ, ν, ε::Union{Real,Nothing}=nothing; kwargs...
)
    # approximate solution `v` of the dual problem
    v = dual_v(rng, c, μ, ν, ε; kwargs...)

    # compute Wasserstein distance from dual solution
    cost = dual_cost(rng, c, v, μ, ν, ε; kwargs...)

    return cost
end

end


text:
using QuantumOptics: NLevelBasis, CompositeBasis, FockBasis
using Test, IonSim
using IonSim.PhysicalConstants: ħ, ca40_qubit_transition_frequency, c, m_ca40
using Suppressor

@suppress_err begin

# setup system
C = Ca40()
L1 = Laser(); L2 = Laser()
chain = LinearChain(
            ions=[C, C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
        )

@testset "traps -- Trap" begin
    T = Trap(configuration=chain, lasers=[L1, L1])

    # test construction of CompositeBasis
    @test T.basis.bases[1] ≡ T.configuration.ions[1]
    @test T.basis.bases[2] ≡ T.configuration.ions[2]
    @test T.basis.bases[3] ≡ T.configuration.vibrational_modes.z[1]

    # test construction of T.δB
    t = 0:1:100
    ones = [1 for _ in t]
    T.δB = 1
    @test T.δB.(t) == ones
    @test T._cnst_δB
    T.δB = sin
    @test T.δB.(t) == sin.(t)
    @test !T._cnst_δB

    # test for warning when lasers=[L, L, L, ...] where L point to the same thing
    warn = "Some lasers point to the same thing. Making copies."
    @test_logs (:warn, warn) Trap(configuration=chain, lasers=[L1, L1, L1])
    # and test that, in this case, copies are made
    T1 = Trap(configuration=chain, lasers=[L1, L1, L1, L1], δB=t->t)
    for i in 1:4, j in i+1:4
        @test !(T1.lasers[i] ≡ T1.lasers[j])
    end

    # test setproperty for Trap
    # setting T.δB to a constant
    T.δB = 100
    @test T.δB(0) == 100
    # setting the :basis directly should have no effect
    basis = T.basis
    T.basis = 100
    @test T.basis == basis
    # changing the configuration should also update the basis
    chain1 = LinearChain(
            ions=[C, C, C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
        )
    @test length(T.basis.bases) ≡ 3
    T.configuration = chain1
    @test T.configuration ≡ chain1
    @test length(T.basis.bases) ≡ 4

    # make sure print/show doesn't throw errors
    print(T); show(T)
end

@testset "traps -- general functions" begin
    Bhat = (x̂ + ŷ + ẑ)/√3
    T = Trap(configuration=chain, lasers=[L1], Bhat=Bhat)

    # get_basis
    @test T.basis == get_basis(T)

    # global_beam!
    L = Laser()
    global_beam!(T, L)
    @test L.pointing == [(1, 1.0), (2, 1.0)]

    # Efield_from_pi_time
    ion_index = 1
    laser_index = 1
    ion = T.configuration.ions[ion_index]
    transition = ("S-1/2", "D-1/2")
    # compare to specific pre-computed value for both methods
    E1 = Efield_from_pi_time(1e-6, Bhat, L1, ion, transition)
    E2 = Efield_from_pi_time(1e-6, T, laser_index, ion_index, transition)
    @test E1 ≈ 153739.56 rtol=1e-2
    @test E1 == E2
    # confirm in-place versions work
    Efield_from_pi_time!(1e-6, Bhat, L1, ion, transition)
    @test L1.E(0) == E1
    Efield_from_pi_time!(1e-6, T, laser_index, ion_index, transition)
    @test L1.E(0) == E1
    # shouldn't be able to have a laser argument where laser.pointing = []
    L = Laser()
    @test_throws AssertionError Efield_from_pi_time(1e-6, Bhat, L, ion, transition)
    L.pointing = [(1, 1.)]
    @test isinf(Efield_from_pi_time(1e-6, x̂, L, ion, transition))

    # Efield_from_rabi_frequency
    ion_index = 1
    laser_index = 1
    ion = T.configuration.ions[ion_index]
    transition = ("S-1/2", "D-1/2")
    # compare to specific pre-computed value for both methods
    E1 = Efield_from_rabi_frequency(5e5, Bhat, L1, ion, transition)
    E2 = Efield_from_rabi_frequency(5e5, T, laser_index, ion_index, transition)
    @test E1 ≈ 153739.56 rtol=1e-2
    @test E1 == E2
    # confirm in-place versions work
    Efield_from_rabi_frequency!(5e5, Bhat, L1, ion, transition)
    @test L1.E(0) == E1
    Efield_from_rabi_frequency!(5e5, T, laser_index, ion_index, transition)
    @test L1.E(0) == E1
    # shouldn't be able to have a laser argument where laser.pointing = []
    L = Laser()
    @test_throws AssertionError Efield_from_rabi_frequency(5e5, Bhat, L, ion, transition)

    # transition_frequency (test against pre-computed values)
    T.B = 4e-4
    f = transition_frequency(T, C, transition)
    @test f ≈ 2.2393e6 rtol=1e-4
    @test transition_frequency(T.B, C, transition) == f
    @test transition_frequency(T, 1, transition) == f

    # set_gradient
    set_gradient!(T, (1, 2), transition, 1e6)
    f1 = transition_frequency(T, T.configuration.ions[1], transition)
    f2 = transition_frequency(T, T.configuration.ions[2], transition)
    @test abs(f1 - f2) ≈ 1e6

    # test :(==)
    chain1 = LinearChain(
            ions=[C, C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(x=[1], z=[1])
        )
    T = Trap(configuration=chain1, lasers=[L1, L2])
    xmode = chain1.vibrational_modes.x[1]
    zmode = chain1.vibrational_modes.z[1]
    cb = get_basis(T)
    @test cb == C ⊗ C ⊗ xmode ⊗ zmode
    @test cb != C ⊗ C ⊗ zmode ⊗ xmode
    @test cb != C ⊗ C ⊗ C ⊗ xmode ⊗ zmode

    # test get_η
    f = ca40_qubit_transition_frequency
    λ = c/f
    # η = |k|cos(θ) * √(ħ / (2M ⋅ N ⋅ 2πν)); cos(θ) ≡ k̂ ⋅ mode_axis; N ≡ number of ions
    η(ν) = (2π/λ) * sqrt(ħ / (2 * m_ca40 * 2 * 2π * ν))
    @test abs(get_η(zmode, L, C)) ≈ η(zmode.ν)
    L.k = (x̂ + ẑ) / √2
    @test abs(get_η(xmode, L, C)) ≈ η(xmode.ν) / √2
end
end  # end suppress


text:
using ClusterValidityIndices
using Test
using Logging

# Set the log level
LogLevel(Logging.Info)

# Include the test utilities
include("utils.jl")

@testset "ClusterValidityIndices.jl" begin
    # Run all of the CVI tests
    include("test_cvis.jl")
end

@testset "Utils" begin
    @info "Sorting CVI data"
    # Load the sample data
    data_paths = readdir("../data", join=true)
    local_data, local_labels = get_cvi_data(data_paths[1])
    # Sort the data (reordering both data and labels monotonically)
    local_data, local_labels = sort_cvi_data(local_data, local_labels)
    # Verify that the labels are monotonic iteratively
    is_monotonic = true
    for i = 1:(length(local_labels) - 1)
        if local_labels[i] > local_labels[i + 1]
            is_monotonic = false
        end
    end
    @test is_monotonic
end

@testset "Label Map" begin
    @info "Label Map Testing"

    # Create the ICVI
    cvi = CH()

    # Load the sample data
    data_paths = readdir("../data", join=true)
    local_data, local_labels = get_cvi_data(data_paths[1])

    # Permute the data
    n_data = length(local_labels)
    indices = randperm(n_data)
    local_data = local_data[:, indices]
    local_labels = local_labels[indices]

    # Incrementally compute the ICVI
    for ix = 1:length(local_labels)
        sample = local_data[:, ix]
        label = local_labels[ix]
        _ = get_icvi!(cvi, sample, label)
    end

    @info "Task map: " cvi.label_map
end

@testset "Example Scripts" begin
    # Switch to the top for execution because our scripts point to the datasets
    # relative to themselves, not relative to the test dir
    @info "Test directory" pwd()
    cd("../")
    @info "Switching working directory to top for running scripts" pwd()

    # Run scripts
    include("../examples/cvi_example.jl")
end


text:
using EmpiricalBayes
using Test
using DelimitedFiles

using Distributions
using Distances

data_folder_path = joinpath(dirname(@__FILE__), "data")

struct TestDist <: ContinuousUnivariateDistribution
    x::Int
end

Distributions.pdf(d::TestDist, n::Real) = d.x
Distributions.mode(d::TestDist) = d.x

############################# Discretization ###########################

@testset "Discretization Tests" begin

# Basic functionality
test_test_statistics = [1, 1.5, 1.2, 2, 2.3, 2.7, 3, 3.3, 3.9, 4]
mids, counts, bin_width = discretize_test_statistics(test_test_statistics, 3)
@test mids == [1.5, 2.5, 3.5]
@test counts == [3, 3, 4]
@test bin_width == 1

@testset "Discretization Random Tests" begin
# Random tests
num_random_tests = 100
for _ in 1:num_random_tests
    bin_width = rand() * 100 + eps()
    num_bins = Int(round(rand()*100)) + 1
    min_ts = rand() * 10
    max_ts = min_ts + bin_width * num_bins

    target_counts = [Int(round(rand()*100)) for i in 1:num_bins]
    target_mids = [min_ts + bin_width/2 + (i-1)*bin_width for i in 1:num_bins]

    rand_test_stats = [[bin_width * (j-1) + min_ts + eps() + rand() * bin_width * 0.9 for i in 1:target_counts[j]] for j in 1:length(target_counts)]
    rand_test_stats = reduce(vcat, rand_test_stats)
    rand_test_stats = [float(x) for x in rand_test_stats]
    push!(rand_test_stats, min_ts)
    push!(rand_test_stats, max_ts)
    target_counts[1] += 1
    target_counts[end] += 1

    rmids, rcounts, rwidths = discretize_test_statistics(rand_test_stats, num_bins)

    @test isapprox(rmids, target_mids, atol=0.00001)
    @test rcounts == target_counts
    @test isapprox(rwidths, bin_width, atol=0.00001)
end
end

# Test silent type conversion
mids, counts, bin_width = discretize_test_statistics(test_test_statistics, 3)
fmids, fcounts, fwidth = discretize_test_statistics(test_test_statistics, 3.0)
@test fmids == mids
@test fcounts == counts
@test fwidth == bin_width

# Test invalid arguments
@test_throws ErrorException discretize_test_statistics(test_test_statistics, 3.3)
@test_throws MethodError discretize_test_statistics("1, 1.5, 1.2, 2, 2.3, 2.7, 3, 3.3, 3.9, 4", 3)

end

############################# Truncating ###########################

@testset "Truncate histogram tests" begin

######## Basic Test
mids, counts, bin_width = discretize_test_statistics(collect(Float64, 1:100), 10)

proportions = [0.89, 0.9, 0.91, 1.0]
targets = [8, 9, 9, 10]

for i in 1:length(proportions)
    @test EmpiricalBayes.keep_bottom_proportion_of_hist(mids, counts, proportions[i], true) == (mids[1:targets[i]], counts[1:targets[i]])
end

end

############################# Fit Null ###############################

@testset "Null Fit Tests" begin

######## Basic test

# Gamma distribution
mids, counts, width = discretize_test_statistics([1, 1.5, 1.2, 2, 2,2,2,2,2, 2.3, 2.7, 3, 3.3, 3.9, 4], 3)

dist = fit_null_distribution(mids, counts, 3, width, 1.0, Gamma, verbose = false)
f(x) = pdf(dist, x)

basic_test_result = f.(0:0.1:10)

basic_test_target = [0.0       ,2.32616e-10,1.22024e-7 ,4.043e-6   ,4.32361e-5 ,0.000248819,0.000967622,0.0028712  ,0.00698951 ,0.0146273  ,0.0271695  ,0.0458241  ,0.0713679  ,0.103957   ,0.143041   ,0.187386   ,0.235202   ,0.284332   ,0.332473   ,0.377397   ,0.417133   ,0.450104   ,0.47521    ,0.491851   ,0.499911   ,0.499696   ,0.491862   ,0.477317   ,0.457138   ,0.43248    ,0.404504   ,0.374322   ,0.342946   ,0.311267   ,0.280033   ,0.249848   ,0.221176   ,0.194349   ,0.169582   ,0.14699,0.126606   ,0.108396   ,0.0922763  ,0.0781275  ,0.0658054  ,0.0551525  ,0.0460054  ,0.0382016  ,0.0315839  ,0.026004   ,0.0213245  ,0.01742    ,0.014178   ,0.0114984  ,0.00929344 ,0.00748661 ,0.00601195 ,0.004813   ,0.00384178 ,0.0030578  ,0.00242711 ,0.00192136 ,0.00151708 ,0.00119488 ,0.000938834,0.000735927,0.000575565,0.000449155,0.000349759,0.000271794,0.000210783,0.000163147,0.000126036,9.71871e-5 ,7.48064e-5 ,5.74787e-5 ,4.40894e-5 ,3.37629e-5 ,2.58132e-5 ,1.97041e-5 ,1.50178e-5,1.14288e-5,8.68487e-6,6.59029e-6,4.9939e-6 ,3.77907e-6,2.85596e-6,2.15554e-6,1.62483e-6,1.22327e-6,9.19837e-7,6.9085e-7 ,5.18267e-7,3.88356e-7,2.90686e-7,2.17344e-7,1.62333e-7,1.2112e-7 ,9.02779e-8,6.72221e-8,5.00055e-8 ]

@test isapprox(basic_test_result, basic_test_target, atol=0.00001)

# Normal distribution

# Test statistics generated from Normal(0, 1) distribution
ts = randn(1000)
mids, counts, width = discretize_test_statistics(ts, 10)
dist = fit_null_distribution(mids, counts, 10, width, 1.0, Normal, verbose = false)
@test isapprox(mean(dist), 0, atol = 0.5)
@test isapprox(std(dist), 1, atol = 0.5)

######## Random Tests
@testset "Null Random Tests" begin

num_random_tests = 300
hellinger_sum = 0
for _ in 1:num_random_tests
    α = Int(round(rand()*10)) + 3
    β = Int(round(rand()*10)) + 0.01
    reference_distr = Gamma(α, β)
    rand_test_stats = rand(reference_distr, 1000)

    num_bins = 3 + Int(round(rand()*100))
    rmids, rcounts, rwidths = discretize_test_statistics(rand_test_stats, num_bins)

    test_distr = fit_null_distribution(rmids, rcounts, num_bins, rwidths, 1.0, Gamma, verbose = false)
    test_pdf(x) = pdf(test_distr, x)

    ref_values = [pdf(reference_distr, x) for x in 0:0.01:20]
    test_values = [test_pdf(x) for x in 0:0.01:20]
    hellinger_dist = hellinger(ref_values, test_values)

    @test hellinger_dist < 0.5

    hellinger_sum += hellinger_dist
end

hellinger_sum = hellinger_sum / num_random_tests
@test hellinger_sum <= 0.3

end

######## Test invalid input error
# Test error out on Gamma inputs where α is small
gamma_distr = Gamma(0.001, 5.0)
ref_test_stats = rand(gamma_distr, 100)
tmids, tcounts, twidths = discretize_test_statistics(ref_test_stats, 10)
@test_throws ErrorException fit_null_distribution(tmids, tcounts, 10, twidths, 1.0, Gamma, verbose = false)

end

############################# Fit Mixture ############################
@testset "Mixture Fit Tests" begin

######## Basic Test - fit spline to standard normal
ref_distr = Normal(0,1);
normal_samples_filepath = joinpath(data_folder_path, "normal_01_1000_samples.txt")
test_stats = readdlm(normal_samples_filepath)
test_mids, test_counts, test_width = discretize_test_statistics(test_stats, 10);
fh = fit_mixture_distribution(test_mids, test_counts, test_width);
test_vals = [fh(x) for x in -5:0.01:5];
ref_vals = [pdf(ref_distr, x) for x in -5:0.01:5];
hellinger_dist = hellinger(test_vals, ref_vals)
@test hellinger_dist < 0.2

end

############################# Calculate Posterior ####################
@testset "Posterior Tests" begin

######## Basic Test
test_post = calculate_posterior([0, 0], [0, 1], TestDist(1), y->1, :two)
@test test_post ≈ [0.5, 0.7310585786300049] atol=0.0001

test_post = calculate_posterior([0, 0], [0, 1], TestDist(2), y->1, :two)
@test test_post ≈ [0.0, 0.4621171572600098] atol=0.0001

test_post = calculate_posterior([1, 2], [0, 1], TestDist(1), y->1, :two)
@test test_post ≈ [0.5, 0.7310585786300049] atol=0.0001

test_post = calculate_posterior([1, 2], [0, 1], TestDist(1), y->0, :two)
@test test_post ≈ [0.0, 0.0] atol=0.0001

test_priors = [0, 1, 1, 0, 1]
test_statistics = [1, 2, 3, 4, 5]
test_null_dist = Normal(2, 5)
test_mix_pdf(x) = pdf(Normal(3, 5), x)
reference_post = [0.469082, 0.725626, 0.736384, 0.529118, 0.756652]
test_post = calculate_posterior(test_statistics, test_priors, test_null_dist, test_mix_pdf, :two)
@test test_post ≈ reference_post atol=0.000001

# No prior
test_post = calculate_posterior([0, 0], TestDist(1), y->1, :two)
@test test_post ≈ [0.5, 0.5] atol=0.0001

# One-tailed test, upper tail
reference_post = [-Inf, -Inf, 0.736384, 0.529118, 0.756652]
test_post = calculate_posterior(test_statistics, test_priors, test_null_dist, test_mix_pdf, :upper)
@test test_post ≈ reference_post atol=0.000001

# One-tailed test, lower tail
reference_post = [0.469082, -Inf, -Inf, -Inf, -Inf]
test_post = calculate_posterior(test_statistics, test_priors, test_null_dist, test_mix_pdf, :lower)
@test test_post ≈ reference_post atol=0.000001

end

############################# empirical_bayes ####################
@testset "Pipeline Tests" begin
# Test that the function at least runs
gamma_stats_filepath = joinpath(data_folder_path, "gamma_51_1000_samples.txt")
gamma_stats = readdlm(gamma_stats_filepath)

priors_filepath = joinpath(data_folder_path, "priors_random_1000_samples.txt")
priors = readdlm(priors_filepath)

eb_at_least_runs = true
try
eb = empirical_bayes(gamma_stats, priors, 10.0, Gamma)
catch
eb_at_least_runs = false
end
@test eb_at_least_runs
# With symbol for distribution
try
eb = empirical_bayes(gamma_stats, priors, 10.0, :Gamma)
catch
eb_at_least_runs = false
end
@test eb_at_least_runs

# Test that it runs without priors
no_priors_eb_at_least_runs = true
try
eb = empirical_bayes(gamma_stats, 10.0, Gamma)
catch
no_priors_eb_at_least_runs = false
end
@test no_priors_eb_at_least_runs
# With symbol for distribution
try
eb = empirical_bayes(gamma_stats, 10.0, :Gamma)
catch
no_priors_eb_at_least_runs = false
end
@test no_priors_eb_at_least_runs

end


text:
export ft_stopintask

function ft_stopintask(ft_handle::Culong)
  ft_status = ccall((:FT_StopInTask, d2xx),Cuint,(Culong,),ft_handle)
  checkstatus(ft_status)
  return nothing
end


text:
function sphere_texture!(vis::Visualizer, mech::Mechanism;
        imagepath = joinpath(@__DIR__, "dojo_texture.png"), name=:robot)
    image = PngImage(imagepath)
    texture = Texture(image=image);
    material = MeshLambertMaterial(map=texture);
    setobject!(vis[name]["bodies"]["body:1"], convert_shape(mech.bodies.values[1].shape), material)
    return nothing
end


text:
include("pulaymix.jl")

function main_debug()
    Npoints = 100
    rho = rand(Npoints)
    rhonew = rand(Npoints)
    beta = 0.5
    MIXDIM = 4
    dvmat = rand(Float64,Npoints,MIXDIM)
    vmat  = rand(Float64,Npoints,MIXDIM)
    iter = 5
    rho = pulaymix!( rho, rhonew, beta, dvmat, vmat, iter, MIXDIM)
end

main_debug()


text:
using VerTeX
using Test

# write your own tests here
@test (x = VerTeX.article("",VerTeX.preamble()*"\\author{x}\n\\date{x}\n\\title{x}\n"); typeof(x) == typeof(VerTeX.dict2tex(VerTeX.tex2dict(x))))


text:
module TensorOperations

export tensorcopy, tensoradd, tensortrace, tensorcontract, tensorproduct, scalar
export tensorcopy!, tensoradd!, tensortrace!, tensorcontract!, tensorproduct!

export @tensor

# Auxiliary functions
#---------------------
include("auxiliary/axpby.jl")
include("auxiliary/error.jl")
include("auxiliary/meta.jl")
include("auxiliary/stridedarray.jl")
include("auxiliary/strideddata.jl")
include("auxiliary/unique2.jl")

# Implementations
#-----------------
include("implementation/indices.jl")
include("implementation/kernels.jl")
include("implementation/recursive.jl")
include("implementation/stridedarray.jl")
include("implementation/strides.jl")

# Index notation
#----------------
include("indexnotation/tensormacro.jl")
include("indexnotation/indexedobject.jl")
include("indexnotation/sum.jl")
include("indexnotation/product.jl")

# Functions
#----------
include("functions/simple.jl")
include("functions/inplace.jl")

end # module


text:
using FDM: Forward, Backward, Central, Nonstandard

@testset "Methods" begin
    for f in [:forward_fdm, :backward_fdm, :central_fdm]
        @eval @test $f(10, 1; bound=1)(sin, 1) ≈ cos(1)
        @eval @test $f(10, 2; bound=1)(sin, 1) ≈ -sin(1)

        @eval @test $f(10, 1; bound=1)(exp, 1) ≈ exp(1)
        @eval @test $f(10, 2; bound=1)(exp, 1) ≈ exp(1)

        @eval @test $f(10, 1; bound=1)(abs2, 1) ≈ 2
        @eval @test $f(10, 2; bound=1)(abs2, 1) ≈ 2

        @eval @test $f(10, 1; bound=1)(sqrt, 1) ≈ .5
        @eval @test $f(10, 2; bound=1)(sqrt, 1) ≈ -.25
    end

    @testset "Adaptation improves estimate" begin
        @test forward_fdm(5, 1)(log, 0.001; adapt=0) ≈ 969.2571703
        @test forward_fdm(5, 1)(log, 0.001; adapt=1) ≈ 1000
    end

    @testset "Limiting step size" begin
        @test !isfinite(central_fdm(5, 1)(abs, 0.001; max_step=0))
        @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0
    end

    @testset "Printing FDMethods" begin
        @test sprint(show, central_fdm(2, 1)) == """
            FDMethod:
              order of method:       2
              order of derivative:   1
              grid:                  [-1, 1]
              coefficients:          [-0.5, 0.5]
            """
        m, _ = fdm(central_fdm(2, 1), sin, 1, Val(true))
        report = sprint(show, m)
        regex_float = r"[\d\.\+-e]+"
        regex_array = r"\[([\d.+-e]+(, )?)+\]"
        @test occursin(Regex(join(map(x -> x.pattern,
            [
                r"FDMethod:",
                r"order of method:", r"\d+",
                r"order of derivative:", r"\d+",
                r"grid:", regex_array,
                r"coefficients:", regex_array,
                r"roundoff error:", regex_float,
                r"bounds on derivatives:", regex_float,
                r"step size:", regex_float,
                r"accuracy:", regex_float,
                r""
            ]
        ), r"\s*".pattern)), report)
    end

    @testset "Breaking deprecations" begin
        @test_throws ErrorException fdm([1,2,3], 4)  # Custom grids need Nonstandard
        for f in (forward_fdm, backward_fdm, central_fdm)
            @test_throws ErrorException f(2, 1; M=1)  # Old kwarg, now misplaced
            @test_throws ErrorException f(2, 1, Val(true))  # Ask fdm for reports instead
        end
    end

    @testset "Types" begin
        @testset "$T" for T in (Forward, Backward, Central)
            @test T(5, 1)(sin, 1; adapt=4) ≈ cos(1)
            @test_throws ArgumentError T(3, 4)
            @test_throws ArgumentError T(40, 5)
            @test_throws ArgumentError T(5, 1)(sin, 1; adapt=200)
            @test_throws ArgumentError T(5, 1)(sin, 1; eps=0.0)
            @test_throws ArgumentError T(5, 1)(sin, 1; bound=0.0)
        end
        @testset "Nonstandard" begin
            @test Nonstandard([-2, -1, 1], 1)(sin, 1) ≈ cos(1)
            @test_throws ArgumentError Nonstandard([-2, -1, 1], 1)(sin, 1; adapt=2)
        end
    end
end


text:
using QuranTree
using Test
using JuliaDB: select
using PrettyTables: pretty_table
using Suppressor: @capture_out

function capture_io(x)
    io = IOBuffer()
    println(io, x)
    return String(take!(io))
end

@testset "QuranTree.jl" begin
    data = QuranData()
    crps, tnzl = load(data)

    # data
    @test crps[end] === "(114:6:3:3)\tn~aAsi\tN\tSTEM|POS:N|LEM:n~aAs|ROOT:nws|MP|GEN"
    @test tnzl[begin] === "1|1|بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ"
    @test typeof(crps) === CorpusRaw
    @test typeof(tnzl) === TanzilRaw

    # indexing
    crpsdata = table(crps)
    tnzldata = table(tnzl)

    @test select(crpsdata[begin].data, :form)[1] === "bi"
    @test select(crpsdata[114][2][1][1].data, :features)[1] === "STEM|POS:N|LEM:malik|ROOT:mlk|MS|GEN"
    @test select(crpsdata[114][2][1:2][1].data, :features)[2] === "PREFIX|Al+"
    @test select(crpsdata[114][4][[1,2,3]][1].data, :features)[2] === "STEM|POS:N|LEM:\$ar~|ROOT:\$rr|MS|GEN"

    @test select(crpsdata[114][1:3][1][1].data, :features)[3] === "STEM|POS:N|LEM:<ila`h|ROOT:Alh|MS|GEN"
    @test select(crpsdata[112][1:3][1:3][2].data, :features)[1] === "STEM|POS:N|LEM:S~amad|ROOT:Smd|MS|NOM"
    @test select(crpsdata[112][1:3][[2,4]][1].data, :features)[3] === "PREFIX|Al+"

    @test select(crpsdata[111][[2,3]][1][1].data, :features)[1] === "STEM|POS:NEG|LEM:maA"
    @test select(crpsdata[111][[2,3]][1:3][2].data, :features)[2] === "STEM|POS:V|IMPF|LEM:yaSolaY|ROOT:Sly|3MS"
    @test select(crpsdata[111][[2,3]][[1,3]][1].data, :features)[4] === "STEM|POS:N|LEM:*uw|FS|ACC"

    @test select(crpsdata[111:112][1][1][1].data, :features)[2] === "STEM|POS:V|IMPV|LEM:qaAla|ROOT:qwl|2MS"
    @test select(crpsdata[111:112][1][1:3][1].data, :features)[4] === "STEM|POS:V|IMPV|LEM:qaAla|ROOT:qwl|2MS"
    @test select(crpsdata[111:112][1][[1,2,3]][1].data, :features)[4] === "STEM|POS:V|IMPV|LEM:qaAla|ROOT:qwl|2MS"

    @test select(crpsdata[111:112][1:3][1][1].data, :features)[5] === "STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|NOM"
    @test select(crpsdata[111:112][1:3][1:3][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"
    @test select(crpsdata[111:112][1:3][[1,2,3]][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"

    @test select(crpsdata[111:112][[1,2,3]][1][1].data, :features)[5] === "STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|NOM"
    @test select(crpsdata[111:112][[1,2,3]][1:3][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"
    @test select(crpsdata[111:112][[1,2,3]][[1,2,3]][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"

    @test select(crpsdata[[111,112]][1][1][1].data, :features)[2] === "STEM|POS:V|IMPV|LEM:qaAla|ROOT:qwl|2MS"
    @test select(crpsdata[[111,112]][1][1:3][1].data, :features)[4] === "STEM|POS:V|IMPV|LEM:qaAla|ROOT:qwl|2MS"
    @test select(crpsdata[[111,112]][1][[1,2,3]][1].data, :features)[4] === "STEM|POS:V|IMPV|LEM:qaAla|ROOT:qwl|2MS"

    @test select(crpsdata[[111,112]][1:3][1][1].data, :features)[5] === "STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|NOM"
    @test select(crpsdata[[111,112]][1:3][1:3][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"
    @test select(crpsdata[[111,112]][1:3][[1,2,3]][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"

    @test select(crpsdata[[111,112]][[1,2,3]][1][1].data, :features)[5] === "STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|NOM"
    @test select(crpsdata[[111,112]][[1,2,3]][1:3][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"
    @test select(crpsdata[[111,112]][[1,2,3]][[1,2,3]][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"

    @test select(crpsdata[end-3:end-2][[1,2,3]][1][1].data, :features)[5] === "STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|NOM"
    @test select(crpsdata[end-3:end-2][[1,2,3]][1:3][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"
    @test select(crpsdata[end-3:end-2][[1,2,3]][[1,2,3]][1].data, :features)[9] === "STEM|POS:N|LEM:*uw|FS|ACC"

    @test select(tnzldata[1].data, :form)[7] === "صِرَٰطَ ٱلَّذِينَ أَنْعَمْتَ عَلَيْهِمْ غَيْرِ ٱلْمَغْضُوبِ عَلَيْهِمْ وَلَا ٱلضَّآلِّينَ"
    @test select(tnzldata[111:113].data, :form)[10] === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ قُلْ أَعُوذُ بِرَبِّ ٱلْفَلَقِ"
    @test select(tnzldata[[111,112,113]].data, :form)[10] === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ قُلْ أَعُوذُ بِرَبِّ ٱلْفَلَقِ"
    @test select(tnzldata[end-3:end-1].data, :form)[10] === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ قُلْ أَعُوذُ بِرَبِّ ٱلْفَلَقِ"

    # verses
    @test verses(crpsdata[1])[7] === "Sira`Ta {l~a*iyna >anoEamota Ealayohimo gayori {lomagoDuwbi Ealayohimo walaA {lD~aA^l~iyna"
    @test verses(crpsdata[1][1:7])[7] === "Sira`Ta {l~a*iyna >anoEamota Ealayohimo gayori {lomagoDuwbi Ealayohimo walaA {lD~aA^l~iyna"
    @test verses(crpsdata[1][7])[1] === "Sira`Ta {l~a*iyna >anoEamota Ealayohimo gayori {lomagoDuwbi Ealayohimo walaA {lD~aA^l~iyna"
    @test verses(tnzldata[1][1])[1] === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ"
    @test verses(tnzldata[1][1:2])[2] === "ٱلْحَمْدُ لِلَّهِ رَبِّ ٱلْعَٰلَمِينَ"
    @test verses(tnzldata)[1] === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ"
    # arabic
    @test arabic(verses(crpsdata[114])[1]) === "قُلْ أَعُوذُ بِرَبِّ ٱلنَّاسِ"
    @test arabic(verses(crpsdata[1][7])[1]) === "صِرَٰطَ ٱلَّذِينَ أَنْعَمْتَ عَلَيْهِمْ غَيْرِ ٱلْمَغْضُوبِ عَلَيْهِمْ وَلَا ٱلضَّآلِّينَ"
    # number
    @test verses(crpsdata[113:114], number=true)[1] === "113:(1,5)"
    @test verses(crpsdata[113:114], number=true, start_end=false)[1] == ([113], [1, 2, 3, 4, 5])

    # words
    @test words(tnzldata[1][7])[1] === "صِرَٰطَ"

    # encoding
    @test encode(arabic(verses(crpsdata[1][7])[1])) === verses(crpsdata[1][7])[1]

    # chapter_name
    @test chapter_name(crpsdata[13][2][1]) === "ٱلرَّعْد"
    @test chapter_name(crpsdata[13][2][1], lang=:english) === "Thunder"
    @test chapter_name(tnzldata[13][2]) === "ٱلرَّعْد"
    @test chapter_name(tnzldata[13][2], lang=:english) === "Thunder"
    @test chapter_name(tnzldata[1], true) === "{lofaAtiHap"
    @test chapter_name(tnzldata[1], true, lang=:english) === "The Opening"

    # dediac
    @test dediac(verses(crpsdata[1][1])[1]) === "bsm {llh {lrHm`n {lrHym"
    @test dediac(arabic(verses(crpsdata[1][1])[1])) === "بسم ٱلله ٱلرحمٰن ٱلرحيم"
    @test dediac(arabic(verses(crpsdata[1][1])[1])) === arabic(dediac(verses(crpsdata[1][1])[1]))

    # normalize
    @test normalize(dediac(verses(crpsdata[1][1])[1])) === "bsm Allh AlrHmAn AlrHym"
    @test normalize(dediac(verses(crpsdata[1][1])[1])) === dediac(normalize(verses(crpsdata[1][1])[1]))
    @test normalize(arabic(verses(crpsdata[1][1])[1]), :alif_khanjareeya) === "بِسْمِ ٱللَّهِ ٱلرَّحْمَانِ ٱلرَّحِيمِ"
    @test normalize(arabic(verses(crpsdata[1][1])[1]), :hamzat_wasl) === "بِسْمِ اللَّهِ الرَّحْمَٰنِ الرَّحِيمِ"
    @test normalize(arabic(verses(crpsdata[2][4])[1]), :alif_maddah) === "وَٱلَّذِينَ يُؤْمِنُونَ بِمَا أُنزِلَ إِلَيْكَ وَمَا أُنزِلَ مِن قَبْلِكَ وَبِٱلْءَاخِرَةِ هُمْ يُوقِنُونَ"
    @test normalize(arabic(verses(crpsdata[2][4])[1]), :alif_hamza_above) === "وَٱلَّذِينَ يُؤْمِنُونَ بِمَآ اُنزِلَ إِلَيْكَ وَمَآ اُنزِلَ مِن قَبْلِكَ وَبِٱلْءَاخِرَةِ هُمْ يُوقِنُونَ"
    @test normalize(arabic(verses(crpsdata[1][5])[1]), :alif_hamza_below) === "اِيَّاكَ نَعْبُدُ وَاِيَّاكَ نَسْتَعِينُ"
    @test normalize(arabic(verses(crpsdata[2][3])[1]), :waw_hamza_above) === "ٱلَّذِينَ يُوْمِنُونَ بِٱلْغَيْبِ وَيُقِيمُونَ ٱلصَّلَوٰةَ وَمِمَّا رَزَقْنَٰهُمْ يُنفِقُونَ"
    @test normalize(arabic(verses(crpsdata[2][15])[1]), :ya_hamza_above) === "ٱللَّهُ يَسْتَهْزِيُ بِهِمْ وَيَمُدُّهُمْ فِى طُغْيَٰنِهِمْ يَعْمَهُونَ"
    @test normalize(arabic(verses(crpsdata[2][2])[1]), :alif_maksura) === "ذَٰلِكَ ٱلْكِتَٰبُ لَا رَيْبَ فِيهِ هُدًي لِّلْمُتَّقِينَ"
    @test normalize(arabic(verses(crpsdata[2][3])[1]), :ta_marbuta) === "ٱلَّذِينَ يُؤْمِنُونَ بِٱلْغَيْبِ وَيُقِيمُونَ ٱلصَّلَوٰهَ وَمِمَّا رَزَقْنَٰهُمْ يُنفِقُونَ"
    @test normalize(arabic(verses(crpsdata[1][1])[1]), [:alif_khanjareeya, :hamzat_wasl]) === "بِسْمِ اللَّهِ الرَّحْمَانِ الرَّحِيمِ"

    @test normalize(verses(tnzldata[1][1])[1], :alif_khanjareeya) === "بِسْمِ ٱللَّهِ ٱلرَّحْمَانِ ٱلرَّحِيمِ"
    @test normalize(verses(tnzldata[1][1])[1], :hamzat_wasl) === "بِسْمِ اللَّهِ الرَّحْمَٰنِ الرَّحِيمِ"
    @test normalize(verses(tnzldata[2][4])[1], :alif_maddah) === "وَٱلَّذِينَ يُؤْمِنُونَ بِمَا أُنزِلَ إِلَيْكَ وَمَا أُنزِلَ مِن قَبْلِكَ وَبِٱلْءَاخِرَةِ هُمْ يُوقِنُونَ"
    @test normalize(verses(tnzldata[2][4])[1], :alif_hamza_above) === "وَٱلَّذِينَ يُؤْمِنُونَ بِمَآ اُنزِلَ إِلَيْكَ وَمَآ اُنزِلَ مِن قَبْلِكَ وَبِٱلْءَاخِرَةِ هُمْ يُوقِنُونَ"
    @test normalize(verses(tnzldata[1][5])[1], :alif_hamza_below) === "اِيَّاكَ نَعْبُدُ وَاِيَّاكَ نَسْتَعِينُ"
    @test normalize(verses(tnzldata[2][3])[1], :waw_hamza_above) === "ٱلَّذِينَ يُوْمِنُونَ بِٱلْغَيْبِ وَيُقِيمُونَ ٱلصَّلَوٰةَ وَمِمَّا رَزَقْنَٰهُمْ يُنفِقُونَ"
    @test normalize(verses(tnzldata[2][15])[1], :ya_hamza_above) === "ٱللَّهُ يَسْتَهْزِيُ بِهِمْ وَيَمُدُّهُمْ فِى طُغْيَٰنِهِمْ يَعْمَهُونَ"
    @test normalize(verses(tnzldata[2][2])[1], :alif_maksura) === "ذَٰلِكَ ٱلْكِتَٰبُ لَا رَيْبَ فِيهِ هُدًي لِّلْمُتَّقِينَ"
    @test normalize(verses(tnzldata[2][3])[1], :ta_marbuta) === "ٱلَّذِينَ يُؤْمِنُونَ بِٱلْغَيْبِ وَيُقِيمُونَ ٱلصَّلَوٰهَ وَمِمَّا رَزَقْنَٰهُمْ يُنفِقُونَ"
    @test normalize(verses(tnzldata[1][1])[1], [:alif_khanjareeya, :hamzat_wasl]) === "بِسْمِ اللَّهِ الرَّحْمَانِ الرَّحِيمِ"

    # setting new transliterator
    my_encoder = Dict(
        Symbol(Char(0x0621)) => Symbol('('),
        Symbol(Char(0x0622)) => Symbol('\''),
        Symbol(Char(0x0623)) => Symbol('&'),
        Symbol(Char(0x0624)) => Symbol('>'),
        Symbol(Char(0x0625)) => Symbol('}'),
        Symbol(Char(0x0626)) => Symbol('<'),
        Symbol(Char(0x0627)) => Symbol('b'),
        Symbol(Char(0x0628)) => Symbol('A'),
        Symbol(Char(0x0629)) => Symbol('t'),
        Symbol(Char(0x062A)) => Symbol('p'),
        Symbol(Char(0x062B)) => Symbol('j'),
        Symbol(Char(0x062C)) => Symbol('v'),
        Symbol(Char(0x062D)) => Symbol('x'),
        Symbol(Char(0x062E)) => Symbol('H'),
        Symbol(Char(0x062F)) => Symbol('*'),
        Symbol(Char(0x0630)) => Symbol('d'),
        Symbol(Char(0x0631)) => Symbol('z'),
        Symbol(Char(0x0632)) => Symbol('r'),
        Symbol(Char(0x0633)) => Symbol('$'),
        Symbol(Char(0x0634)) => Symbol('s'),
        Symbol(Char(0x0635)) => Symbol('D'),
        Symbol(Char(0x0636)) => Symbol('S'),
        Symbol(Char(0x0637)) => Symbol('Z'),
        Symbol(Char(0x0638)) => Symbol('T'),
        Symbol(Char(0x0639)) => Symbol('g'),
        Symbol(Char(0x063A)) => Symbol('E'),
        Symbol(Char(0x0640)) => Symbol('f'),
        Symbol(Char(0x0641)) => Symbol('_'),
        Symbol(Char(0x0642)) => Symbol('k'),
        Symbol(Char(0x0643)) => Symbol('q'),
        Symbol(Char(0x0644)) => Symbol('m'),
        Symbol(Char(0x0645)) => Symbol('l'),
        Symbol(Char(0x0646)) => Symbol('h'),
        Symbol(Char(0x0647)) => Symbol('n'),
        Symbol(Char(0x0648)) => Symbol('Y'),
        Symbol(Char(0x0649)) => Symbol('w'),
        Symbol(Char(0x064A)) => Symbol('F'),
        Symbol(Char(0x064B)) => Symbol('y'),
        Symbol(Char(0x064C)) => Symbol('K'),
        Symbol(Char(0x064D)) => Symbol('N'),
        Symbol(Char(0x064E)) => Symbol('u'),
        Symbol(Char(0x064F)) => Symbol('a'),
        Symbol(Char(0x0650)) => Symbol('~'),
        Symbol(Char(0x0651)) => Symbol('i'),
        Symbol(Char(0x0652)) => Symbol('^'),
        Symbol(Char(0x0653)) => Symbol('o'),
        Symbol(Char(0x0654)) => Symbol('`'),
        Symbol(Char(0x0670)) => Symbol('#'),
        Symbol(Char(0x0671)) => Symbol(':'),
        Symbol(Char(0x06DC)) => Symbol('{'),
        Symbol(Char(0x06DF)) => Symbol('\"'),
        Symbol(Char(0x06E0)) => Symbol('@'),
        Symbol(Char(0x06E2)) => Symbol(';'),
        Symbol(Char(0x06E3)) => Symbol('['),
        Symbol(Char(0x06E5)) => Symbol('.'),
        Symbol(Char(0x06E6)) => Symbol(','),
        Symbol(Char(0x06E8)) => Symbol('-'),
        Symbol(Char(0x06EA)) => Symbol('!'),
        Symbol(Char(0x06EB)) => Symbol('%'),
        Symbol(Char(0x06EC)) => Symbol('+'),
        Symbol(Char(0x06ED)) => Symbol(']')
    );

    basmala = arabic(verses(crpsdata[1][1])[1])

    @transliterator my_encoder "MyEncoder"
    @test encode(basmala) === "A~\$^l~ :mmiun~ :mziux^lu#h~ :mziux~Fl~"
    @test arabic(encode(basmala)) === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ"
    @test dediac(encode(basmala)) === "A\$l :mmn :mzxl#h :mzxFl"
    @test normalize(encode(basmala)) === "A~\$^l~ bmmiun~ bmziux^lubh~ bmziux~Fl~"

    @transliterator :default
    @test encode(basmala) === "bisomi {ll~ahi {lr~aHoma`ni {lr~aHiymi"
    @test arabic(encode(basmala)) === "بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ"
    @test dediac(encode(basmala)) === "bsm {llh {lrHm`n {lrHym"
    @test normalize(encode(basmala)) === "bisomi All~ahi Alr~aHomaAni Alr~aHiymi"

    # features
    @test root(parse(Features, select(crpsdata[112].data, :features)[1])) === "qwl"
    @test lemma(parse(Features, select(crpsdata[112].data, :features)[1])) === "qaAla"
    @test special(parse(Features, select(crpsdata.data, :features)[53])) === "<in~"
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[2]), Stem) === true
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[end-4]), Suffix) === true
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[end-3]), Prefix) === true
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[2]), Noun) === true
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[3]), ProperNoun) === true
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[5]), Adjective) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[68054]), ImperativeVerbalNoun) === true
    @test isfeature(parse(Features, select(crpsdata[1].data, :features)[23]), Personal) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[50]), Demonstrative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[35]), Relative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[210]), Time) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[291]), Location) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[14]), Plural) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[38]), Preposition) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[419]), EmphaticLam) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[5717]), ImperativeLam) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[2131]), PurposeLam) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[997]), EmphaticNun) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[25]), Coordinating) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[245]), Subordinating) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[31]), Accusative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[233]), Amendment) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[2252]), Answer) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[66305]), Aversion) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[926]), Cause) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[1544]), Certainty) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[166]), Circumstantial) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[14139]), Comitative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[416]), Conditional) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[120]), Equalization) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[3532]), Exhortation) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[618]), Explanation) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[821]), Exceptive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[17319]), Future) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[226]), Inceptive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[3758]), Interpretation) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[640]), Interogative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[46]), Negative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[223]), Preventive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[214]), Prohibition) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[72]), Resumption) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[182]), Restriction) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[2542]), Retraction) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[504]), Result) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[232]), Supplemental) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[18810]), Surprise) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[10975]), Vocative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[49]), DisconnectedLetters) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[24]), FirstPerson) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[23]), SecondPerson) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[39]), ThirdPerson) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[2]), Masculine) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[71]), Feminine) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[5]), Singular) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[2303]), Dual) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[14]), Plural) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[24]), Verb) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[36]), Perfect) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[24]), Imperfect) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[28]), Imperative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[531]), Subjunctive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[126]), Jussive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[86]), Passive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[499]), VerbFormII) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[171]), VerbFormIII) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[36]), VerbFormIV) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[970]), VerbFormV) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[591]), VerbFormVI) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[1538]), VerbFormVII) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[60]), VerbFormVIII) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[13216]), VerbFormIX) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[27]), VerbFormX) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[114155]), VerbFormXI) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[9063]), VerbFormXII) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[19]), ActiveParticle) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[42]), PassiveParticle) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[305]), VerbalNoun) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[57]), Indefinite) === true # no verse with DEF feature
    @test isfeature(parse(Features, select(crpsdata.data, :features)[9]), Nominative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[2]), Genetive) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[24]), VerbFormI) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[27]), VerbFormI) === false
    @test isfeature(parse(Features, select(crpsdata.data, :features)[27]), Indicative) === true
    @test isfeature(parse(Features, select(crpsdata.data, :features)[126]), Indicative) === false
    @test isfeature(parse(Features, select(crpsdata.data, :features)[86]), Active) === false
    @test isfeature(parse(Features, select(crpsdata.data, :features)[36]), Active) === true

    @test encode(SimpleEncoder, basmala) === "Ba+Kasra | Seen+Sukun | Meem+Kasra | <space> | HamzatWasl | Lam | Lam+Shadda+Fatha | Ha+Kasra | <space> | HamzatWasl | Lam | Ra+Shadda+Fatha | HHa+Sukun | Meem+Fatha | AlifKhanjareeya | Noon+Kasra | <space> | HamzatWasl | Lam | Ra+Shadda+Fatha | HHa+Kasra | Ya | Meem+Kasra"

    # printing
    out = capture_io(crpsdata)
    @test out === "Quranic Arabic Corpus (morphology)\n(C) 2011 Kais Dukes\n\nTable with 128219 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(tnzldata)
    @test out === "Tanzil Quran Text (Uthmani)\n(C) 2008-2010 Tanzil.net\n\nTable with 6236 rows, 3 columns:\nchapter  verse  form\n─────────────────────────────────────────────────────────────────────\n1        1      \"بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ\"\n1        2      \"ٱلْحَمْدُ لِلَّهِ رَبِّ ٱلْعَٰلَمِينَ\"\n1        3      \"ٱلرَّحْمَٰنِ ٱلرَّحِيمِ\"\n1        4      \"مَٰلِكِ يَوْمِ ٱلدِّينِ\"\n1        5      \"إِيَّاكَ نَعْبُدُ وَإِيَّاكَ نَسْتَعِينُ\"\n1        6      \"ٱهْدِنَا ٱلصِّرَٰطَ ٱلْمُسْتَقِيمَ\"\n1        7      \"صِرَٰطَ ٱلَّذِينَ أَنْعَمْتَ عَلَيْهِمْ غَيْرِ ٱلْمَغْضُوبِ عَلَيْهِمْ وَلَا ٱلضَّآلِّينَ\"\n2        1      \"بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ الٓمٓ\"\n2        2      \"ذَٰلِكَ ٱلْكِتَٰبُ لَا رَيْبَ فِيهِ هُدًى لِّلْمُتَّقِينَ\"\n⋮\n113      4      \"وَمِن شَرِّ ٱلنَّفَّٰثَٰتِ فِى ٱلْعُقَدِ\"\n113      5      \"وَمِن شَرِّ حَاسِدٍ إِذَا حَسَدَ\"\n114      1      \"بِسْمِ ٱللَّهِ ٱلرَّحْمَٰنِ ٱلرَّحِيمِ قُلْ أَعُوذُ بِرَبِّ ٱلنَّاسِ\"\n114      2      \"مَلِكِ ٱلنَّاسِ\"\n114      3      \"إِلَٰهِ ٱلنَّاسِ\"\n114      4      \"مِن شَرِّ ٱلْوَسْوَاسِ ٱلْخَنَّاسِ\"\n114      5      \"ٱلَّذِى يُوَسْوِسُ فِى صُدُورِ ٱلنَّاسِ\"\n114      6      \"مِنَ ٱلْجِنَّةِ وَٱلنَّاسِ\"\n\n"

    out = capture_io(Noun())
    @test out === "N\n"

    out = capture_io(crps);
    @test out[5000:6000] === "ka\\tDEM\\tSTEM|POS:DEM|LEM:*a`lik|MS\", \"(2:2:2:1)\\t{lo\\tDET\\tPREFIX|Al+\", \"(2:2:2:2)\\tkita`bu\\tN\\tSTEM|POS:N|LEM:kita`b|ROOT:ktb|M|NOM\", \"(2:2:3:1)\\tlaA\\tNEG\\tSTEM|POS:NEG|LEM:laA|SP:<in~\", \"(2:2:4:1)\\trayoba\\tN\\tSTEM|POS:N|LEM:rayob|ROOT:ryb|M|ACC\", \"(2:2:5:1)\\tfiy\\tP\\tSTEM|POS:P|LEM:fiY\", \"(2:2:5:2)\\thi\\tPRON\\tSUFFIX|PRON:3MS\", \"(2:2:6:1)\\thudFY\\tN\\tSTEM|POS:N|LEM:hudFY|ROOT:hdy|M|INDEF|NOM\", \"(2:2:7:1)\\tl~i\\tP\\tPREFIX|l:P+\", \"(2:2:7:2)\\tlo\\tDET\\tPREFIX|Al+\", \"(2:2:7:3)\\tmut~aqiyna\\tN\\tSTEM|POS:N|ACT|PCPL|(VIII)|LEM:mut~aqiyn|ROOT:wqy|MP|GEN\", \"(2:3:1:1)\\t{l~a*iyna\\tREL\\tSTEM|POS:REL|LEM:{l~a*iY|MP\", \"(2:3:2:1)\\tyu&ominu\\tV\\tSTEM|POS:V|IMPF|(IV)|LEM:'aAmana|ROOT:Amn|3MP\", \"(2:3:2:2)\\twna\\tPRON\\tSUFFIX|PRON:3MP\", \"(2:3:3:1)\\tbi\\tP\\tPREFIX|bi+\", \"(2:3:3:2)\\t{lo\\tDET\\tPREFIX|Al+\", \"(2:3:3:3)\\tgayobi\\tN\\tSTEM|POS:N|LEM:gayob|ROOT:gyb|M|GEN\", \"(2:3:4:1)\\twa\\tCONJ\\tPREFIX|w:CONJ+\", \"(2:3:4:2)\\tyuqiymu\\tV\\tSTEM|POS:V|IMPF|(IV)|LEM:>aqaAma|ROOT:qwm|3MP\", \"(2:3:4:3)\\twna\\tPRON\\tSUFFIX|PRON:3MP\""

    out = capture_io(tnzl);
    @test out[5000:5110] === "ن قَبْلِكُمْ لَعَلَّكُمْ تَتَّقُونَ\", \"2|22|ٱلَّذِى جَعَلَ لَكُ"

    out = capture_io(crpsdata[1])
    @test out === "Chapter 1: ٱلْفَاتِحَة (The Opening)\n\nTable with 48 rows, 6 columns:\nColumns:\n#  colname   type\n───────────────────\n1  verse     Int64\n2  word      Int64\n3  part      Int64\n4  form      String\n5  tag       String\n6  features  String\n\n"

    out = capture_io(crpsdata[[112,113]])
    @test out === "Chapters: \n ├112 (ٱلْإِخْلَاص-Purity of Faith) \n └113 (ٱلْفَلَق-Daybreak)\n\nTable with 49 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[end-1:end])
    @test out === "Chapter 113-114: ٱلْفَلَق-ٱلنَّاس (Daybreak-People)\n\nTable with 60 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[1][1])
    @test out === "Chapter 1 ٱلْفَاتِحَة (The Opening)\nVerse 1\n\nTable with 7 rows, 5 columns:\nword  part  form          tag    features\n─────────────────────────────────────────────────────────────────────────────\n1     1     \"bi\"          \"P\"    \"PREFIX|bi+\"\n1     2     \"somi\"        \"N\"    \"STEM|POS:N|LEM:{som|ROOT:smw|M|GEN\"\n2     1     \"{ll~ahi\"     \"PN\"   \"STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|GEN\"\n3     1     \"{l\"          \"DET\"  \"PREFIX|Al+\"\n3     2     \"r~aHoma`ni\"  \"ADJ\"  \"STEM|POS:ADJ|LEM:r~aHoma`n|ROOT:rHm|MS|GEN\"\n4     1     \"{l\"          \"DET\"  \"PREFIX|Al+\"\n4     2     \"r~aHiymi\"    \"ADJ\"  \"STEM|POS:ADJ|LEM:r~aHiym|ROOT:rHm|MS|GEN\"\n\n"

    out = capture_io(crpsdata[1][1:2])
    @test out === "Chapter 1 ٱلْفَاتِحَة (The Opening)\nVerses 1-2\n\nTable with 14 rows, 6 columns:\nColumns:\n#  colname   type\n───────────────────\n1  verse     Int64\n2  word      Int64\n3  part      Int64\n4  form      String\n5  tag       String\n6  features  String\n\n"

    out = capture_io(crpsdata[1][[1,2]])
    @test out === "Chapter 1 ٱلْفَاتِحَة (The Opening)\nVerses 1, 2\n\nTable with 14 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[end-1:end][1])
    @test out === "Chapters 113-114: ٱلْفَلَق-ٱلنَّاس (Daybreak-People)\nVerse 1\n\nTable with 12 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[end-1:end][[1,2]])
    @test out === "Chapters 113-114: ٱلْفَلَق-ٱلنَّاس (Daybreak-People)\nVerses 1, 2\n\nTable with 19 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[end-1:end][1:2])
    @test out === "Chapters 113-114: ٱلْفَلَق-ٱلنَّاس (Daybreak-People)\nVerses 1-2\n\nTable with 19 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[[112,113]][1])
    @test out === "Chapters: \n ├112 (ٱلْإِخْلَاص-Purity of Faith) \n └113 (ٱلْفَلَق-Daybreak)\nVerse 1\n\nTable with 10 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = capture_io(crpsdata[[112,113]][[1,2]])
    @test out === "Chapters: \n ├112 (ٱلْإِخْلَاص-Purity of Faith) \n └113 (ٱلْفَلَق-Daybreak)\nVerses 1, 2\n\nTable with 17 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    meta = MetaData(
        "Quranic Arabic Corpus (morphology)",
        "Kais Dukes",
        "The Quranic Arabic Corpus includes syntactic and morphological\nannotation of the Qur'an, and builds on the verified Arabic text\ndistributed by the Tanzil project.",
        "http://corpus.quran.com/",
        "English",
        "2011",
        "GNU General Public License",
        "0.4"
    )
    out = capture_io(meta)
    @test out === "Quranic Arabic Corpus (morphology) v0.4\nCopyright (C) 2011 Kais Dukes\nGNU General Public License\nhttp://corpus.quran.com/\n\nThe Quranic Arabic Corpus includes syntactic and morphological\nannotation of the Qur'an, and builds on the verified Arabic text\ndistributed by the Tanzil project.\n\n"

    out = capture_io(crpsdata[[112,113]][1:2])
    @test out === "Chapters: \n ├112 (ٱلْإِخْلَاص-Purity of Faith) \n └113 (ٱلْفَلَق-Daybreak)\nVerses 1-2\n\nTable with 17 rows, 7 columns:\nColumns:\n#  colname   type\n───────────────────\n1  chapter   Int64\n2  verse     Int64\n3  word      Int64\n4  part      Int64\n5  form      String\n6  tag       String\n7  features  String\n\n"

    out = @capture_out begin
        description(parse(Features, select(crpsdata[1].data, :features)[1]))
    end;
    @test out === """Prefix
    ──────
    Preposition:
     ├ data: P
     ├ desc: Preposition prefix ('by', 'with', 'in')
     └ ar_label: حرف جر\n"""

    out = @capture_out begin
        @desc parse(Features, select(crpsdata[1].data, :features)[end])
    end;
    @test out === """Stem
    ────
    Noun:
     ├ data: N
     ├ desc: Noun
     └ ar_label: اسم
    QuranTree.Lemma:
     └ data: DaA^l~
    QuranTree.Root:
     └ data: Dll
    ActiveParticle:
     ├ data: ACT PCPL
     ├ desc: Active particle
     └ ar_label: اسم فاعل
    Active:
     ├ data: ACT
     ├ desc: Active voice (default)
     └ ar_label: مبني للمعلوم
    Masculine:
     ├ data: M
     ├ desc: Masculine
     └ ar_label: الجنس
    Plural:
     ├ data: P
     ├ desc: Plural
     └ ar_label: العدد
    Genetive:
     ├ data: GEN
     ├ desc: Genetive case
     └ ar_label: مجرور\n"""

    out = @capture_out begin crpsdata[1][1][1][1] end;
    out === """Chapter 1 ٱلْفَاتِحَة (The Opening)
    Verse 1

    Table with 1 rows, 5 columns:
    word  part  form  tag  features
    ─────────────────────────────────────────────
    1     1     "bi"  "P"  Features("PREFIX|bi+")\n""";

    out = @desc(1)
    @test out === missing;

    # out = @capture_out(pretty_table(crpsdata));
    # @test out[1000:2000] === "────────────────────────────────────────────────────────\n        1       1       1       1                 bi        P                                                             Features(\"PREFIX|bi+\")\n        1       1       1       2               somi        N                                     Features(\"STEM|POS:N|LEM:{som|ROOT:smw|M|GEN\")\n        1       1       2       1            {ll~ahi       PN                                    Features(\"STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|GEN\")\n        1       1       3       1                 {l      DET                                                             Features(\"PREFIX|Al+\")\n        1       1       3       2         r~aHoma`ni      ADJ                             Features(\"STEM|POS:ADJ|LEM:r~aHoma`n|ROOT:rHm|MS|GEN\")\n        1       1       4       1                 {l      DET                                              "

    # out = @capture_out(pretty_table(crpsdata[1]));
    # @test out[1000:2000] === "                                                Features(\"PREFIX|bi+\")\n      1       1       2          somi        N                         Features(\"STEM|POS:N|LEM:{som|ROOT:smw|M|GEN\")\n      1       2       1       {ll~ahi       PN                        Features(\"STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|GEN\")\n      1       3       1            {l      DET                                                 Features(\"PREFIX|Al+\")\n      1       3       2    r~aHoma`ni      ADJ                 Features(\"STEM|POS:ADJ|LEM:r~aHoma`n|ROOT:rHm|MS|GEN\")\n      1       4       1            {l      DET                                                 Features(\"PREFIX|Al+\")\n      1       4       2      r~aHiymi      ADJ                   Features(\"STEM|POS:ADJ|LEM:r~aHiym|ROOT:rHm|MS|GEN\")\n      2       1       1           {lo      DET                                                 Features(\"PREFIX|Al+\")\n      2       1       2        Hamodu        N                        Features(\"STEM|POS:N|LEM:Hamod|ROO"

    # out = @capture_out(pretty_table(crpsdata[1][1]))
    # @test out === "────────────────────────────────────────────────────────────────────────────────────────────────\n   word    part         form      tag                                                 features\n  Int64   Int64       String   String                                                 Features\n────────────────────────────────────────────────────────────────────────────────────────────────\n      1       1           bi        P                                   Features(\"PREFIX|bi+\")\n      1       2         somi        N           Features(\"STEM|POS:N|LEM:{som|ROOT:smw|M|GEN\")\n      2       1      {ll~ahi       PN          Features(\"STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|GEN\")\n      3       1           {l      DET                                   Features(\"PREFIX|Al+\")\n      3       2   r~aHoma`ni      ADJ   Features(\"STEM|POS:ADJ|LEM:r~aHoma`n|ROOT:rHm|MS|GEN\")\n      4       1           {l      DET                                   Features(\"PREFIX|Al+\")\n      4       2     r~aHiymi      ADJ     Features(\"STEM|POS:ADJ|LEM:r~aHiym|ROOT:rHm|MS|GEN\")\n────────────────────────────────────────────────────────────────────────────────────────────────\n"

    # out = @capture_out(pretty_table(crpsdata[1][1:2]))
    # @test out[end-1000:end] === "{lo      DET                                   Features(\"PREFIX|Al+\")\n      2       1       2       Hamodu        N          Features(\"STEM|POS:N|LEM:Hamod|ROOT:Hmd|M|NOM\")\n      2       2       1           li        P                                  Features(\"PREFIX|l:P+\")\n      2       2       2        l~ahi       PN          Features(\"STEM|POS:PN|LEM:{ll~ah|ROOT:Alh|GEN\")\n      2       3       1        rab~i        N           Features(\"STEM|POS:N|LEM:rab~|ROOT:rbb|M|GEN\")\n      2       4       1          {lo      DET                                   Features(\"PREFIX|Al+\")\n      2       4       2   Ea`lamiyna        N     Features(\"STEM|POS:N|LEM:Ea`lamiyn|ROOT:Elm|MP|GEN\")\n────────────────────────────────────────────────────────────────────────────────────────────────────────\n"

    # remaining todo:
    #   check pretty_table
    #   test display/show
    #   use suppressor

end


text:
# Scripts for field line postprocessing.
#
# Create seeds for field line tracing
# Identify boundary points from Paraview polyline tracing output
# Integrate electric fields along the magnetopause boundary
#
# Hongyang Zhou, hyzhou@umich.edu 11/11/2019

using DelimitedFiles, PyPlot, NumericalIntegration, LinearAlgebra
using LazySets
#using MATLAB
#using ConcaveHull

const Rg = 2634000        # [m], radius of Ganymede
const e  = 1.60217662e-19 # [C], electron charge

function findBoundary(filename, DoPlot::Bool=false)

   fieldline = readdlm(filename, ',', header=true)
   data = fieldline[1]

   # count the total number of field lines
   ilinetotal = 0
   start = []
   for i=1:size(data,1)
      if data[i,4] == 0.0
         ilinetotal += 1
         push!(start, i)
      end
   end

   line = Array{Float64, 2}(undef, 6, ilinetotal)
   iline = 1

   for i=1:ilinetotal-1
      iStart, iEnd = start[i], start[i+1]-1
      line[1,i] = data[iStart,5]
      line[2,i] = data[iStart,6]
      line[3,i] = data[iStart,7]
      line[4,i] = data[iEnd,5]
      line[5,i] = data[iEnd,6]
      line[6,i] = data[iEnd,7]
   end

   line[1,end] = data[start[end],5]
   line[2,end] = data[start[end],6]
   line[3,end] = data[start[end],7]
   line[4,end] = data[end,5]
   line[5,end] = data[end,6]
   line[6,end] = data[end,7]

   # status
   # (0: open, 1: closed along B, 2: closed along -B 3: fully closed)
   # (-1: cells inside body, -2: loop ray within block, -3: strange)
   status = fill(0, ilinetotal)

   # Cluster seed points based on their end points
   for i=1:ilinetotal
      if line[4,i]^2 + line[5,i]^2 + line[6,i]^2 < 1.1^2
         status[i] = 1
         #elseif abs(line[4,i]) ≥ 2.9 || abs(line[5,i]) ≥ 2.9 || abs(line[6,i]) ≥ 2.9
         #   status[i] = 0
      end
   end

   pts_closed = line[1:3, status .== 1]
   pts_open = line[1:3, status .== 0]

   if DoPlot
      scatter3D(pts_closed[1,:], pts_closed[2,:], pts_closed[3,:])
      scatter3D(pts_open[1,:], pts_open[2,:], pts_open[3,:])
   end

   x = pts_closed[1,:]
   y = pts_closed[2,:]
   z = pts_closed[3,:]

   # Matlab
   #=
   k, area = mxcall(:boundary, 2, x, y)
   k = Int.(k)
   figure()
   plot(x[k], y[k])
   return x[k], y[k], z[k]
   =#

   # LazySets
   v = [[x[i], y[i]] for i in 1:length(x)]
   hull = convex_hull(v) # LazySets
   xy = hcat(hull...)
   x = [xy[1,:]; xy[1,1]]
   y = [xy[2,:]; xy[2,1]]
   z = fill(2.0, length(x))
   return x, y, z

   # ConcaveHull
   #=
   v = [[x[i], y[i]] for i in 1:length(x)]
   hull = concave_hull(v,200) # ConcaveHull
   x = [p[1] for p in hull.vertices]
   y = [p[2] for p in hull.vertices]
   z = fill(2.0, length(x))
   return x,y,z
   =#
end

function saveBoundary(x,y,z, filename::String="boundary.txt")
   open(filename, "w") do io
      write(io, "\"X\",\"Y\",\"Z\"\n")
      writedlm(io, [x y z], ',')
   end
end

function createSeeds(filename::String="seeds.txt", dir::String="./", Verbose::Bool=false)

   xMin, xMax = -2.0, 4.3
   yMin, yMax = -3.0, 3.0
   z = 2.0

   dx, dy = 5e-2, 5e-2

   xRange = xMin:dx:xMax
   yRange = yMin:dy:yMax

   X = [x for x in xRange, y in yRange if (x-1)^2 + y^2 > 2.4^2 && (x-1)^2 + y^2 < 3.5^2]
   Y = [y for x in xRange, y in yRange if (x-1)^2 + y^2 > 2.4^2 && (x-1)^2 + y^2 < 3.5^2]
   Z = fill(z, size(X))

   open(dir*filename, "w") do io
      write(io, "\"X\",\"Y\",\"Z\"\n")
      writedlm(io, [X Y Z], ',')
   end

   Verbose && println("Seeds created.")
end

"""
    integrate_along_boundary(filename)

Integrate variables along the 2D boundary line.

May add Hall velocity later!
"""
function integrate_along_boundary(filename, DoPlot::Bool=false)
   bc = readdlm(filename, ',', header=true)
   data = bc[1]

   x = @view data[:,7]
   y = @view data[:,8]

   # Get E field: E = -uxB
   Ex = Vector{Float64}(undef,size(data,1))
   Ey = Vector{Float64}(undef,size(data,1))
   Ez = Vector{Float64}(undef,size(data,1))

   for i = 1:size(data,1)
      Ex[i], Ey[i], Ez[i] = -cross(data[i,1:3], data[i,4:6])
   end

   ϕx = cumul_integrate(x, Ex)
   ϕy = cumul_integrate(y, Ey)
   Coef = 1e3 * 1e-9 * Rg * 1e-3 # [km/s]*[nT]*[Rg] --> [kV]
   ϕ  = (ϕx .+ ϕy) .* Coef # [kV]

   if DoPlot
      figure()
      plot(ϕ)
   end

   maximum(ϕ) - minimum(ϕ)
end

function integrate_along_boundary_hall(filename, DoPlot::Bool=false)
   bc = readdlm(filename, ',', header=true)
   data = bc[1]

   ρ = @view data[:,1]
   ux= @view data[:,2]
   uy= @view data[:,3]
   uz= @view data[:,4]
   bx= @view data[:,5]
   by= @view data[:,6]
   bz= @view data[:,7]
   jx= @view data[:,8]
   jy= @view data[:,9]
   jz= @view data[:,10]
   x = @view data[:,11]
   y = @view data[:,12]

   # Get E field: E = -uxB
   Ex = Vector{Float64}(undef,size(data,1))
   Ey = Vector{Float64}(undef,size(data,1))
   Ez = Vector{Float64}(undef,size(data,1))

   # Calculate electron bulk velocity in Hall MHD
   # Transform the hall velocity into [km/s]
   Coef = 1e-6 /e * 1e-3 * 1e-6
   uex = @. ux - jx/ρ * Coef
   uey = @. uy - jy/ρ * Coef
   uez = @. uz - jz/ρ * Coef

   for i = 1:size(data,1)
      Ex[i], Ey[i], Ez[i] =
         -cross([uex[i], uey[i], uez[i]], [bx[i], by[i], bz[i]])
   end

   ϕx = cumul_integrate(x, Ex)
   ϕy = cumul_integrate(y, Ey)
   Coef = 1e3 * 1e-9 * Rg * 1e-3 # [km/s]*[nT]*[Rg] --> [kV]
   ϕ  = (ϕx .+ ϕy) .* Coef # [kV]

   if DoPlot
      figure()
      plot(ϕ)
   end

   #maximum(ϕ) - minimum(ϕ) # This would return the CPCP on the night side
   maximum(ϕ) - minimum(ϕ) + ϕ[end] # Return the CPCP on the day side
end