text:
module SimradRaw

export datagrams, readencapsulateddatagram, readdatagramblock,
    readdatagrambody, readdatagram, readdatagramheader, Datagram,
    filetime

# References
# Simrad EK60 Context sensitive on-line help, Release 2.4.X,
# http://www.simrad.net/ek60_ref_english/default.htm

# Read low level C types

function readlong(stream::IO)
    read(stream, Int32)
end

function readshort(stream::IO)
    read(stream, Int16)
end

function readfloat(stream::IO)
    reinterpret(Float32, read(stream, 4))[1]
end

function readbytes(stream::IO, n::Integer)
    read(stream, n)
end

function readchars(stream::IO, n::Integer)
    String(read(stream, n))
end

function readstring(stream::IO, n::Integer)
    rstrip(readchars(stream, n), '\0')
end

function readfloats(stream::IO, n::Integer)
    reinterpret(Float32, read(stream, 4 * n))
end

function readshorts(stream::IO, n::Integer)
    b = Array{Int16}(undef,n)
    read!(stream, b)
    return b
end

function readint8s(stream::IO, n::Integer)
    b = Array{Int8}(undef,n)
    read!(stream, b)
    return b
end

function readushorts(stream::IO, n::Integer)
    read(stream, UInt16, n)
end

function readdword(stream::IO)
    read(stream, UInt32)
end

#

struct DateTime
    lowdatetime::UInt32
    highdatetime::UInt32
end

function readdatetime(stream::IO)
    lowdatetime = readdword(stream)
    highdatetime = readdword(stream)
    DateTime(lowdatetime, highdatetime)
end

"
Returns the filetime being the number of 100 nanosecond intervals
since January 1, 1601 of a given DateTime.
"
function filetime(d::DateTime)
    d.highdatetime * 4294967296 + d.lowdatetime
end

#

"""
    datagrams(filename, datagramreader=readdatagram)

`datagrams` returns an iterator over the datagrams in the RAW file
designated by `filename`.

"""
function datagrams(filename::AbstractString;
                   datagramreader=readdatagram::Function, handleexceptions=false)

    datagrams([filename], datagramreader=datagramreader, handleexceptions=handleexceptions)
end

function datagrams(filenames::Vector{String};
                   datagramreader=readdatagram::Function,
                   handleexceptions=false)
    function _it(chn1)
        for filename in filenames

            # FIXME: The spec requires us to check endianness by comparing
            # length fields, but in practice, everyone is using PC/Windows

            open(filename) do f
                er = false
                while !eof(f) && !er
                    try
                        datagram = readencapsulateddatagram(f, datagramreader=datagramreader)
                        put!(chn1, datagram)
                    catch ex
                        if !handleexceptions
                            rethrow(ex)
                        end
                        # We regularly see corruption towards the end of RAW files
                        p = position(f)
                        @warn("Can't read datagram from $filename ($p), skipping")
                        er = true
                    end
                end
            end
        end
    end

    return Channel(_it, ctype=Datagram)
end

#

function load(filename::AbstractString)
    collect(datagrams(filename))
end

#

function readencapsulateddatagram(stream::IO;
                                  datagramreader=readdatagram::Function)
    length = readlong(stream)
    datagram = datagramreader(stream, length)
    length2 = length = readlong(stream)

    if length != length2
        error("Invalid datagram")
    end
    datagram
end

#

abstract type Datagram end

struct DatagramHeader
    datagramtype::String
    datetime::DateTime
end

"""
    filetime(d::DatagramHeader)

Returns the filetime being the number of 100 nanosecond intervals
since January 1, 1601 of a given DatagramHeader.
"""
function filetime(d::DatagramHeader)
    filetime(d.datetime)
end

function readdatagramheader(stream::IO, length::Integer)
    headerlength = 12
    datagramtype = readchars(stream, 4)
    datetime = readdatetime(stream)
    l = length - headerlength
    dgheader = DatagramHeader(datagramtype, datetime)
    body = readdatagramblock(stream, l)

    return dgheader, body
end

function readdatagrambody(stream::IO, length::Integer, dgheader::DatagramHeader)

    datagramtype = dgheader.datagramtype

    if datagramtype == "XML0"
        return readxmldatagram(stream, length, dgheader)
    elseif datagramtype == "FIL1"
        return readbinarydatagram(stream, length, dgheader)
    elseif datagramtype == "CON0"
        return readconfigurationdatagram(stream, dgheader)
    elseif datagramtype == "NME0"
        return readtextdatagram(stream, length, dgheader)
    elseif datagramtype == "RAW0"
        return readsamplebinarydatagram0(stream, dgheader)
    elseif datagramtype == "RAW3"
        return readsamplebinarydatagram3(stream, dgheader)
    elseif datagramtype == "MRU0"
        return readmrudatagram(stream, dgheader)
    elseif datagramtype == "TAG0"
        return readtextdatagram(stream, length, dgheader)
    else
        d = readbinarydatagram(stream, length, dgheader)
        @warn("Unrecognised datagram type, treating as binary, ", datagramtype)
        return d
    end
end

function readdatagram(stream::IO, length::Integer)

    headerlength = 12
    datagramtype = readchars(stream, 4)
    datetime = readdatetime(stream)

    dgheader = DatagramHeader(datagramtype, datetime)

    l = length - headerlength

    return readdatagrambody(stream, l, dgheader)

end

#

function readdatagramblock(stream::IO, length::Integer)
    read(stream, length)
end

#

struct XMLDatagram <: Datagram
    dgheader::DatagramHeader
    text::String
end

function readxmldatagram(stream::IO, l::Integer, dgheader::DatagramHeader)
    text = readchars(stream, l)
    XMLDatagram(dgheader, text)
end

#

struct BinaryDatagram <: Datagram
    dgheader::DatagramHeader
    bytes
end

function readbinarydatagram(stream::IO, l::Integer, dgheader::DatagramHeader)
    bytes = read(stream, l)
    BinaryDatagram(dgheader, bytes)
end

# The MRU binary datagram contains motion sensor data

struct MRUDatagram <: Datagram
    dgheader::DatagramHeader
    heave::Float32
    roll::Float32
    pitch::Float32
    heading::Float32
end

function readmrudatagram(stream::IO, dgheader::DatagramHeader)
    heave =  readfloat(stream)
    roll = readfloat(stream)
    pitch = readfloat(stream)
    heading = readfloat(stream)

    MRUDatagram(dgheader, heave, roll, pitch, heading)
end

#

struct TextDatagram <: Datagram
    dgheader::DatagramHeader
    text::String
end

function readtextdatagram(stream::IO, length::Integer, dgheader::DatagramHeader)
    text = readchars(stream, length)
    TextDatagram(dgheader, text)
end

# EK80 sample datagram

struct SampleDatagram3 <: Datagram
    dgheader::DatagramHeader
    channelid::String
    datatype::Int16
    offset::Int32
    count::Int32
    samples # TODO
end

# The sample datagram 3 contains sample data from each "ping". The
# datagram may have a different size and contain different kind of
# data, depending on the DataType parameter.

function readsamplebinarydatagram3(stream::IO, dgheader::DatagramHeader)
    channelid = readstring(stream, 128)
    datatype = readshort(stream)
    spare = readchars(stream, 2)
    offset = readlong(stream);
    count = readlong(stream);

    # The number of values in Samples[] depends on the value of Count
    # and the Datatype.  As an example a DataType decimal value of
    # 1032 means that Samples[] contains ComplexFloat32 samples and
    # that each sample consists of 4 complex numbers (one from each of
    # the 4 transducer quadrants).

    if datatype == 1032
        samples = readbytes(stream, 4 * 2 * count * 4)
    else
        error("Datatype $(datatype) not yet implemented.")
    end

    SampleDatagram3(dgheader, channelid, datatype, offset, count, samples)
end

# EK60 sample datagram

struct SampleDatagram0 <: Datagram
    dgheader::DatagramHeader
    channel::Int16 # Channel number
    mode::Int16 # Datatype
    transducerdepth::Float32 # [m]
    frequency::Float32 # [Hz]
    transmitpower::Float32 # [W]
    pulselength::Float32 # [s]
    bandwidth::Float32 # [Hz]
    sampleinterval::Float32 # [s]
    soundvelocity::Float32 # [m/s]
    absorptioncoefficient::Float32 # [dB/m]
    heave::Float32 # [m]
    txroll::Float32 # [deg]
    txpitch::Float32 # [deg]
    temperature::Float32 # [C]
    rxroll::Float32 # [Deg]
    rxpitch::Float32 # [Deg]
    offset::Int32 #First sample
    count::Int32 # Number of samples
    power::Array{Int16} # Compressed format - See Remark 1!

    # 1. Power: The power data contained in the sample datagram is
    # compressed. In order to restore the correct value(s), you must
    # decompress the value according to the equation below.
    # y = x (10    # * log10(2) / 256) where: x = power value derived
    # from the datagram, y = converted value (in dB)

    athwartshipangle::Array{Int8}
    alongshipangle::Array{Int8} # See Remark 2 below!

    # 2. Angle: The fore-and-aft (alongship) and athwartship
    # electrical angles are output as one 16-bit word. The alongship
    # angle is the most significant byte while the athwartship angle
    # is the least significant byte. Angle data is expressed in 2's
    # complement format, and the resolution is given in steps of
    # 180/128 electrical degrees per unit. Positive numbers denotes
    # the fore and starboard directions.

end

function readsamplebinarydatagram0(stream::IO, dgheader::DatagramHeader)
    channel = readshort(stream) # Channel number
    mode = readshort(stream) # Datatype
    transducerdepth = readfloat(stream) # [m]
    frequency = readfloat(stream) # [Hz]
    transmitpower = readfloat(stream) # [W]
    pulselength = readfloat(stream) # [s]
    bandwidth = readfloat(stream) # [Hz]
    sampleinterval = readfloat(stream) # [s]
    soundvelocity = readfloat(stream) # [m/s]
    absorptioncoefficient = readfloat(stream) # [dB/m]
    heave = readfloat(stream) # [m]
    txroll = readfloat(stream) # [deg]
    txpitch = readfloat(stream) # [deg]
    temperature= readfloat(stream) # [C]
    spare = readshort(stream)
    spare = readshort(stream)
    rxroll = readfloat(stream) # [Deg]
    rxpitch = readfloat(stream) # [Deg]
    offset = readlong(stream) # First sample
    count = readlong(stream) # Number of samples

    power = readshorts(stream, count) # Compressed format - See Remark 1!

    athwartshipangle = []
    alongshipangle = []

    if mode != 1
        angles = readint8s(stream, count * 2)
        athwartshipangle = angles[1:2:end]
        alongshipangle = angles[2:2:end]
    end

    SampleDatagram0(dgheader, channel, mode, transducerdepth,
                    frequency, transmitpower, pulselength, bandwidth,
                    sampleinterval, soundvelocity,
                    absorptioncoefficient, heave, txroll, txpitch,
                    temperature, rxroll, rxpitch, offset, count,
                    power, athwartshipangle, alongshipangle)

end

# EK60 configuration datagram

struct ConfigurationHeader
    surveyname::String
    transectname::String
    soundername::String
    version::String
    transducercount::Int32
end

function readconfigurationheader(stream::IO)

    surveyname = readstring(stream, 128) # "Loch Ness"
    transectname = readstring(stream, 128)
    soundername = readstring(stream, 128) # "ER60"
    version = readstring(stream, 30)
    spare = readchars(stream, 98)
    transducercount = readlong(stream) # 1 to 7

    ConfigurationHeader(surveyname, transectname,
                        soundername, version, transducercount)

end

#

struct ConfigurationTransducer
    channelid::String
    beamtype # 0 = Single, 1 = Split
    frequency # [Hz]
    gain # [dB]
    equivalentbeamangle # [dB]
    beamwidthalongship # [degree]
    beamwidthathwartship # [degree]
    anglesensitibityalongship
    anglesensitivityathwartship
    angleoffsetalongship # [degree]
    angleoffsetathwartship # [degree]
    posx # future use
    posy # future use
    posz # future use
    dirx # future use
    diry # future use
    dirz # future use
    pulselengthtable # Available pulse lengths for the channel [s]
    gaintable # Gain for each pulse length in the PulseLengthTable [dB]
    sacorrectiontable # Sa correction for each pulse length in the
                      # PulseLengthTable [dB]
    gptsoftwareversion
end

function readconfigurationtransducer(stream::IO)
    channelid = readstring(stream, 128) # Channel identification
    beamtype = readlong(stream) # 0 = Single, 1 = Split
    frequency = readfloat(stream) # [Hz]
    gain = readfloat(stream) # [dB] - See note below
    equivalentbeamangle = readfloat(stream) # [dB]
    beamwidthalongship = readfloat(stream) # [degree]
    beamwidthathwartship = readfloat(stream) # [degree]
    anglesensitivityalongship= readfloat(stream)
    anglesensitivityathwartship = readfloat(stream)
    angleoffsetalongship = readfloat(stream) # [degree]
    angleoffsetathwartship = readfloat(stream) # [degree]
    posx = readfloat(stream) # future use
    posy = readfloat(stream) # future use
    posz = readfloat(stream) # future use
    dirx = readfloat(stream) # future use
    diry = readfloat(stream) # future use
    dirz = readfloat(stream) # future use
    pulselengthtable = readfloats(stream, 5) # Available pulse lengths for the channel [s]
    spare1 = readchars(stream, 8) # future use
    gaintable = readfloats(stream, 5) # Gain for each pulse length in the PulseLengthTable [dB]
    spare2 = readchars(stream, 8) # future use
    sacorrectiontable = readfloats(stream, 5) # Sa correction for each pulse length in the PulseLengthTable [dB]
    spare3 = readchars(stream, 8);
    gptsoftwareversion =readstring(stream, 16);
    spare4 = readchars(stream, 28);

    ConfigurationTransducer(channelid, beamtype, frequency, gain,
                            equivalentbeamangle, beamwidthalongship,
                            beamwidthathwartship,
                            anglesensitivityalongship,
                            anglesensitivityathwartship,
                            angleoffsetalongship,
                            angleoffsetathwartship, posx, posy, posz,
                            dirx, diry, dirz, pulselengthtable,
                            gaintable,
                            sacorrectiontable,
                            gptsoftwareversion)

end

#

struct ConfigurationDatagram <: Datagram
    dgheader::DatagramHeader
    configurationheader::ConfigurationHeader
    configurationtransducers::Array{ConfigurationTransducer}
end

function readconfigurationdatagram(stream::IO, dgheader::DatagramHeader)
    configurationheader = readconfigurationheader(stream);

    n =  configurationheader.transducercount
    configurationtransducers = []
    for i = 1:n
        configurationtransducer = readconfigurationtransducer(stream)
        push!(configurationtransducers, configurationtransducer)
    end

    ConfigurationDatagram(dgheader,
                          configurationheader,
                          configurationtransducers)
end

end # module


text:
@testset "selecttransform" begin
    rng = MersenneTwister(123456)

    select = [1, 3, 5]
    t = SelectTransform(select)

    x_cols = ColVecs(randn(rng, maximum(select), 6))
    x_rows = RowVecs(randn(rng, 4, maximum(select)))

    @testset "$(typeof(x))" for x in [x_cols, x_rows]
        x′ = map(t, x)
        @test all([t(x[n]) == x[n][select] for n in eachindex(x)])
        @test all([t(x[n]) == x′[n] for n in eachindex(x)])
    end

    select2 = [2, 3, 5]
    KernelFunctions.set!(t, select2)
    @test t.select == select2

    @test repr(t) == "Select Transform (dims: $(select2))"
    test_ADs(()->transform(SEKernel(), SelectTransform([1,2])))
end

ext: jl


text:
__precompile__()

module OrthographicProjection

	using Common
	import Common.AABB, Common.LAR, Common.Points
	import Common.getmodel
	using FileManager
	using Printf
	using Images

	#
	include("struct.jl")
	include("trie_traversal.jl")
	include("saves.jl")
	include("main.jl")
	include("process_point.jl")
	include("util.jl")

	export Common, FileManager
end # module


text:
# ---
# title: 414. Third Maximum Number
# id: problem414
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/third-maximum-number/description/>
# hidden: true
# ---
#
# Given a **non-empty** array of integers, return the **third** maximum number
# in this array. If it does not exist, return the maximum number. The time
# complexity must be in O(n).
#
# **Example 1:**
#
#
#
#     Input: [3, 2, 1]
#
#     Output: 1
#
#     Explanation: The third maximum is 1.
#
#
# **Example 2:**
#
#
#
#     Input: [1, 2]
#
#     Output: 2
#
#     Explanation: The third maximum does not exist, so the maximum (2) is returned instead.
#
#
# **Example 3:**
#
#
#
#     Input: [2, 2, 3, 1]
#
#     Output: 1
#
#     Explanation: Note that the third maximum here means the third maximum distinct number.
#     Both numbers with value 2 are both considered as second maximum.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using IndependentHypothesisWeighting
using Test

@testset "IndependentHypothesisWeighting.jl" begin
    include("test_weighted_multiple_testing.jl")
    include("test_ihw.jl")
    include("test_grenander.jl")
end


text:
# Test for nep helper creators
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "NEPTypes: creat help" begin
    # Mder_NEP
    function my_Mder(s,der)
        A0=ones(Float64,3,3);
        A1=ones(Float64,3,3)*3; A1[2,3]=0;
        return (3^der)*exp(3*s)*A0+factorial(der)*((-1)^der)*A1/(3.3+s)^(der+1)
    end
    AA0=ones(Float64,3,3);
    AA1=ones(Float64,3,3)*3; AA1[2,3]=0;
    nep_ref=SPMF_NEP([AA0,AA1],[s->exp(3*s), S->inv(3.3*I+S)]);
    nep=Mder_NEP(3,my_Mder)

    λ=3.3;
    norm(compute_Mder(nep,λ,4)-compute_Mder(nep_ref,λ,4))

    λ=-1.2+0.2im;
    X=[ -0.715845   0.865534   0.254796 ; -0.856405  -0.482516   0.0265129;
  1.02593   -0.62892   -2.09615  ]
    @test norm(compute_Mlincomb(nep,λ,X)-compute_Mlincomb(nep_ref,λ,X))<eps()*100

    # Mder_Mlincomb_NEP
    function my_Mlincomb(s,X)
        A0=ones(Float64,3,3);
        A1=ones(Float64,3,3)*3; A1[2,3]=0;
        z=zeros(size(A0,1))
        for k=1:size(X,2)
            der=k-1;
            z+=(3^der)*exp(3*s)*(A0*X[:,k]);
            z+=(factorial(der)*((-1)^der)/(3.3+s)^(der+1))*(A1*X[:,k]);
        end
        return z;
    end
    nep_mlincomb=Mder_Mlincomb_NEP(3,my_Mder, my_Mlincomb)
    @test norm(compute_Mlincomb(nep,λ,X)-compute_Mlincomb(nep_mlincomb,λ,X))<eps()*100
end


text:
#= Formatting template compose functions =#
const LATEX_BACKSLASH_ESCAPE = Set("&%\$#_{}")
const LATEX_OTHER_ESCAPE = Dict(
    '~' => "{\\textasciitilde}",
    '^' => "{\\textasciicircum}",
    '\\' => "{\\textbackslash}")
function escapelatex(text)
    result = Char[]
    for c in text
        if c ∈ LATEX_BACKSLASH_ESCAPE
            push!(result, '\\')
            push!(result, c)
        elseif haskey(LATEX_OTHER_ESCAPE, c)
            append!(result, LATEX_OTHER_ESCAPE[c])
        else
            push!(result, c)
        end
    end
    String(result)
end

romanfont(s) = "\\mathrm{$s}"

#= Available formatting templates =#
const REQUIREMENTS = Dict(
    :finance => FormatSpecification([
        ParenthesizeNegative()]),
    :us => FormatSpecification([
        DigitSeparator(","),
        DecimalSeparator("."),
        CurrencySymbol(location=:before)]),
    :european => FormatSpecification([
        DigitSeparator("."),
        DecimalSeparator(","),
        CurrencySymbol(location=:after)]),
    :indian => FormatSpecification([
        DigitSeparator(",", (3, 2)),
        DecimalSeparator("."),
        CurrencySymbol(location=:before)]),
    :brief => FormatSpecification([
        CurrencySymbol(symtype=:short, spacing=:none, glued=:require)]),
    :latex => FormatSpecification([
        CurrencySymbol(compose=[escapelatex, romanfont]),
        RenderAs(:minus_sign, Dict("-" => 0)),
        RenderAs(:zero_dash, Dict("\\textrm{---}" => 0)),
        RenderAs(:thin_space, Dict("\\," => 0))]),
    :plain => FormatSpecification([]),
    :local => FormatSpecification([
        CurrencySymbol(location=:dependent)]))


text:
using FLPDevelopment
using Test

@testset "FLPDevelopment.jl" begin
    # Write your own tests here.
end


text:
using DynaMos
using Base.Test
using SymPy

# Test 2D ground
g = ground2D()
@assert g.mass.val == [Inf]
@assert g.position.val == [0,0]
@assert g == ground2D()

# Test connections
name = "testbody"
free2D(name)
hinge2D(name)
hingeslider2D(name)
c1 = slider2D(name,Sym(0))
c2 = sliderx2D(name)
slidery2D(name)

@assert c1 == c2

#Test bodies
b = Body2D()
@assert b == Body2D()


text:
ENV["PYTHON"] = "/usr/bin/python3"
using PyCall
using PartsBasedDetector: make_estimator, estimate

url = "http://images0.chictopia.com/photos/districtofchic/9466442894/" *
      "black-skinny-jeans-dl1961-jeans-sky-blue-pixie-market-shirt_400.jpg"

filename = basename(url)
run(`wget --continue $url -O $filename`)

root = Pkg.dir("PartsBasedDetector")
model_path = joinpath(root, "deps", "src", "pbd", "models", "mat", "lsp.mat")

estimator = make_estimator(model_path)
candidates = estimate(estimator, filename, n_candidates = 3)

@pyimport matplotlib.pyplot as plt
@pyimport scipy.misc as misc

parts = candidates[1].parts
image = misc.imread(filename)
plt.imshow(image)
height, width = size(image)
plt.xlim([width, 0])
plt.ylim([height, 0])
plt.plot(parts[:, 1], parts[:, 2], "ko")
plt.show()


text:
using Test
using ITensors
using ITensors.NDTensors

@testset "Allow general mixtures of collections of indices" begin
  d = 2
  is = Index.((d, d, d, d, d), ("i", "j", "k", "l", "m"))
  i, j, k, l, m = is
  is1 = ([i, j], k, (l, m))
  is2 = [[i, j], k, (l, m)]
  A = randn(dims(is))
  D = randn(minimum(dims(is)))
  x = randn()
  @test hassameinds(ITensor(i), (i,))
  @test hassameinds(ITensor(Float64, i), (i,))
  @test hassameinds(ITensor(is1), is)
  @test hassameinds(ITensor(is2), is)
  @test hassameinds(ITensor(is1...), is)
  @test hassameinds(ITensor(Float64, is1), is)
  @test hassameinds(ITensor(Float64, is2), is)
  @test hassameinds(ITensor(Float64, is1...), is)
  @test hassameinds(randomITensor(is1), is)
  @test hassameinds(randomITensor(is2), is)
  @test hassameinds(randomITensor(is1...), is)
  @test hassameinds(randomITensor(Float64, is1), is)
  @test hassameinds(randomITensor(Float64, is2), is)
  @test hassameinds(randomITensor(Float64, is1...), is)
  @test hassameinds(ITensor(x, is1), is)
  @test hassameinds(ITensor(x, is2), is)
  @test hassameinds(ITensor(x, is1...), is)
  @test hassameinds(ITensor(Float64, x, is1), is)
  @test hassameinds(ITensor(Float64, x, is2), is)
  @test hassameinds(ITensor(Float64, x, is1...), is)
  @test hassameinds(ITensor(Float64, undef, is1), is)
  @test hassameinds(ITensor(Float64, undef, is2), is)
  @test hassameinds(ITensor(Float64, undef, is1...), is)
  @test hassameinds(ITensor(undef, is1), is)
  @test hassameinds(ITensor(undef, is2), is)
  @test hassameinds(ITensor(undef, is1...), is)
  @test hassameinds(emptyITensor(is1), is)
  @test hassameinds(emptyITensor(is2), is)
  @test hassameinds(emptyITensor(is1...), is)
  @test hassameinds(emptyITensor(Float64, is1), is)
  @test hassameinds(emptyITensor(Float64, is2), is)
  @test hassameinds(emptyITensor(Float64, is1...), is)
  @test hassameinds(ITensor(A, is1), is)
  @test hassameinds(ITensor(A, is2), is)
  @test hassameinds(ITensor(A, is1...), is)
  @test hassameinds(itensor(A, is1), is)
  @test hassameinds(itensor(A, is2), is)
  @test hassameinds(itensor(A, is1...), is)
  @test hassameinds(ITensor(Float64, A, is1), is)
  @test hassameinds(ITensor(Float64, A, is2), is)
  @test hassameinds(ITensor(Float64, A, is1...), is)
  @test hassameinds(itensor(Float64, A, is1), is)
  @test hassameinds(itensor(Float64, A, is2), is)
  @test hassameinds(itensor(Float64, A, is1...), is)
  @test hassameinds(diagITensor(is1), is)
  @test hassameinds(diagITensor(is2), is)
  @test hassameinds(diagITensor(is1...), is)
  @test hassameinds(diagITensor(Float64, is1), is)
  @test hassameinds(diagITensor(Float64, is2), is)
  @test hassameinds(diagITensor(Float64, is1...), is)
  @test hassameinds(diagITensor(D, is1), is)
  @test hassameinds(diagITensor(D, is2), is)
  @test hassameinds(diagITensor(D, is1...), is)
  @test hassameinds(diagITensor(Float64, D, is1), is)
  @test hassameinds(diagITensor(Float64, D, is2), is)
  @test hassameinds(diagITensor(Float64, D, is1...), is)
  @test hassameinds(diagITensor(x, is1), is)
  @test hassameinds(diagITensor(x, is2), is)
  @test hassameinds(diagITensor(x, is1...), is)
  @test hassameinds(diagITensor(Float64, x, is1), is)
  @test hassameinds(diagITensor(Float64, x, is2), is)
  @test hassameinds(diagITensor(Float64, x, is1...), is)
  @test hassameinds(diagitensor(D, is1), is)
  @test hassameinds(diagitensor(D, is2), is)
  @test hassameinds(diagitensor(D, is1...), is)
  @test hassameinds(diagitensor(Float64, D, is1), is)
  @test hassameinds(diagitensor(Float64, D, is2), is)
  @test hassameinds(diagitensor(Float64, D, is1...), is)
  @test hassameinds(delta(is1), is)
  @test hassameinds(delta(is2), is)
  @test hassameinds(delta(is1...), is)
  @test hassameinds(delta(Float64, is1), is)
  @test hassameinds(delta(Float64, is2), is)
  @test hassameinds(delta(Float64, is1...), is)
  @test hasinds(combiner(is1), is)
  @test hasinds(combiner(is2), is)
  @test hasinds(combiner(is1...), is)
end

@testset "Allow general mixtures of collections of QN indices" begin
  d = [QN() => 2]
  is = Index.((d, d, d, d, d), ("i", "j", "k", "l", "m"))
  i, j, k, l, m = is
  is1 = ([i, j], k, (l, m))
  is2 = [[i, j], k, (l, m)]
  A = randn(dims(is))
  D = randn(minimum(dims(is)))
  x = randn()
  @test hassameinds(ITensor(i), (i,))
  @test hassameinds(ITensor(Float64, i), (i,))
  @test hassameinds(ITensor(is1), is)
  @test hassameinds(ITensor(is2), is)
  @test hassameinds(ITensor(is1...), is)
  @test hassameinds(ITensor(Float64, is1), is)
  @test hassameinds(ITensor(Float64, is2), is)
  @test hassameinds(ITensor(Float64, is1...), is)
  @test hassameinds(randomITensor(is1), is)
  @test hassameinds(randomITensor(is2), is)
  @test hassameinds(randomITensor(is1...), is)
  @test hassameinds(randomITensor(Float64, is1), is)
  @test hassameinds(randomITensor(Float64, is2), is)
  @test hassameinds(randomITensor(Float64, is1...), is)
  @test hassameinds(ITensor(x, is1), is)
  @test hassameinds(ITensor(x, is2), is)
  @test hassameinds(ITensor(x, is1...), is)
  @test hassameinds(ITensor(Float64, x, is1), is)
  @test hassameinds(ITensor(Float64, x, is2), is)
  @test hassameinds(ITensor(Float64, x, is1...), is)
  @test hassameinds(ITensor(Float64, undef, is1), is)
  @test hassameinds(ITensor(Float64, undef, is2), is)
  @test hassameinds(ITensor(Float64, undef, is1...), is)
  @test hassameinds(ITensor(undef, is1), is)
  @test hassameinds(ITensor(undef, is2), is)
  @test hassameinds(ITensor(undef, is1...), is)
  @test hassameinds(emptyITensor(is1), is)
  @test hassameinds(emptyITensor(is2), is)
  @test hassameinds(emptyITensor(is1...), is)
  @test hassameinds(emptyITensor(Float64, is1), is)
  @test hassameinds(emptyITensor(Float64, is2), is)
  @test hassameinds(emptyITensor(Float64, is1...), is)
  @test hassameinds(ITensor(A, is1), is)
  @test hassameinds(ITensor(A, is2), is)
  @test hassameinds(ITensor(A, is1...), is)
  @test hassameinds(itensor(A, is1), is)
  @test hassameinds(itensor(A, is2), is)
  @test hassameinds(itensor(A, is1...), is)
  @test hassameinds(ITensor(Float64, A, is1), is)
  @test hassameinds(ITensor(Float64, A, is2), is)
  @test hassameinds(ITensor(Float64, A, is1...), is)
  @test hassameinds(itensor(Float64, A, is1), is)
  @test hassameinds(itensor(Float64, A, is2), is)
  @test hassameinds(itensor(Float64, A, is1...), is)
  @test hassameinds(diagITensor(is1), is)
  @test hassameinds(diagITensor(is2), is)
  @test hassameinds(diagITensor(is1...), is)
  @test hassameinds(diagITensor(Float64, is1), is)
  @test hassameinds(diagITensor(Float64, is2), is)
  @test hassameinds(diagITensor(Float64, is1...), is)
  @test hassameinds(diagITensor(D, is1), is)
  @test hassameinds(diagITensor(D, is2), is)
  @test hassameinds(diagITensor(D, is1...), is)
  @test hassameinds(diagITensor(Float64, D, is1), is)
  @test hassameinds(diagITensor(Float64, D, is2), is)
  @test hassameinds(diagITensor(Float64, D, is1...), is)
  @test hassameinds(diagITensor(x, is1), is)
  @test hassameinds(diagITensor(x, is2), is)
  @test hassameinds(diagITensor(x, is1...), is)
  @test hassameinds(diagITensor(Float64, x, is1), is)
  @test hassameinds(diagITensor(Float64, x, is2), is)
  @test hassameinds(diagITensor(Float64, x, is1...), is)
  @test hassameinds(diagitensor(D, is1), is)
  @test hassameinds(diagitensor(D, is2), is)
  @test hassameinds(diagitensor(D, is1...), is)
  @test hassameinds(diagitensor(Float64, D, is1), is)
  @test hassameinds(diagitensor(Float64, D, is2), is)
  @test hassameinds(diagitensor(Float64, D, is1...), is)
  @test hassameinds(delta(is1), is)
  @test hassameinds(delta(is2), is)
  @test hassameinds(delta(is1...), is)
  @test hassameinds(delta(Float64, is1), is)
  @test hassameinds(delta(Float64, is2), is)
  @test hassameinds(delta(Float64, is1...), is)
  @test hasinds(combiner(is1), is)
  @test hasinds(combiner(is2), is)
  @test hasinds(combiner(is1...), is)

  # With flux
  @test hassameinds(ITensor(QN(), i), (i,))
  @test hassameinds(ITensor(Float64, QN(), i), (i,))
  @test hassameinds(ITensor(QN(), is1), is)
  @test hassameinds(ITensor(QN(), is2), is)
  @test hassameinds(ITensor(QN(), is1...), is)
  @test hassameinds(ITensor(Float64, QN(), is1), is)
  @test hassameinds(ITensor(Float64, QN(), is2), is)
  @test hassameinds(ITensor(Float64, QN(), is1...), is)
  @test hassameinds(randomITensor(QN(), is1), is)
  @test hassameinds(randomITensor(QN(), is2), is)
  @test hassameinds(randomITensor(QN(), is1...), is)
  @test hassameinds(randomITensor(Float64, QN(), is1), is)
  @test hassameinds(randomITensor(Float64, QN(), is2), is)
  @test hassameinds(randomITensor(Float64, QN(), is1...), is)
  @test hassameinds(ITensor(x, QN(), is1), is)
  @test hassameinds(ITensor(x, QN(), is2), is)
  @test hassameinds(ITensor(x, QN(), is1...), is)
  @test hassameinds(ITensor(Float64, x, QN(), is1), is)
  @test hassameinds(ITensor(Float64, x, QN(), is2), is)
  @test hassameinds(ITensor(Float64, x, QN(), is1...), is)
  @test hassameinds(ITensor(Float64, undef, QN(), is1), is)
  @test hassameinds(ITensor(Float64, undef, QN(), is2), is)
  @test hassameinds(ITensor(Float64, undef, QN(), is1...), is)
  @test hassameinds(ITensor(undef, QN(), is1), is)
  @test hassameinds(ITensor(undef, QN(), is2), is)
  @test hassameinds(ITensor(undef, QN(), is1...), is)
  @test_throws ErrorException hassameinds(emptyITensor(QN(), is1), is)
  @test_throws ErrorException hassameinds(emptyITensor(QN(), is2), is)
  @test_throws ErrorException hassameinds(emptyITensor(QN(), is1...), is)
  @test_throws ErrorException hassameinds(emptyITensor(Float64, QN(), is1), is)
  @test_throws ErrorException hassameinds(emptyITensor(Float64, QN(), is2), is)
  @test_throws ErrorException hassameinds(emptyITensor(Float64, QN(), is1...), is)
end

@testset "Test Index collection as Vector of abstract type" begin
  d = 2
  i = Index(d)
  A = randn(d, d)
  T = itensor(A, Index[i', dag(i)])
  @test storage(T) isa NDTensors.Dense{Float64}
  T = itensor(A, Any[i', dag(i)])
  @test storage(T) isa NDTensors.Dense{Float64}

  i = Index([QN() => d])
  A = randn(d, d)
  T = itensor(A, Index[i', dag(i)])
  @test storage(T) isa NDTensors.BlockSparse{Float64}
  T = itensor(A, Any[i', dag(i)])
  @test storage(T) isa NDTensors.BlockSparse{Float64}
end

@testset "Test output types of ITensors.indices" begin
  i = Index(2)
  @test ITensors.indices([i'', i', i]) == Index{Int}[i'', i', i]
  @test ITensors.indices((i'', i', i)) == (i'', i', i)
  @test ITensors.indices(((i'',), (i',), i)) == (i'', i', i)
  @test ITensors.indices(((i'', i'), (i,))) == (i'', i', i)
  @test ITensors.indices([(i'',), (i',), (i,)]) == Index{Int}[i'', i', i]
  @test ITensors.indices(Any[(i'',), (i',), (i,)]) == Index{Int}[i'', i', i]
  @test ITensors.indices([(i'',), (i',), [i]]) == Index{Int}[i'', i', i]
  @test ITensors.indices([(i'',), i', [i]]) == Index{Int}[i'', i', i]
  @test ITensors.indices(Any[(i'',), i', [i]]) == Index{Int}[i'', i', i]
  @test ITensors.indices(((i'',), i', [i])) == Index{Int}[i'', i', i]
end


text:
@testset "771.jewels-and-stones.jl" begin
    @test num_jewels_in_stones("aA", "aAAbbbb") == 3
    @test num_jewels_in_stones("z", "ZZ") == 0
end


text:
module StorageMirrorServer

using Random

using Glob

using Pkg
using Pkg: TOML
using Pkg.Artifacts: artifact_names

using Tar
using CodecZlib: GzipCompressor, GzipDecompressor
using TranscodingStreams: TranscodingStream

using HTTP

using ThreadPools
using IterTools

using Dates
using ProgressMeter

compress(io::IO) = TranscodingStream(GzipCompressor(level = 9), io)
decompress(io::IO) = TranscodingStream(GzipDecompressor(), io)

const default_http_parameters = Dict(
    :retry => true,
    :retries => 2,
    :timeout => 7200,
)

export mirror_tarball, read_packages

include("types.jl")
include("utils/utils.jl")
include("utils/server_utils.jl")
include("utils/registry_utils.jl")

include("mirror_tarball.jl")

end # module


text:
function skip_indent(fst::FST)
    if fst.typ == LITERAL && fst.val == ""
        return true
    elseif fst.typ === NEWLINE || fst.typ === NOTCODE
        return true
    end
    false
end

function walk(f, nodes::Vector{FST}, s::State, indent::Int)
    for (i, n) in enumerate(nodes)
        if n.typ === NEWLINE && i < length(nodes)
            if is_closer(nodes[i+1])
                s.line_offset = nodes[i+1].indent
            elseif !skip_indent(nodes[i+1])
                s.line_offset = indent
            end
        else
            walk(f, n, s)
        end
    end
end

"""
    walk(f, fst::FST, s::State)

Walks `fst` calling `f` on each node.

In situations where descending further into a subtree is not desirable `f`
should return a value other than `nothing`.
"""
function walk(f, fst::FST, s::State)
    stop = f(fst, s)
    (stop != nothing || is_leaf(fst)) && return
    walk(f, fst.nodes, s, fst.indent)
end

function increment_line_offset!(fst::FST, s::State)
    is_leaf(fst) || return
    s.line_offset += length(fst)
    return nothing
end

function add_indent!(fst::FST, s::State, indent)
    indent == 0 && return
    lo = s.line_offset
    f = (fst::FST, s::State) -> begin
        fst.indent += indent
        return nothing
    end
    walk(f, fst, s)
    s.line_offset = lo
end

# unnest, converts newlines to whitespace
function nl_to_ws!(fst::FST, nl_inds::Vector{Int})
    for (i, ind) in enumerate(nl_inds)
        fst[ind] = Whitespace(fst[ind].len)
        i == length(nl_inds) || continue
        pn = fst[ind-1]
        if pn.typ === TRAILINGCOMMA || pn.typ === TRAILINGSEMICOLON
            pn.val = ""
            pn.len = 0
        elseif pn.typ === INVERSETRAILINGSEMICOLON
            pn.val = ";"
            pn.len = 1
        elseif fst.typ === Binary && fst[ind+1].typ === WHITESPACE
            # remove additional indent
            fst[ind+1] = Whitespace(0)
        end
    end
end

function nl_to_ws!(fst::FST, s::State)
    nl_inds = findall(n -> n.typ === NEWLINE && can_nest(n), fst.nodes)
    length(nl_inds) > 0 || return
    margin = s.line_offset + fst.extra_margin + length(fst)
    margin <= s.opts.margin && nl_to_ws!(fst, nl_inds)
    return
end

function dedent!(style::S, fst::FST, s::State) where {S<:AbstractStyle}
    if is_closer(fst) || fst.typ === NOTCODE
        fst.indent -= s.opts.indent
    elseif is_leaf(fst) || fst.typ === StringN
        return
    else
        fst.indent -= s.opts.indent
    end
end

function dedent!(style::YASStyle, fst::FST, s::State)
    if is_closer(fst) || fst.typ === NOTCODE
        fst.indent -= s.opts.indent
    elseif is_leaf(fst) || fst.typ === StringN
        return
    elseif is_unnamed_iterable(fst)
        fst.indent = s.line_offset
        if is_opener(fst[1])
            fst.indent += 1
        end
    elseif is_named_iterable(fst)
        fst.indent = s.line_offset + length(fst[1]) + length(fst[2])
    else
        fst.indent = s.line_offset
    end
end

function unnest!(style::S, fst::FST, s::State) where {S<:AbstractStyle}
    if is_leaf(fst)
        s.line_offset += length(fst)
    end

    dedent!(style, fst, s)

    if is_leaf(fst) || fst.typ === StringN || !can_nest(fst)
        return
    end

    nl_to_ws!(fst, s)

    return
end

function unnest!(style::S) where {S<:AbstractStyle}
    (fst::FST, s::State) -> begin
        unnest!(style, fst, s)
    end
end

"""
    nest_if_over_margin!(
        style,
        fst::FST,
        s::State,
        idx::Int;
        stop_idx::Union{Int,Nothing} = nothing,
    )::Bool

Converts the node at `idx` to a `NEWLINE` if the current margin plus the additional margin
from `fst[idx:stop_idx-1]` is greater than the allowed margin.

If `stop_idx == nothing` the range is `fst[idx:end]`.

Returns whether nesting occurred.
"""
function nest_if_over_margin!(
    style,
    fst::FST,
    s::State,
    idx::Int;
    stop_idx::Union{Int,Nothing} = nothing,
)::Bool
    @assert fst[idx].typ == PLACEHOLDER
    margin = s.line_offset
    if stop_idx === nothing
        margin += sum(length.(fst[idx:end])) + fst.extra_margin
    else
        margin += sum(length.(fst[idx:stop_idx-1]))
    end

    if margin > s.opts.margin ||
       (idx < length(fst.nodes) && is_comment(fst[idx+1])) ||
       (idx > 1 && is_comment(fst[idx-1]))
        fst[idx] = Newline(length = fst[idx].len)
        s.line_offset = fst.indent
        return true
    end

    nest!(style, fst[idx], s)
    return false
end


text:
module HelpersTest

using Test
using ForneyLab: ensureMatrix, isApproxEqual, isRoundedPosDef, huge, tiny, format, leaftypes, cholinv, diageye, *, ^, @symmetrical
using LinearAlgebra: Diagonal, isposdef, I, Hermitian

@testset "Helpers" begin
    @testset "ensureMatrix" begin
        # should convert input argument to a Matrix or Nothing
        @test ensureMatrix([1.0, 2.0]) == Diagonal([1.0, 2.0])
        @test ensureMatrix(Diagonal([1.0, 2.0])) == Diagonal([1.0, 2.0])
        @test ensureMatrix(Matrix(1.0I,2,2)) == Matrix(1.0I,2,2)
        @test ensureMatrix(1.0) == Matrix(1.0I,1,1)
        @test ensureMatrix(nothing) == nothing
    end

    @testset "isApproxEqual" begin
        # should work for scalars, vectors and matrices
        @test isApproxEqual(1.0, 1.0+1e-15) == true
        @test isApproxEqual(1.0, 1.0+1e-9) == false
        @test isApproxEqual([1.0, 1.0], [1.0, 1.0].+1e-15) == true
        @test isApproxEqual([1.0, 1.0], [1.0, 1.0].+1e-9) == false
        @test isApproxEqual(Matrix(1.0I,3,3), Matrix(1.0I,3,3).+1e-15) == true
        @test isApproxEqual(Matrix(1.0I,3,3), Matrix(1.0I,3,3).+1e-9) == false
    end

    @testset "isRoundedPosDef" begin
        # should check positive definiteness with robustness for numerical errors
        K = inv([4.0 3.0 2.0;
                 3.0 4.0 3.0;
                 2.0 3.0 4.0])
        @test isposdef(K) == false
        @test isRoundedPosDef(K) == true
        @test isRoundedPosDef(Diagonal([1.0, 2.0, 3.0])) == true
    end

    @testset "cholinv" begin
        # should perform a matrix inversion on a positive (semi)definite matrix
        A = [2.0 1.0; 1.0 2.0]
        @test cholinv(A) ≈ inv(A)
        B = [2.0 1.0; 1.0 1.0] #PD
        E = [1.0 1.0+1e8; 1.0+1e8 1.0] #non-PD not due to numerical precision
        @test isApproxEqual(cholinv(B), [1.0 -1.0; -1.0 2.0])
        @test_throws Exception cholesky(Nothing, E)
        F = Diagonal([2.0, 3.0])
        @test cholinv(F) == inv(F)
        G = [Diagonal([2.0]) reshape([1.0],1,1); reshape([1.0],1,1) Diagonal([2.0])]
        @test isApproxEqual(cholinv(G), inv(Matrix(G)))
    end

    @testset "diageye" begin
        # should be shorthand for Diagonal(eye(M))
        M = diageye(3)
        @test typeof(M) == Diagonal{Float64, Array{Float64,1}}
        @test M == Diagonal(ones(3))
    end

    @testset "operators for diagonal matrices" begin
        # should yield Diagonal
        D = Diagonal([1.0, 2.0])
        M = [2.0 1.0; 1.0 2.0]
        @test D.*D == Diagonal([1.0, 4.0])
        @test D.*M == Diagonal([2.0, 4.0])
        @test M.*D == Diagonal([2.0, 4.0])
        @test D^0.5 == Diagonal([1.0, sqrt(2.0)])
        @test sqrt(D) == Diagonal([1.0, sqrt(2.0)])
    end

    @testset "format" begin
        # should return formatted strings for various input types
        @test format(0.0000001) == "1.00e-07"
        @test format(0.0) == "0.00"
        @test format(true) == "true"
        @test format(:a) == "a"
        @test format([7.345345456456456464564645645645, 0.00005345, -0.000145, -108.0]) == "[7.35, 5.34e-05, -1.45e-04, -1.08e+02]"
        @test format([7.345345456456456464564645645645 0.00005345; -0.000145 -108.0]) == "[[7.35, 5.34e-05][-1.45e-04, -1.08e+02]]"
    end

    @testset "symmetrical" begin
        @symmetrical foo(a::Int, b::Float64) = 1

        @test foo(1, 1.0) === foo(1.0, 1)

        @symmetrical bar(a::Int, b::Tuple{L, R}) where { L <: Real, R <: Real } = 1

        @test bar(1, Tuple{Float64, Float32}((1.0, 1.0f0))) === bar(Tuple{Float64, Float32}((1.0, 1.0f0)), 1)
        @test bar(1, Tuple{Float32, Float64}((1.0f0, 1.0))) === bar(Tuple{Float32, Float64}((1.0f0, 1.0)), 1)

        @symmetrical function baz(a::Int, b::Float64, c::String) where A where B where C
            return 1
        end

        @test baz(1, 1.0, "1") === baz(1.0, 1, "1")
    end

    @testset "leaftypes" begin
        # should return all subtypes that are leafs on the type tree
        @test Set(leaftypes(Integer)) == Set([BigInt, Bool, UInt128, UInt16, UInt32, UInt64, UInt8, Int128, Int16, Int32, Int64, Int8])
        @test Set(leaftypes(AbstractFloat)) == Set([BigFloat, Float16, Float32, Float64])
    end
end

end # module


text:
import Base: ==, -, *, +, append!, convert, copy, delete!, deleteat!,
  firstindex, getindex, hash, in, isempty, isequal, lastindex, length, merge, merge!,
  push!, read, setindex!, show, size, sizeof, sort!, sort, summary, write
import DSP: filtfilt
import ..SeisIO: BUF,
  FDSNget!,
  FiltDefs,
  FloatArray,
  InstrumentPosition,
  KW,
  SeisData,
  checkbuf!,
  checkbuf_8!,
  checkbuf_strict!,
  clear_notes!,
  code2loctyp,
  code2resptyp,
  code2typ,
  datafields,
  dtconst,
  fdsn_chp,
  fdsn_uhead,
  get_http_req,
  loctyp2code,
  note!,
  read_misc,
  read_string_vec,
  reset_sacbuf,
  resptyp2code,
  sac_nul_c,
  show_os,
  split_id,
  str_trunc,
  sμ,
  tnote,
  typ2code,
  unindexed_fields,
  webhdr,
  write_misc,
  write_sac_channel,
  write_string_vec,
  writesac,
  μs


text:
module ShaleDrillingLikelihood_DrillingDataStructure_Test

using ShaleDrillingLikelihood
using Test
using Random
using StatsBase
using StatsFuns
using InteractiveUtils
using Dates
using BenchmarkTools

using ShaleDrillingLikelihood: SimulationDraws, _u, _v, SimulationDrawsMatrix, SimulationDrawsVector,
    AbstractDrillModel, TestDrillModel,
    ExogTimeVars, _timestamp, _timevars, Quarter,
    j1ptr, j2ptr, tptr, zchars, jtstart, ichars, _y, _x, j1chars, hasj1ptr, maxj1length,
    j1start, j1stop, j1length, j1_range,
    tstart, tstop, trange, tlength,
    zcharsvec,
    AbstractDataDrill, DataDrill,
    ObservationGroup,
    Observation, ObservationDrill,
    DrillUnit,
    AbstractRegimeType, InitialDrilling, DevelopmentDrilling, FinishedDrilling,
    AbstractDrillRegime, DrillInitial, DrillDevelopment, DrillLease,
    _i, _data,
    ObservationGroup,
    # action,
    state,
    uniti,
    simulate_lease,
    randsumtoone,
    num_initial_leases,
    DataDrillStartOnly

@testset "Drilling Data Structure" begin

    @testset "initial vs dev't" begin
        @test InitialDrilling()+1 == InitialDrilling()+100 == DevelopmentDrilling()
        @test InitialDrilling()-1 == nothing
        @test DevelopmentDrilling()+1 == FinishedDrilling()
        @test DevelopmentDrilling()-1 == InitialDrilling()
        @test FinishedDrilling()+1 == nothing
        @test FinishedDrilling()-1 == DevelopmentDrilling()
        @test InitialDrilling() in (InitialDrilling(), DevelopmentDrilling())
    end

    @testset "DataDrill Construction" begin
        nt = 12
        _j1ptr   = [1, 1]
        _j2ptr   = 1:1
        _tptr    = [1, 5]
        _jtstart = [1,]
        _jchars  = zeros(Int,0)
        _ichars  = [(1.0,), ]
        y = collect(1:4)
        x = collect(1:4)
        _zchars  = ExogTimeVars([(x,) for x in randn(nt)], range(Date(2003,10); step=Quarter(1), length=nt))
        _zchars_short  = ExogTimeVars([(x,) for x in randn(2)], range(Date(2003,10); step=Quarter(1), length=2))

        DataDrill(TestDrillModel(), _j1ptr, _j2ptr, _tptr, _jtstart, _jchars, _ichars, y, x, _zchars)
        @test_throws DimensionMismatch DataDrill(TestDrillModel(), [1,], _j2ptr, _tptr, _jtstart, _jchars, _ichars, y, x, _zchars)
        @test_throws ErrorException    DataDrill(TestDrillModel(),       _j1ptr, _j2ptr, _tptr, _jtstart, _jchars, _ichars, y, x, _zchars_short)
    end

    @testset "Iteration over DataDrill" begin
        nt = 12
        _j1ptr   = [1, 2]
        _j2ptr   = 2:2
        _tptr    = [1, 5, 6]
        _jtstart = [1, 5]
        _jchars  = [1.0, ]
        _ichars  = [(1.0,), ]
        _tchars  = [(x,x,) for x in 1:5]
        y = collect(1:5)
        x = copy(y)
        _zchars  = ExogTimeVars([(x,) for x in randn(nt)], range(Date(2003,10); step=Quarter(1), length=nt))

        data = DataDrill(TestDrillModel(), _j1ptr, _j2ptr, _tptr, _jtstart, _jchars, _ichars, y, x, _zchars)

        @test length(data) == 1
        @test maxj1length(data) == 1

        @test j1ptr(data, 1) == 1
        @test j1_range(data,1) == 1:1
        @test j1length(data,1) == 1
        @test j1stop(data,1) == 1
        @test hasj1ptr(data) == true

        @test tstart(data,1) == 1
        @test tstop(data,1) == 4
        @test trange(data,1) == 1:4
        @test tlength(data,1) == 4

        @test tstart(data,2) == 5
        @test tstop(data,2) == 5
        @test trange(data,2) == 5:5
        @test tlength(data,2) == 1

        @test zcharsvec(data, 2) == _timevars(_zchars)[2:end]

        # iterate through all
        i,r,l,t = (0,0,0,0)

        for unit in data
            i+= 1
            @test _i(unit) == i
            @test isa(unit, DrillUnit)
            @test DataDrill(unit) === data
            @test length(unit) == DevelopmentDrilling()
            @test num_initial_leases(unit) == j1length(unit)

            for regime in unit
                r += 1
                @test isa(regime, AbstractDrillRegime)
                @test DataDrill(regime) === data

                for lease in regime
                    @test isa(lease, DrillLease)
                    @test DataDrill(lease) === data
                    l += 1

                    for obs in lease
                        t+= 1
                        @test isa(obs, ObservationDrill)
                        @test ichars(obs) == (1.0,)
                        @test _x(obs) == t
                        @test _y(obs) == t
                        @test zchars(obs) == _zchars[t]
                    end

                end
            end
        end
        @test (i,r,l,t) == (1,2,2,5)
    end

    @testset "Create random DataDrill" begin

        theta = [1.0, 1.0, 1.0, -1.0, 0.5]
        data = DataDrill(TestDrillModel(), theta)

        # go through leases
        j = 0
        t = 0
        for (i,unit) in enumerate(data)
            @test isa(unit, DrillUnit)

            for (r,regimes) in enumerate(unit)
                @test isa(regimes, AbstractDrillRegime)
                @test r in 1:2
                length(regimes) > 0 && (@test sum(j1chars(regimes)) ≈ 1)
                if r == 1
                    @test isa(regimes, DrillInitial)
                else
                    @test isa(regimes, DrillDevelopment)
                end

                for lease in regimes
                    j += 1
                    @test uniti(lease) == i
                    @test isa(lease, DrillLease)

                    for obs in lease
                        t += 1
                    end

                end
            end
        end
        @test j == length(tptr(data))-1
        @test t == length(_y(data))

        @testset "Drilling Start Only" begin
            dso = DataDrillStartOnly(data)
            for (ud, us) in zip(data, dso)
                for (rd, rs) in zip(ud, us)
                    for (ld, ls) in zip(rd, rs)

                        if length(ld) > 0
                            @test length(ls) == 1
                            @test first(_x(ld)) == first(_x(ls))
                            @test first(_y(ld)) == first(_y(ls))
                        else
                            @test length(ls) == 0
                        end

                    end
                end
            end
        end

    end # testset: random drilling data

end # overall testset

end # module


text:
# PCI Utility Functions
function pcipath(bus, device, fn)
    bus_str    = string(bus; base = 16, pad = 2)
    device_str = string(device; base = 16, pad = 2)
    fn_str     = string(fn; base = 16, pad = 1)
    path = joinpath(
        "/proc/bus/pci",
        bus_str,
        "$device_str.$fn_str"
    )

    return path
end

function pcipath(group, bus, device, fn)
    iszero(group) && return pcipath(bus, device, fn)
    error("Cannot deal with non-zero PCI Groups yet")
end

#####
##### PCI Handles
#####

abstract type AbstractPCIHandle end
mutable struct Handle <: AbstractPCIHandle
    fd::IOStream

    # Open the os file - attach a finalizer so the file is appropriately closed.
    function Handle(str::String; read = true, write = true)
        handle = new(open(str; read = read, write = write, lock = false))
        finalizer(close, handle)
        return handle
    end
end
Base.close(P::Handle) = close(P.fd)

# Construct the Handle from a group, bus, device, and offset
#
# This is largely reverse engineered from `pci.cpp` from `https://github.com/opcm/pcm`
Handle(args::Integer...) = Handle(pcipath(args...))
Base.seek(P::Handle, offset::IndexZero) = seek(P.fd, value(offset))

function Base.write(h::Handle, v, offset::IndexZero)
    seek(h, offset)
    bytes = write(h.fd, v)
    flush(h.fd)
    return bytes
end

# Wrapping the `pread` calls is MUCH faster because it only requires a single
# transition into kernel code.
function Base.read(P::Handle, ::Type{T}, offset::IndexZero) where {T}
    # Make a direct `pread` call
    object = Ref{T}()
    _ = ccall(
        :pread,
        Csize_t,
        (Cint, Ptr{Cvoid}, Csize_t, Cint),
        fd(P.fd),
        object,
        sizeof(T),
        value(offset)
    )

    return object[]
end


text:
using MPT
using Test

@testset "MPT.jl" begin
    # Write your tests here.
end


text:
using IntroToJulia
using Base.test

@test 1==1


text:
module CompositionalNetworks

# Usings
using ConstraintDomains
using Dictionaries
import Distances
using Evolutionary
using JuliaFormatter
using Memoization
using OrderedCollections
using Random
using ThreadPools
using ThreadSafeDicts
using Unrolled

export Composition
export ICN
export aggregation_layer
export arithmetic_layer
export code
export comparison_layer
export compose
export compose_to_file!
export composition
export composition_to_file!
export explore
export explore_learn_compose
export hamming
export lazy
export lazy_param
export learn_compose
export manhattan
export max_icn_length
export minkowski
export optimize!
export regularization
export show_layers
export symbols
export transformation_layer

# Include utils
include("utils.jl")
include("metrics.jl")

# Includes layers
include("layer.jl")
include("transformation.jl")
include("arithmetic.jl")
include("aggregation.jl")
include("comparison.jl")

# Include ICN
include("icn.jl")
include("composition.jl")

# Genetic Algorithm and learning
include("genetic.jl")
include("explore.jl")
include("learn.jl")

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to convert osculating elements to mean elements using the J2
#   osculating algorithm.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A; Crawford, P (2008). SGP4 orbit determination. AIAA/AAS
#       Astrodynamics Specialist Conference, Honoulu, HI.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export rv_to_mean_elements_j2osc

"""
    rv_to_mean_elements_j2osc(vJD::AbstractVector{T}, vr::AbstractVector{Tv}, vv::AbstractVector{Tv}, W = I; mean_elements_epoch::Symbol = :end, max_it::Int = 50, j2_gc = sgp4_gc_wgs84, atol::Number = 2e-4, rtol::Number = 2e-4) where {T,Tv<:AbstractVector}

Compute the mean elements for SGP4 based on the position `vr` and velocity
vectors `vr` represented in TEME reference frame. The epoch of those
measurements [Julian Day] must be in `vJD`.

The matrix `W` defined the weights for the least-square algorithm.

# Keywords

* `estimate_bstar`: If `true`, then the BSTAR parameters of the TLE will be
                    estimated.
* `mean_elements_epoch`: If it is  `:end`, the epoch of the mean elements will
                         be equal to the last value in `vJD`. Otherwise, if it
                         is `:begin`, the epoch will be selected as the first
                         value in `vJD`.
* `max_it`: The maximum allowed number of iterations.
* `sgp4_gc`: SPG4 constants (see `SGP4_GravCte`).
* `atol`: The tolerance for the absolute value of the residue. If, at any
          iteration, the residue is lower than `atol`, then the iterations stop.
* `rtol`: The tolerance for the relative difference between the residues. If, at
          any iteration, the relative difference between the residues in two
          consecutive iterations is lower than `rtol`, then the iterations stop.

# Returns

* The epoch of the elements [Julian Day].
* The mean elements for SGP4 algorithm:
    * Semi-major axis [m];
    * Eccentricity [ ];
    * Inclination [rad];
    * Right ascension of the ascending node [rad];
    * Argument of perigee [rad];
    * True anomaly [rad];
    * BSTAR (0 if `estimate_bstar` is `false`).
* The covariance matrix of the mean elements estimation.

"""
function rv_to_mean_elements_j2osc(vJD::AbstractVector{T},
                                   vr::AbstractVector{Tv},
                                   vv::AbstractVector{Tv},
                                   W = I;
                                   mean_elements_epoch::Symbol = :end,
                                   max_it::Int = 50,
                                   j2_gc = j2_gc_egm08,
                                   atol::Number = 2e-4,
                                   rtol::Number = 2e-4) where {T, Tv<:AbstractVector}

    # Number of measurements.
    num_meas = length(vr)

    # Check if the orbit epoch must be the first or the last element.
    if mean_elements_epoch == :end
        r₁ = last(vr)
        v₁ = last(vv)
        epoch = last(vJD)
    else
        r₁ = first(vr)
        v₁ = first(vv)
        epoch = first(vJD)
    end

    # Initial guess of the mean elements.
    #
    # NOTE: x₁ is the previous estimate and x₂ is the current estimate.
    x₁ = SVector{6,T}(r₁[1], r₁[2], r₁[3], v₁[1], v₁[2], v₁[3])
    x₂ = x₁

    # Number of states in the input vector.
    num_states = length(x₁)

    # Covariance matrix.
    P = SMatrix{num_states, num_states, T}(I)

    # Variable to store the last residue.
    σ_i₋₁ = T(NaN)

    # Variable to store how many iterations the residue increased. This is used
    # to account for divergence.
    Δd = 0

    # Header
    @printf("          %10s %20s %20s\n", "Iter.", "Residue", "Res. variation")

    # Loop until the maximum allowed iteration.
    @inbounds for it = 1:max_it
        x₁ = x₂

        # Variables to store the summations to compute the least square fitting
        # algorithm.
        ΣAᵀWA = @SMatrix zeros(num_states, num_states)
        ΣAᵀWb = @SVector zeros(num_states)

        # Variable to store the residue in this iteration.
        σ_i = T(0)

        @views for k = 1:num_meas
            # Obtain the measured ephemerides.
            y = vcat(vr[k], vv[k])

            # Obtain the computed ephemerides considering the current estimate
            # of the mean elements.
            Δt = (vJD[k] - epoch)*86400

            r̂, v̂ = _j2osc_sv(Δt, j2_gc, epoch,
                             x₁[1], x₁[2], x₁[3],
                             x₁[4], x₁[5], x₁[6])

            ŷ = vcat(r̂, v̂)

            # Compute the error.
            b = y - ŷ

            # Compute the Jacobian.
            A = _j2osc_jacobian(Δt, epoch, x₁, ŷ)

            # Accumulation.
            ΣAᵀWA += A'*W*A
            ΣAᵀWb += A'*W*b
            σ_i   += sum(W*(b.^2))
        end

        # Normalize the residue.
        σ_i /= num_meas

        # Update the estimate.
        P  = SMatrix{num_states, num_states, T}(pinv(ΣAᵀWA))
        δx = P*ΣAᵀWb

        # Limit the correction to avoid divergence.
        for i = 1:num_states
            abs(δx[i] / x₁[i]) > 0.01 &&
                (δx = setindex(δx, 0.01 * abs(x₁[i]) * sign(δx[i]), i))
        end

        x₂ = x₁ + δx

        # Compute the residue variation.
        σ_p = (σ_i - σ_i₋₁)/σ_i₋₁

        if it != 1
            @printf("PROGRESS: %10d %20g %20g %%\n", it, σ_i, 100σ_p)
        else
            @printf("PROGRESS: %10d %20g %20s\n", it, σ_i, "---")
        end

        # Check if the residue is increasing.
        if σ_i < σ_i₋₁
            Δd = 0
        else
            Δd += 1
        end

        # If the residue increased by three iterations and the residue is higher
        # than 5e8, then we abort because the iterations are diverging.
        ( (Δd ≥ 3) && (σ_i > 5e11) ) && error("The iterations diverged!")

        # Check if the condition to stop has been reached.
        ((abs(σ_p) < rtol) || (σ_i < atol) || (it ≥ max_it)) && break

        σ_i₋₁ = σ_i
    end

    # Convert the state vector to mean elements.
    orb = rv_to_kepler(x₂[1], x₂[2], x₂[3], x₂[4], x₂[5], x₂[6], epoch)

    # Return the mean elements for J2 osculating algorithm and the covariance
    # matrix.
    return orb, P
end

################################################################################
#                              Private functions
################################################################################

# Compute the J2 osculating algorithm considering all variables in a state
# vector.
function _j2osc_sv(Δt::Number,
                   j2_gc::J2_GravCte,
                   epoch::Number,
                   rx_i::Number,
                   ry_i::Number,
                   rz_i::Number,
                   vx_i::Number,
                   vy_i::Number,
                   vz_i::Number)

    r_i = @SVector [rx_i, ry_i, rz_i]
    v_i = @SVector [vx_i, vy_i, vz_i]

    # Obtain the required mean elements to initialize the J2 osculating
    # propagator.
    orb_i = rv_to_kepler(r_i, v_i, epoch)

    j2oscd = j2osc_init(orb_i.t, orb_i.a, orb_i.e, orb_i.i, orb_i.Ω, orb_i.ω,
                        orb_i.f, 0, 0; j2_gc = j2_gc)
    r, v = j2osc!(j2oscd, Δt)

    # Return the elements using SI units.
    return r, v
end

function _j2osc_jacobian(Δt::T,
                         epoch::T,
                         x₁::SVector{NS, T},
                         y₁::SVector{NO, T};
                         pert::T = 1e-3,
                         perttol::T = 1e-5,
                         j2_gc::J2_GravCte = j2_gc_egm08) where {NS, NO, T}

    num_states = NS
    dim_obs    = NO
    M          = zeros(T, dim_obs, num_states)

    # Auxiliary variables.
    x₂ = copy(x₁)

    @inbounds for j = 1:num_states
        # State that will be perturbed.
        α = x₂[j]

        # Obtain the perturbation, taking care to avoid small values.
        ϵ = T(0)
        pert_i = pert

        for _ = 1:5
            ϵ = α * pert_i
            abs(ϵ) > perttol && break
            pert_i *= 1.4
        end

        α += ϵ

        # Avoid division by zero in cases that α is very small. In this
        # situation, we force |α| = perttol.
        abs(α) < perttol && (α = sign(α) * perttol)

        # Modify the perturbed state.
        x₂ = setindex(x₂, α, j)

        # Obtain the Jacobian by finite differentiation.
        r, v = _j2osc_sv(Δt, j2_gc, epoch, x₂...)
        y₂   = vcat(r, v)

        M[:,j] .= (y₂ .- y₁)./ϵ

        # Restore the value of the perturbed state for the next iteration.
        x₂ = setindex(x₂, x₁[j], j)
    end

    return M
end


text:
using Test
using Random
using StaticBitSets

@testset "Out of Bounds" begin

    for N in 1:10, T in (UInt8,UInt16,UInt32,UInt64,UInt128)
        x = SBitSet{N,T}()
        min = 1
        max = N*8*sizeof(T)
        @test_throws AssertionError push(x,max+1)
        @test_throws AssertionError push(x,max+10)
        @test_throws AssertionError push(x,min-1)
        @test_throws AssertionError push(x,-10)
    end

end

@testset "Type compatibilities" begin

    combinations = [(1,UInt128),(2,UInt64),(4,UInt32),(8,UInt16),(16,UInt8)]
    max = 128

    data = Vector{Vector{UInt64}}(collect(eachrow(rand(1:max,100,max÷2))))
    bitsets = [SBitSet{N,T}.(data) for (N,T) in combinations]

    for i in 1:length(combinations)-1
        @test collect.(bitsets[i]) == collect.(bitsets[i+1])
    end

end

@testset "SBitSet{$N,$T}, to, from, push, pop and membership" for N in 1:4, T in (UInt8,UInt16,UInt32,UInt64,UInt128)

    max = N*8*sizeof(T)
    for round in 1:100
        content = Set()
        bs = SBitSet{N,T}()
        for i in rand(1:max,N*N)
            push!(content,i)
            bs = push(bs,i)
            @test i ∈ bs
            @test Set(collect(bs)) == content
        end
        for i in rand(1:max,N*N)
            pop!(content,i,nothing)
            bs = pop(bs,i)
            @test i ∉ bs
            @test Set(collect(bs)) == content
        end
    end
end

@testset "SBitSet{$N,$T}, bitwise ops, empty, isdisjoint and length" for N in 1:4, T in (UInt8,UInt16,UInt32,UInt64,UInt128)

    max = N*8*sizeof(T)

    data = Vector{Vector{UInt64}}(collect(eachrow(rand(1:max,100,4*N))))
    sets = Set.(data)
    bitsets = SBitSet{N,T}.(data)

    @test sets == Set.(collect.(bitsets))

    for (i,j) in collect(Iterators.product(1:length(data), 1:length(data)))
        @test sets[i] ∩ sets[j] == Set(bitsets[i] ∩ bitsets[j])
        @test setdiff(sets[i],sets[j]) ∪ setdiff(sets[j],sets[i]) == Set(xor(bitsets[i],bitsets[j]))
        @test isempty(sets[i] ∩ sets[j]) == isempty(Set(bitsets[i] ∩ bitsets[j]))
        @test isempty(sets[i] ∩ sets[j]) == isdisjoint(sets[i],sets[j])
        @test (length(sets[i] ∩ sets[j]) == 1) == singleton_intersection(bitsets[i],bitsets[j])
        @test length(sets[i] ∩ sets[j]) == length(Set(bitsets[i] ∩ bitsets[j]))
        @test sets[i] ∪ sets[j] == Set(bitsets[i] ∪ bitsets[j])
        @test isempty(sets[i] ∪ sets[j]) == isempty(Set(bitsets[i] ∪ bitsets[j]))
        @test length(sets[i] ∪ sets[j]) == length(Set(bitsets[i] ∪ bitsets[j]))
        @test setdiff(sets[i],sets[j]) == Set(bitsets[i] ~ bitsets[j])
        @test (sets[i] ⊆ sets[j]) == (bitsets[i] ⊆  bitsets[j])
    end

end


text:
import SymbolicUtils.Rewriters: RestartedChain
using DataStructures

export semipolynomial_form, semilinear_form, semiquadratic_form, polynomial_coeffs

## BoundedDegreeMonomial helper type

"""
    BoundedDegreeMonomial

# Attrtibutes
- p -- a monomial
- coeff -- the coefficient (monomial means p * coeff)
- overdegree -- boolean flag shows if degree of p * coeff
  is too high (this is marked by the semi-polynomialization process)

"""
struct BoundedDegreeMonomial
    p::Union{Mul, Pow, Int, Sym, Term}
    coeff::Any
    overdegree::Bool
end

highdegree(x) = BoundedDegreeMonomial(1, x, true)

highdegree(x::BoundedDegreeMonomial) = (@assert(x.overdegree); x)

SymbolicUtils.symtype(b::BoundedDegreeMonomial) = symtype(b.p)

isop(x, op) = istree(x) && operation(x) === op

isop(op) = x -> isop(x, op)

pdegree(x::Mul) = sum(values(x.dict))
pdegree(x::Union{Sym, Term}) = 1
pdegree(x::Pow) = pdegree(x.base) * x.exp
pdegree(x::Number) = 0

# required by unsorted_arguments etc.
SymbolicUtils.unstable_pow(x::BoundedDegreeMonomial, i::Integer) = (@assert(i==1); x)

_degree(x::BoundedDegreeMonomial) = x.overdegree ? Inf : pdegree(x.p)

_degree(x) = isop(x, *) ? sum(_degree, unsorted_arguments(x)) : 0

## Semi-polynomial form

isboundedmonom(x) = x isa BoundedDegreeMonomial && !x.overdegree
bareterm(x, f, args;kw...) = Term{symtype(x)}(f, args)

function mark_and_exponentiate(expr, vars, deg)

    # Step 1
    # Mark all the interesting variables -- substitute without recursing into nl forms

    expr′ = mark_vars(expr, vars)

    # Step 2
    # Construct and propagate BoundedDegreeMonomial for ^ and *

    rules = [@rule (~a::isboundedmonom) ^ (~b::(x-> x isa Integer)) =>
             BoundedDegreeMonomial(((~a).p)^(~b), (~a).coeff ^ (~b), ~b > deg)
             @rule (~a::isop(+)) ^ (~b::(x -> x isa Integer)) => pow_of_add(~a, ~b, deg, vars)

             @rule *(~~x) => mul_bounded(~~x, deg)]

    expr′ = Postwalk(RestartedChain(rules), similarterm=bareterm)(expr′)
end

function semipolyform_terms(expr, vars::OrderedSet, deg)
    if deg <= 0
        throw(ArgumentError("Degree for semi-polynomial form must be > 0"))
    end

    expr′ = mark_and_exponentiate(expr, vars, deg)

    # Step 3: every term now show have at most one valid monomial -- find the coefficient for it.
    Postwalk(Chain([@rule *(~~x, ~a::isboundedmonom, ~~y) =>
                                BoundedDegreeMonomial((~a).p,
                                                      (~a).coeff * _mul(~~x) * _mul(~~y),
                                                      (~a).overdegree)]),
             similarterm=bareterm)(expr′)

end

function mark_vars(expr, vars)
    if expr in vars
        return BoundedDegreeMonomial(expr, 1, false)
    end

    if !istree(expr)
        return expr
    else
        op = operation(expr)
        args = unsorted_arguments(expr)

        if op === (+) || op === (*)
            return Term{symtype(expr)}(op, map(x -> mark_vars(x, vars), args))
        elseif op === (^)
            base, exp = arguments(expr)
            Term{symtype(expr)}(^, [mark_vars(base, vars), exp])
        elseif length(args) == 1
            if linearity_1(op)
                return Term{symtype(expr)}(op, map(x -> mark_vars(x, vars), args))
            else
                return expr
            end
        else
            return expr
        end
    end
end

function bifurcate_terms(expr)
    # Step 4: Bifurcate polynomial and nonlinear parts:

    if expr isa BoundedDegreeMonomial
        return isboundedmonom(expr) ?
            (Dict(expr.p => expr.coeff), 0) :
            (Dict(), expr.p * expr.coeff)
    elseif istree(expr) && operation(expr) == (+)
        args = collect(all_terms(expr))
        polys = filter(isboundedmonom, args)

        poly = Dict()
        for p in polys
            if haskey(poly, p.p)
                poly[p.p] += p.coeff
            else
                poly[p.p] = p.coeff
            end
        end

        nls = filter(!isboundedmonom, args)
        nl = cautious_sum(nls)

        return poly, nl
    else
        return Dict(), expr
    end
end

function init_semipoly_vars(vars)
    set = OrderedSet(unwrap.(vars))
    @assert length(set) == length(vars) "vars passed to semi-polynomial form must be unique"
    set
end

## API:
#

"""
    semipolynomial_form(expr, vars, degree)

Semi-polynomial form. Returns a tuple of two objects:

1. A dictionary of coefficients keyed by monomials in `vars` upto the given `degree`,
2. A residual expression which has all terms not represented as a product of monomial and a coefficient

    semipolynomial_form(exprs::AbstractArray, vars, degree)

For every expression in `exprs` computes the semi-polynomial form as above and
returns a tuple of two objects -- a vector of coefficient dictionaries,
and a vector of residual terms.
"""
function semipolynomial_form(expr, vars, degree)
    expr = unwrap(expr)
    vars = init_semipoly_vars(vars)

    bifurcate_terms(semipolyform_terms(expr, vars, degree))
end

function semipolynomial_form(exprs::AbstractArray, vars, degree)
    exprs = unwrap.(exprs)
    vars = init_semipoly_vars(vars)

    matches = map(x -> semipolyform_terms(x, vars, degree), exprs)
    tmp = map(bifurcate_terms, matches)
    dicts, nls = map(first, tmp), map(last, tmp)
end

"""
    polynomial_coeffs(expr, vars)

Find coefficients of a polynomial in `vars`.

Returns a tuple of two elements:
1. A dictionary of coefficients keyed by monomials in `vars`
2. A residual expression which is the constant term

(Same as `semipolynomial_form(expr, vars, Inf)`)
"""
polynomial_coeffs(expr, vars) = semipolynomial_form(expr, vars, Inf)

"""
    semilinear_form(exprs::AbstractVector, vars::AbstractVector)

Returns a tuple of a sparse matrix `A`, and a residual vector `c` such that,

`A * vars + c` is the same as `exprs`.
"""
function semilinear_form(exprs::AbstractArray, vars)
    exprs = unwrap.(exprs)
    vars = init_semipoly_vars(vars)
    ds, nls = semipolynomial_form(exprs, vars, 1)

    idxmap = Dict(v=>i for (i, v) in enumerate(vars))

    I = Int[]
    J = Int[]
    V = Num[]

    for (i, d) in enumerate(ds)
        for (k, v) in d
            push!(I, i)
            push!(J, idxmap[k])
            push!(V, v)
        end
    end

    sparse(I,J,V, length(exprs), length(vars)), wrap.(nls)
end

"""
    semilinear_form(exprs::AbstractVector, vars::AbstractVector)

Returns a tuple of 4 objects:

1. a matrix `A` of dimensions (m x n)
2. a matrix `B` of dimensions (m x (n+1)*n/2)
3. a vector `v2` of length (n+1)*n/2 containing monomials of `vars` upto degree 2 and zero where they are not required.
4. a residual vector `c` of length m.

where `n == length(exprs)` and `m == length(vars)`.

The result is arranged such that, `A * vars + B * v2 + c` is the same as `exprs`.
"""
function semiquadratic_form(exprs, vars)
    exprs = unwrap.(exprs)
    vars = init_semipoly_vars(vars)
    ds, nls = semipolynomial_form(exprs, vars, 2)

    idxmap = Dict(v=>i for (i, v) in enumerate(vars))

    m, n = length(exprs), length(vars)
    I1 = Int[]
    J1 = Int[]
    V1 = Num[]

    I2 = Int[]
    J2 = Int[]
    V2 = Num[]

    v2_I = Int[]
    v2_V = Num[]

    for (i, d) in enumerate(ds)
        for (k, v) in d
            if pdegree(k) == 1
                push!(I1, i)
                push!(J1, idxmap[k])
                push!(V1, v)
            elseif pdegree(k) == 2
                push!(I2, i)
                if isop(k, ^)
                    b, e = arguments(k)
                    @assert e == 2
                    q = idxmap[b]
                    j = div(q*(q+1), 2)
                    push!(J2, j) # or div(q*(q-1), 2) + q
                    push!(V2, v)
                else
                    @assert isop(k, *)
                    a, b = unsorted_arguments(k)
                    p, q = extrema((idxmap[a], idxmap[b]))
                    j = div(q*(q-1), 2) + p
                    push!(J2, j)
                    push!(V2, v)
                end
                push!(v2_I, j)
                push!(v2_V, k)
            else
                error("This should never happen")
            end
        end
    end

    #v2 = SparseVector(div(n * (n + 1), 2), v2_I, v2_V) # When it works in the future
    # until then
    v2 = zeros(Num, div(n * (n + 1), 2))
    v2[v2_I] .= v2_V

    tuple(sparse(I1,J1,V1, m, n),
          sparse(I2,J2,V2, m, div(n * (n + 1), 2)),
          v2,
          wrap.(nls))
end

## Utilities

function partition(n, parts)
    parts == 1 && return n
    [[[i, p...] for p in partition(n-i, parts-1)]
     for i=0:n] |> Iterators.flatten |> collect
end

all_terms(x) = istree(x) && operation(x) == (+) ? collect(Iterators.flatten(map(all_terms, unsorted_arguments(x)))) : (x,)

unwrap_bp(x::BoundedDegreeMonomial) = x.p * x.coeff
function unwrap_bp(x)
    x = unwrap(x)
    istree(x) ? similarterm(x, operation(x), map(unwrap_bp, unsorted_arguments(x))) : x
end

cautious_sum(nls) = isempty(nls) ? 0 : isone(length(nls)) ? unwrap_bp(first(nls)) : sum(unwrap_bp, nls)

_mul(x) = isempty(x) ? 1 : isone(length(x)) ? first(x) : prod(x)

function mul(a, b, deg)
    if isop(a, +)
        return Term{symtype(a)}(+, map(x->mul(x, b, deg), unsorted_arguments(a)))
    elseif isop(b, +)
        return Term{symtype(a)}(+, mul.((a,), unsorted_arguments(b), deg))
    elseif a isa BoundedDegreeMonomial
        return BoundedDegreeMonomial(a.p, a.coeff * b, a.overdegree)
    elseif b isa BoundedDegreeMonomial
        return BoundedDegreeMonomial(b.p, a * b.coeff, b.overdegree)
    else
        return a * b
    end
end

function mul(a::BoundedDegreeMonomial, b::BoundedDegreeMonomial, deg)
    if a.overdegree || b.overdegree || pdegree(a.p) + pdegree(b.p) > deg
        highdegree(a.p * b.p * a.coeff * b.coeff)
    else
        BoundedDegreeMonomial(a.p * b.p, a.coeff * b.coeff, false)
    end
end

function mul_bounded(xs, deg)
    length(xs) == 1 ?
        first(xs) :
        reduce((x,y) -> mul(x,y, deg), xs)
end

function pow(a::BoundedDegreeMonomial, b, deg)
    return BoundedDegreeMonomial((a.p)^b, a.coeff ^ b, b > deg)
end

pow(a, b, deg) = a^b

function pow_of_add(a, b, deg, vars)
    b == 0 && return 1
    b == 1 && return a

    @assert isop(a, +)
    within_deg(x) = (d=_degree(a);d <= deg)

    a = mark_and_exponentiate(a, vars, deg)
    args = all_terms(a)
    mindeg = minimum(_degree, args)
    maxdeg = maximum(_degree, args)

   #if maxdeg * b <= deg
   #    return mark_and_exponentiate(expand(unwrap_bp(a^b)), vars, deg)
   #end

    if mindeg * b > deg
        return highdegree(unwrap_bp(a^b)) # That's pretty high degree
    end

    interesting = filter(within_deg, args)
    nls = filter(!within_deg, args)

    if isempty(interesting)
        a^b # Don't need to do anything!
    else
        q = partial_multinomial_expansion(interesting, b, deg)
        if maxdeg * b <= deg
            q # q is the whole enchilada
        else
            Term{Real}(+, [all_terms(q)...,
                           map(highdegree, all_terms(unwrap_bp(a^b - q)))...])
        end
    end
end

function partial_multinomial_expansion(xs, exp, deg)
    zs = filter(iszero∘_degree, xs)
    nzs = filter((!iszero)∘_degree, xs)
    if isempty(zs)
        terms = nzs
    else
        terms = [+(zs...), nzs...]
    end
    degs = map(_degree, terms)

    q = []
    nfact = factorial(exp)
    for ks in partition(exp, length(terms))
        td = sum(degs .* ks)

        coeff = div(nfact, prod(factorial, ks))
        if td == 0
            push!(q, *(coeff, (x^y for (x, y) in zip(terms, ks) if !iszero(y))...))
        elseif td <= deg
            push!(q,
                  mul_bounded(vcat(coeff,
                                   [pow(x,y, deg) for (x, y) in zip(terms, ks) if !iszero(y)]),
                              deg))
        end
    end
    return Term{Real}(+, q)
end


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "PIPotential"  begin

#---

using ACE1, JuLIP, Test
using ACE1: combine

#---

@info("Testing :Any Potential to evaluate energy, forces, virial")
basis = ACE1.Utils.rpi_basis(species = :Any, N = 3, r0 = 2.7, rcut = 5.0, maxdeg = 8)
c = ACE1.Random.randcoeffs(basis)
V = combine(basis, c)

#---

at1 = rattle!(bulk(:Fe, cubic=true) * (2,2,2), 0.1)
at2 = deepcopy(at1); at2.Z[:] .= AtomicNumber(:Al)
at3 = deepcopy(at1); at2.Z[1:3:end] .= AtomicNumber(:Al)

for f in (energy, forces, virial)
   println(@test f(V, at1) ≈ f(V, at2) ≈ f(V, at3))
end

#---

@info(" ... check that a non-Any potential will fail this test")
basis1 = ACE1.Utils.rpi_basis(species = :Fe, N = 3, r0 = 2.7, rcut = 5.0, maxdeg = 8)
V1 = combine(basis1, c)   # same coefficients

println(@test energy(V1, at1) ≈ energy(V, at1))
for at in (at2, at3)
   println(@test (
      try
         energy(V1, at2)
         false
      catch
         true
      end
   ))
end

#---

end


text:
"""
    step particle system
        solves 1-step feasibility problem
"""
function step(q1, q2, u1, h;
    tol = 1.0e-8, max_iter = 100)
    # 1-step optimization problem:
    #     find z
    #     s.t. r(z) = 0
    #
    # z = (q, n, sϕ, sb, db)
    #     s are slack variables for convenience
    #     sb[2:3] is the friction force

    # initialize
    z = 1.0 * ones(11)
    z[1:3] = copy(q2)

    # initialize soc variables
    sb, _ = κ_so(z[6:8])
    sb[1] += 1.0
    db, _ = κ_so(z[9:11])
    db[1] += 1.0

    z[6:8] = sb
    z[9:11] = db

    μ = 1.0 # barrier parameter
    flag = false

    θ = [q1; q2; u1]

    e = [1.0; 0.0; 0.0]

    function r(z, θ)
        # system variables
        q3 = view(z, 1:3)
        n = z[4]
        sϕ = z[5]
        sb = view(z, 6:8)
        db = view(z, 9:11)

        q1 = view(θ, 1:3)
        q2 = view(θ, 4:6)
        u1 = view(θ, 7:9)

        λ = [view(sb, 2:3); n] # contact forces
        ϕ = signed_distance(model, q3) # signed-distance function
        vT = (view(q3, 1:2) - view(q2, 1:2)) ./ h

        # action optimality conditions
        [dynamics(model, q1, q2, q3, u1, λ, h);
         sϕ - ϕ;
         n * sϕ - μ;

         # maximum dissipation optimality conditions
         vT - view(db, 2:3);
         sb[1] - model.friction_coeff * n;
         cone_product(db, sb) - μ * e]
    end

    # Jacobian
    function Rz(z, θ)
        # differentiate r
        _r(w) = r(w, θ)
        _R = ForwardDiff.jacobian(_r, z)
        return _R
    end

    function Rθ(z, θ)
        # differentiate r
        _r(w) = r(z, w)
        _R = ForwardDiff.jacobian(_r, θ)

        return _R
    end

    function check_variables(z)
        # system variables
        q3 = view(z, 1:3)
        n = z[4]
        sϕ = z[5]
        sb = view(z, 6:8)
        db = view(z, 9:11)

        if n <= 0.0
            return true
        end

        if sϕ <= 0.0
            return true
        end

        if !κ_so(sb)[2]
            # println("bs not in cone")
            return true
        end

        if !κ_so(db)[2]
            # println("bz not in cone")
            return true
        end

        return false
    end

    for k = 1:4
        extra_iters = 0
        for i = 1:max_iter
            # compute residual, residual Jacobian
            res = r(z, θ)
            if norm(res) < tol
                # println("   iter ($i) - norm: $(norm(res))")
                # return z, true
                flag = true
                continue
            end

            jac = Rz(z, θ)

            # compute step
            Δ = jac \ res

            # line search the step direction
            α = 1.0

            iter = 0
            while check_variables(z - α * Δ) # backtrack inequalities
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > 50
                    @error "backtracking line search fail"
                    flag = false
                    return z, false
                end
            end

            while norm(r(z - α * Δ, θ))^2.0 >= (1.0 - 0.001 * α) * norm(res)^2.0
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > 50
                    @error "line search fail"
                    flag = false
                    return z, false
                end
            end

            # update
            z .-= α * Δ
        end

        μ = 0.1 * μ
        # println("μ: $μ")
    end

    # Δz = -1.0 * (Rz(z, θ)' * Rz(z, θ) + 1.0e-6 * I) \ (Rz(z, θ)' * Rθ(z, θ)) # damped least-squares direction
    Δz = -1.0 * Rz(z, θ) \ Rθ(z, θ)

    q3 = view(z, 1:3)
    n = z[4]
    b = view(z, 6:8)
    Δq1 = view(Δz, 1:3, 1:3)
    Δq2 = view(Δz, 1:3, 3 .+ (1:3))
    Δu1 = view(Δz, 1:3, 6 .+ (1:3))

    return q3, n, b, Δq1, Δq2, Δu1, flag
end

"""
    step particle system
        solves 1-step feasibility problem
"""
function step_slack(q1, q2, u1, h;
    tol = 1.0e-8, max_iter = 100)
    # 1-step optimization problem:
    #     find z
    #     s.t. r(z) = 0
    #
    # z = (q, n, sϕ, sb, db)
    #     s are slack variables for convenience
    #     sb[2:3] is the friction force

    # initialize
    z = 1.0 * ones(11)
    z[1:3] = copy(q2)

    # initialize soc variables
    sb, _ = κ_so(z[6:8])
    sb[1] += 1.0
    db, _ = κ_so(z[9:11])
    db[1] += 1.0

    z[6:8] = sb
    z[9:11] = db

    # μ = 1.0 # barrier parameter
    flag = false

    θ = [q1; q2; u1]

    e = [1.0; 0.0; 0.0]

    function r(z, θ)
        # system variables
        q3 = view(z, 1:3)
        n = z[4]
        sϕ = z[5]
        sb = view(z, 6:8)
        db = view(z, 9:11)

        q1 = view(θ, 1:3)
        q2 = view(θ, 4:6)
        u1 = view(θ, 7:9)
        μ = θ[10]

        λ = [view(sb, 2:3); n] # contact forces
        ϕ = signed_distance(model, q3) # signed-distance function
        vT = (view(q3, 1:2) - view(q2, 1:2)) ./ h

        # @show u1
        # action optimality conditions
        [dynamics(model, q1, q2, q3, u1, λ, h);
         sϕ - ϕ;
         n * sϕ - μ;

         # maximum dissipation optimality conditions
         vT - view(db, 2:3);
         sb[1] - model.friction_coeff * n;
         cone_product(db, sb) - μ * e]
    end

    # Jacobian
    function Rz(z, θ)
        # differentiate r
        _r(w) = r(w, θ)
        _R = ForwardDiff.jacobian(_r, z)
        return _R
    end

    function Rθ(z, θ)
        # differentiate r
        _r(w) = r(z, w)
        _R = ForwardDiff.jacobian(_r, θ)

        return _R
    end

    function check_variables(z)
        # system variables
        q3 = view(z, 1:3)
        n = z[4]
        sϕ = z[5]
        sb = view(z, 6:8)
        db = view(z, 9:11)

        if n <= 0.0
            return true
        end

        if sϕ <= 0.0
            return true
        end

        if !κ_so(sb)[2]
            # println("bs not in cone")
            return true
        end

        if !κ_so(db)[2]
            # println("bz not in cone")
            return true
        end

        return false
    end

    for i = 1:max_iter
        # compute residual, residual Jacobian
        res = r(z, θ)
        if norm(res) < tol
            # println("   iter ($i) - norm: $(norm(res))")
            # return z, true
            flag = true
            continue
        end

        jac = Rz(z, θ)

        # compute step
        Δ = jac \ res

        # line search the step direction
        α = 1.0

        iter = 0
        while check_variables(z - α * Δ) # backtrack inequalities
            α = 0.5 * α
            # println("   α = $α")
            iter += 1
            if iter > 50
                @error "backtracking line search fail"
                flag = false
                return z, false
            end
        end

        while norm(r(z - α * Δ, θ))^2.0 >= (1.0 - 0.001 * α) * norm(res)^2.0
            α = 0.5 * α
            # println("   α = $α")
            iter += 1
            if iter > 50
                @error "line search fail"
                flag = false
                return z, false
            end
        end

        # update
        z .-= α * Δ
    end

    # Δz = -1.0 * (Rz(z, θ)' * Rz(z, θ) + 1.0e-6 * I) \ (Rz(z, θ)' * Rθ(z, θ)) # damped least-squares direction
    Δz = -1.0 * Rz(z, θ) \ Rθ(z, θ)

    q3 = view(z, 1:3)
    n = z[4]
    b = view(z, 6:8)
    Δq1 = view(Δz, 1:3, 1:3)
    Δq2 = view(Δz, 1:3, 3 .+ (1:3))
    Δu1 = view(Δz, 1:3, 6 .+ (1:4))

    return q3, n, b, Δq1, Δq2, Δu1, flag
end


text:
module XGBoost

include("xgboost_lib.jl")
export DMatrix, Booster
export xgboost, predict, save, nfold_cv, slice, get_info, set_info, dump_model, importance

end # module


text:
using WordTokenizers: rulebased_split_sentences, nltk_word_tokenize

"""
abstract type for tokenizers.

Each tokenizer is link with a tokenization (by
 defining `tokenization(::Tokenizer) = Tokenization()`).
 The overall framework dispatch on both tokenizer and
 tokenization, but most of the time we only add methods
 for tokenization. This allow further composability and
 can interfere the tokenization process with given
 tokenizer.
"""
abstract type AbstractTokenizer end

"""
abstract type for tokenization.

The tokenization procedure is separate into multiple
 `TokenStages` and recursive calls of `splitting`, `wrap`,
 and `tokenize`. `splitting` break string into substrings,
 `wrap` mark the substrings with new `TokenStages`, and
 `tokenize` is responsible for the tokenization.
"""
abstract type AbstractTokenization end

"""
abstract type for type that wrap input into specific stage for control tokenization.

There are six builtin stages in TextEncodeBase (all abstract XStage <: TokenStages):

    1. Document <: DocumentStage: the input string is a full document,
     and thus need to be splitted into multiple sentence.
    2. Sentence <: SentenceStage: the input string is a full string,
     and thus need to be splitted into multiple part (SubSentence/Word/Token).
    3. SubSentence <: SubSentenceStage: special wrapper for case where the tokenizer
     does not directly break sentence all into words/tokens and these pieces contain
     multiple words/tokens, but you need the information that they are not full sentence.
    4. Word <: WordStage: the input string is a single word.
    5. SubWord <: SubWordStage: similar to SubSentence, but for word.
    6. Token <: TokenStage: the final piece of the tokenization process.
     Generally, it's used to specify the end of this piece and should
     never be splitted.

Each wrapper have two field: `x` for the input, `meta` for extra information (`nothing` if not provided).
"""
abstract type TokenStages end
abstract type DocumentStage    <: TokenStages end
abstract type SentenceStage    <: TokenStages end
abstract type SubSentenceStage <: TokenStages end
abstract type WordStage        <: TokenStages end
abstract type SubWordStage     <: TokenStages end
abstract type TokenStage       <: TokenStages end

struct Document{T, M}    <: DocumentStage    ; x::T; meta::M; end
struct Sentence{T, M}    <: SentenceStage    ; x::T; meta::M; end
struct SubSentence{T, M} <: SubSentenceStage ; x::T; meta::M; end
struct Word{T, M}        <: WordStage        ; x::T; meta::M; end
struct SubWord{T, M}     <: SubWordStage     ; x::T; meta::M; end
struct Token{T, M}       <: TokenStage       ; x::T; meta::M; end

Document(x) = Document(x, nothing)
Sentence(x) = Sentence(x, nothing)
SubSentence(x) = SubSentence(x, nothing)
Word(x) = Word(x, nothing)
SubWord(x) = SubWord(x, nothing)
Token(x) = Token(x, nothing)

getvalue(x::TokenStages) = x.x
getmeta(x::TokenStages) = x.meta
hasmeta(x::TokenStages) = !isnothing(getmeta(x))

setmeta(x::Document, meta) = Document(x.x, meta)
setmeta(x::Sentence, meta) = Sentence(x.x, meta)
setmeta(x::SubSentence, meta) = SubSentence(x.x, meta)
setmeta(x::Word, meta) = Word(x.x, meta)
setmeta(x::SubWord, meta) = SubWord(x.x, meta)
setmeta(x::Token, meta) = Token(x.x, meta)

setvalue(x::Document, y) = Document(y, x.meta)
setvalue(x::Sentence, y) = Sentence(y, x.meta)
setvalue(x::SubSentence, y) = SubSentence(y, x.meta)
setvalue(x::Word, y) = Word(y, x.meta)
setvalue(x::SubWord, y) = SubWord(y, x.meta)
setvalue(x::Token, y) = Token(y, x.meta)

updatemeta(::Nothing, meta) = meta
updatemeta(a::NamedTuple, meta::NamedTuple) = merge(a, meta)

updatevalue(f, x::TokenStages) = setvalue(x, f(getvalue(x)))
updatemeta(x::TokenStages, meta) = setmeta(x, updatemeta(getmeta(x), meta))

function Base.show(io::IO, t::TokenStages)
    print(io, typeof(t).name.name)
    vs = filter(!isnothing, ntuple(i->getfield(t, i), fieldcount(typeof(t))))
    if length(vs) == 1
        print(io, '(')
        show(io, vs[1])
        print(io, ')')
    else
        print(io, vs)
    end
end

const ParentStages = Union{Nothing, TokenStages}

"""
splittability trait

The splittability trait decide whether the given combination (tokenizer x tokenization x stage) is
 splittable or not (`Splittable` or `UnSplittable`). For example, `DefaultTokenization` and `SentenceStage`
 is splittable (i.e. `splittability(::DefaultTokenization, ::SentenceStage) = Splittable()`). The splittability
 change the behavior of `tokenize`: if it's splittable, `tokenize` will try to call `splitting` on the input,
 `wrap` each splitting result and recurse. Otherwise, it will directly call `wrap` and then recurse into `tokenize`.
"""
abstract type Splittability end
struct Splittable <: Splittability end
struct UnSplittable <: Splittability end

"""
    splittability(args...)

Return the splittability (`Splittable`/`UnSplittable`) of given argument combination.
 Overload to make a `TokenStages` splittable.
"""
function splittability end

"""
    splittable(args...)

Return `true` if the splittability of given argument combination is `Splittable()`.
"""
splittable(args...) = splittable(splittability(args...))
splittable(::Splittable)   = true
splittable(::UnSplittable) = false

splitting(::typeof(splittability), args...) = splitting(splittability(args...), args...)
splitting(::Splittable, args...) = splitting(args...)
splitting(::UnSplittable, args...) = error("Argument is unsplittable: ", args)

# dispatch: tokenizer -> parent stage -> tokenization -> token stage
let ATR = AbstractTokenizer, AT = AbstractTokenization
    # splittability: overload to make specific combination splittable
    global @inline splittability(tkr::ATR, t::AT, x::TokenStages) = splittability(tkr, nothing, t, x)
    global @inline splittability(tkr::ATR, s::ParentStages, t::AT, x::TokenStages) = splittability(s, t, x)
    global @inline splittability(::ParentStages, t::AT, x::TokenStages) = splittability(t, x)
    global @inline splittability(t::AT, x::TokenStages) = UnSplittable()
    # after `Splittable`, define how to split it
    global @inline splitting(::ATR, p::ParentStages, t::AT, x::TokenStages) = splitting(p, t, x)
    global @inline splitting(p::ParentStages, t::AT, x::TokenStages) = splitting(t, x)
    # a callback for `splitting`, `x` is the result of `splitting(::ATR, ::ParentStages,  ::TokenStages)`
    global @inline splitting(::ATR, p::ParentStages, t::AT, s::TokenStages, x) = splitting(p, t, s, x)
    global @inline splitting(p::ParentStages, t::AT, s::TokenStages, x) = splitting(t, s, x)
    global @inline splitting(::AT, ::TokenStages, x) = x
    # splittable (4-arg): wrap the splitting result into specific `TokenStages`, e.g. "word" => Word("word")
    global @inline wrap(::ATR, p::ParentStages, t::AT, s::TokenStages, x) = wrap(p, t, s, x)
    global @inline wrap(p::ParentStages, t::AT, s::TokenStages, x) = wrap(t, s, x)
    global @inline wrap(::AT, ::TokenStages, x::TokenStages) = x # already wrapped
    # unsplittable (3-arg): transform the input into next `TokenStages`, e.g. Word("word") => Token("word")
    global @inline wrap(::ATR, p::ParentStages, t::AT, x::TokenStages)    = wrap(p, t, x)
    global @inline wrap(p::ParentStages, t::AT, s::TokenStages)    = wrap(t, s)
    # the outer-most api, splitting input and recursively tokenize the result. ignore if input is empty
    global @inline tokenize(tkr::ATR, t::AT, x::TokenStages) = tokenize(tkr, nothing, t, x)
    global @inline tokenize(tkr::ATR, s::ParentStages, t::AT, x::TokenStages) = tokenize_procedure(tkr, s, t, x)
    global @inline tokenize(tkr::ATR, s::ParentStages, t::AT, x::TokenStage) = isempty(getvalue(x)) ? TokenStage[] : TokenStage[wrap(tkr, s, t, x)]
end

"""
    tokenization_procedure(tokenizer, tokenizaton, stage)

The procedure of tokenization (`splitting` + `wrap` + `tokenize`).
"""
@inline tokenize_procedure(tkr, t, x) = tokenize_procedure(tkr, t, nothing, x)
@inline tokenize_procedure(tkr, s, t, x) = tokenize_procedure!(append!, splittability, TokenStage[], tkr, s, t, x)
@inline tokenize_procedure!(op, v, tkr, s, t, x) = tokenize_procedure!(op, splittability, v, tkr, s, t, x)
@inline tokenize_procedure!(op, ::typeof(splittability), v, tkr, s, t, x) = tokenize_procedure!(op, splittability(tkr, s, t, x), v, tkr, s, t, x)

function tokenize_procedure!(op, ::Splittable, v, tkr, s, t, x)
    isempty(getvalue(x)) && return v
    for sp in splitting(tkr, s, t, x, splitting(splittability, tkr, s, t, x))
        v1 = tokenize(tkr, x, t, wrap(tkr, s, t, x, sp))
        op(v, v1)
    end
    return v
end

function tokenize_procedure!(op, ::UnSplittable, v, tkr, s, t, x)
    isempty(getvalue(x)) && return v
    op(v, tokenize(tkr, x, t, wrap(tkr, s, t, x)))
    return v
end

"""
    splitting(t::AbstractTokenization, x::TokenStages)

Split `x` given its tokenization stage. For example,
 the default behavior of a document stage is splitting into
 sentences (with `WordTokenizers.split_sentences`).

Overload this method for custom tokenization.
"""
function splitting end

"""
    wrap(t::AbstractTokenization, s::TokenStages, x)

Mark the tokenization stage of `x`, which is part of the splitting result of `s`.
 For example, if we are doing simple whitespace tokenization and at the sentence stage,
 then `x` is just single word of `s` and thus return `Word(x)` (or `Token(x)`).
 Skip if `x` is already a `TokenStages`. (this method only apply to splittable stages)

Overload this method to control the tokenization process.
"""
function wrap end

@eval $((@macroexpand @doc """
    wrap(t::AbstractTokenization, x::TokenStages)

A handler for unsplittable stages (token/word/...).

Overload this method for custom transform.
"""
function wrap(t::AbstractTokenization, x::TokenStages) end
).args[2])

# abstract type for convenience

abstract type BaseTokenization <: AbstractTokenization end

struct DefaultTokenization <: BaseTokenization end

splittability(::BaseTokenization, x::Union{DocumentStage, SentenceStage, SubSentenceStage}) = Splittable()
splittability(::BaseTokenization, x::Union{WordStage, SubWordStage}) = UnSplittable()

splitting(::BaseTokenization, d::DocumentStage)    = rulebased_split_sentences(getvalue(d))
splitting(::BaseTokenization, s::SentenceStage)    = nltk_word_tokenize(getvalue(s))
splitting(::BaseTokenization, s::SubSentenceStage) = nltk_word_tokenize(getvalue(s))

wrap(::BaseTokenization, d::DocumentStage, x) = Sentence(x, getmeta(d))
wrap(::BaseTokenization, s::SentenceStage, x) = Word(x, getmeta(s))
wrap(::BaseTokenization, s::SubSentenceStage, x) = Word(x, getmeta(s))
wrap(::BaseTokenization, w::WordStage, x) = SubWord(x, getmeta(w))

wrap(::BaseTokenization, w::WordStage)    = Token(getvalue(w), getmeta(w))
wrap(::BaseTokenization, w::SubWordStage) = Token(getvalue(w), getmeta(w))
wrap(::BaseTokenization, t::TokenStage)   = t

abstract type WrappedTokenization{T<:AbstractTokenization} <: AbstractTokenization end

base(t::WrappedTokenization) = t.base

splittability(p::Union{Nothing, TokenStages}, t::WrappedTokenization, x::TokenStages) = splittability(p, base(t), x)
splitting(p::Union{Nothing, TokenStages}, t::WrappedTokenization, x::TokenStages)    = splitting(p, base(t), x)
splitting(p::Union{Nothing, TokenStages}, t::WrappedTokenization, s::TokenStages, x) = splitting(p, base(t), s, x)

wrap(p::Union{Nothing, TokenStages}, t::WrappedTokenization, s::TokenStages, x) = wrap(p, base(t), s, x)
wrap(p::Union{Nothing, TokenStages}, t::WrappedTokenization, s::TokenStages)    = wrap(p, base(t), s)

# tokenizer api

"""
    tokenization(::AbstractTokenizer) :: AbstractTokenization

Return the tokenization object of given tokenizer.
"""
tokenization(::AbstractTokenizer) = DefaultTokenization()

"""
    preprocess(tkr::AbstractTokenizer, x)

Preprocess the input `x`. This is only called during `tkr(x)`.
"""
preprocess(t::AbstractTokenizer, x::TokenStages) = updatevalue(Base.Fix1(preprocess, t), x)
preprocess(t::AbstractTokenizer, x) = x

(t::AbstractTokenizer)(x::TS) where {TS <: TokenStages} = tokenize(t, nothing, tokenization(t), preprocess(t, x))

function Base.show(io::IO, t::AbstractTokenizer)
    T = typeof(t)
    n = fieldcount(T)
    print(io, nameof(T))
    if n != 0
        base = hasfield(T, :tokenization)
        basei = base ? findfirst(==(:tokenization), fieldnames(T)) : 0
        print(io,  '(')
        base && show(io, t.tokenization)
        for i = 1:n
            i == basei && continue
            (!base && i == 1) || print(io, ", ")
            print(io, fieldname(T, i))
            print(io, " = ")
            show(IOContext(io, :limit=>true), getfield(t, i))
        end
        print(io, ')')
    end
end

function Base.show(io::IO, t::AbstractTokenization)
    T = typeof(t)
    n = fieldcount(T)
    print(io, nameof(T))
    if n != 0
        base = hasfield(T, :base)
        basei = base ? findfirst(==(:base), fieldnames(T)) : 0
        print(io,  '(')
        base && show(io, t.base)
        for i = 1:n
            i == basei && continue
            (!base && i == 1) || print(io, ", ")
            print(io, fieldname(T, i))
            print(io, " = ")
            show(IOContext(io, :limit=>true), getfield(t, i))
        end
        print(io, ')')
    end
end

Base.show(io::IO, ::DefaultTokenization) = print(io, :default)