text:
# File defining the preconditioner types

type SparsifyingPreconditioner
    Msp::SparseMatrixCSC{Complex{Float64},Int64}
    As::SparseMatrixCSC{Complex{Float64},Int64} # tranposition matrix
    MspInv
    solverType::String
    function SparsifyingPreconditioner(Msp, As; solverType::String="UMPFACK")
        tic();

        if solverType == "UMPFACK"
            MspInv = lufact(Msp)

        elseif solverType == "MKLPardiso"
            # using MKLPardiso from Julia Sparse (only shared memory)
            println("factorizing the matrix using MKL Pardiso")
            MspInv = MKLPardisoSolver();
            set_nprocs!(MspInv, 4)
            #setting the type of the matrix
            set_matrixtype!(MspInv,3)
            # setting we are using a transpose
            set_iparm!(MspInv,12,2)
            # setting the factoriation phase
            set_phase!(MspInv, 12)
            X = zeros(Complex128, size(Msp)[1],1)
            # factorizing the matrix
            pardiso(MspInv,X, Msp,X)
            # setting phase and parameters to solve and transposing the matrix
            # this needs to be done given the different C and Fortran convention
            # used by Pardiso (C convention) and Julia (Fortran Convention)
            set_phase!(MspInv, 33)
            set_iparm!(MspInv,12,2)
        end

        println("time for the factorization was ", toc() )
        new(Msp,As, MspInv,solverType) #
    end
end

# Encapsulation of the preconditioner in order to use preconditioned GMRES
import Base.\
import Base.A_ldiv_B!

function \(M::SparsifyingPreconditioner, b::Array{Complex128,1})
    # we apply the Sparsifying preconditioner
    rhs = (M.As*b)
    if M.solverType == "UMPFACK"
        return M.MspInv\rhs
    elseif M.solverType == "MKLPardiso"
        set_phase!(M.MspInv, 33)
        u = zeros(Complex128,length(rhs),1)
        pardiso(M.MspInv, u, M.Msp, rhs)
    end
end

# function A_ldiv_B!(M::SparsifyingPreconditioner, v)
# # we apply the Sparsifying preconditioner
#     rhs = (M.As*v)
#     if M.solverType == "UMPFACK"
#         v =  M.MspInv\rhs
#     elseif M.solverType == "MKLPardiso"
#         set_phase!(M.MspInv, 33)
#         v = zeros(Complex128,length(rhs),1)
#         pardiso(M.MspInv, v, M.Msp, rhs)
#     end
# end

function A_ldiv_B!(M::SparsifyingPreconditioner, v::SubArray{Complex{Float64},1,Array{Complex{Float64},2}})
    # we apply the Sparsifying preconditioner
    rhs = (M.As*v)
    if M.solverType == "UMPFACK"
        v[:] =  M.MspInv\rhs
    elseif M.solverType == "MKLPardiso"
        set_phase!(M.MspInv, 33)
        #v = zeros(Complex128,length(rhs),1)
        pardiso(M.MspInv, v, M.Msp, rhs)
    end
end

# function A_ldiv_B!(u::SubArray{Complex{Float64},1,Array{Complex{Float64},2}},
#                    M::SparsifyingPreconditioner,
#                    v::SubArray{Complex{Float64},1,Array{Complex{Float64},2}})
# # we apply the Sparsifying preconditioner
#     @assert length(u) == length(v)
#     rhs = (M.As*v)
#     if M.solverType == "UMPFACK"
#         u =  M.MspInv\rhs
#     elseif M.solverType == "MKLPardiso"
#         set_phase!(M.MspInv, 33)
#         #u = zeros(Complex128,length(rhs),1)
#         pardiso(M.MspInv, u, M.Msp, rhs)
#     end
# end

# function A_ldiv_B!(M::SparsifyingPreconditioner, v)
#     # we apply the Sparsifying preconditioner
#     rhs = (M.As*v)
#     if M.solverType == "UMPFACK"
#         v =  M.MspInv\rhs
#     elseif M.solverType == "MKLPardiso"
#         set_phase!(M.MspInv, 33)
#         v = zeros(Complex128,length(rhs),1)
#         pardiso(M.MspInv, v, M.Msp, rhs)
#     end
#     println("hey!")
# end

function A_ldiv_B!(u,
                   M::SparsifyingPreconditioner,
                   v)
# we apply the Sparsifying preconditioner
    @assert length(u) == length(v)
    rhs = (M.As*v)
    if M.solverType == "UMPFACK"
        u =  M.MspInv\rhs
    elseif M.solverType == "MKLPardiso"
        set_phase!(M.MspInv, 33)
        #u = zeros(Complex128,length(rhs),1)
        pardiso(M.MspInv, u, M.Msp, rhs)
    end
end

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
using Flux: trainable, IdSet

export namedparams, namedparams!

namedparams!(p, x::AbstractArray{<:Number}, prefix, seen = IdSet()) = p[x] = strip(prefix, ['/'])

function namedparams!(p, x, prefix, seen = IdSet())
    x in seen && return
    push!(seen, x)
    for (name, child) in pairs(trainable(x))
        namedparams!(p, child, prefix * string(name) * "/", seen)
    end
end

function namedparams(m...)
    ps = IdDict{Any, String}()
    namedparams!(ps, m, "")
    return ps
end

ext: jl


text:
# --- convolution

# http://developer.amd.com/resources/documentation-articles/articles-whitepapers/tiled-convolution-fast-image-filtering/
# has a detailed argument about why and how to do this better on a
# real gpu (although this case is simpler as it fits in a single tile)

"""convolve a 3x3 cross-shaped kernel."""
function convolve(input, border, kernel, output, n)
    @forall i j n begin
        a = 0
        # cannot use forneighbours because we need di, dj
        for (di, dj) in ((0, 1), (0, -1), (1, 0), (-1, 0))
            ii, jj = i+di, j+dj
            if ii > 0 && ii <= n && jj > 0 && jj <= n
                a = a + kernel[di+2, dj+2] * input[ii, jj]
            else
                a = a + kernel[di+2, dj+2] * border
            end
            output[i, j] = a
        end
    end
end

function blank!{N}(data, p::Position{N})
    @forall i j N begin
        if point(p, i, j) != empty
            data[i, j] = 0
        end
    end
end


text:
using BinaryBuilder

# zlib version
name = "Zlib"
version = v"1.2.11"

# Collection of sources required to build zlib
sources = [
    ArchiveSource("https://zlib.net/zlib-$(version).tar.gz",
                  "c3e5e9fdd5004dcb542feda5ee4f0ff0744628baf8ed2dd5d66f8ca1197cb1a1"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/zlib-*
mkdir build && cd build

# We use `-DUNIX=true` to ensure that it is always named `libz` instead of `libzlib` or something ridiculous like that.
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DUNIX=true \
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON \
    ..
make install -j${nproc}
install_license ../README
"""

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libz", :libz),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Note: we explicitly lie about this because we don't have the new
# versioning APIs worked out in BB yet.
version = v"1.2.12"
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using PyTorchTraining
using Test

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


text:
input = [1, 1, 2, 2]

push!(input, input[1])

println(input)

lastVal = 0
sum = 0
for val in input

    if val == lastVal
        sum += val
        lastVal = val
    end
end
println(sum)


text:
# This is used on Julia version that have the Base.Ryu module.

using Base.Ryu

function plain_precision_heuristic(xs::AbstractArray{<:AbstractFloat})
    ys = filter(isfinite, xs)
    e10max = -(e10min = typemax(Int))
    for y in ys
        if isapprox(y, 0, atol=1e-16)
            continue
        end
        _, e10 = Ryu.reduce_shortest(convert(Float32, y))
        e10min = min(e10min, e10)
        e10max = max(e10max, e10)
    end
    return precision = min(-e10min, -e10max+16)
end

# Print a floating point number at fixed precision. Pretty much equivalent to
# @sprintf("%0.$(precision)f", x), without the macro issues.
function format_fixed(x::AbstractFloat, precision::Integer)
    if x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    return Ryu.writefixed(x, precision)
end

# Print a floating point number in scientific notation at fixed precision. Sort of equivalent
# to @sprintf("%0.$(precision)e", x), but prettier printing.
function format_fixed_scientific(x::AbstractFloat, precision::Integer,
                                 engineering::Bool)
    if iszero(x)
        return "0"
    elseif isinf(x)
        return signbit(x) ? "-∞" : "∞"
    elseif isnan(x)
        return "NaN"
    end

    if engineering
        # precision applies before convertic to engineering format
        # that is, number of digits would be the same as the with non-engineering format
        base_digits, power = get_engineering_string(x, precision)
    else
        e_format_number = Ryu.writeexp(x, precision)
        base_digits, power = split(e_format_number, 'e')
    end

    buf = IOBuffer()

    print(buf, base_digits)
    print(buf, "×10")

    if power[1] == '-'
        print(buf, '⁻')
    end
    leading_index = findfirst(c -> '1' <= c <= '9', power)

    if leading_index === nothing
        print(buf, superscript_numerals[1])
        return String(take!(buf))
    end

    for digit in power[leading_index:end]
        if digit == '-'
            print(buf, '⁻')
        elseif '0' <= digit <= '9'
            print(buf, superscript_numerals[digit - '0' + 1])
        end

    end

    return String(take!(buf))
end

function get_engineering_string(x::AbstractFloat, precision::Integer)
    e_format_number = Ryu.writeexp(x, precision)
    base_digits, power = split(e_format_number, 'e')

    int_power = parse(Int, power)
    positive = int_power >= 0

    # round the power to the nearest multiple of 3
    # positive power -> move the "." to the right by mode, round the power to the higher power
    # negative power -> move the "." to the right by mode, round the power to the lower power
    # ex:
    # 1.2334e5 = 123.334e3
    # 1.2334-5 = 12.3334e-6

    if positive
        indices_to_move = int_power % 3
    else
        indices_to_move = 3 - abs(int_power) % 3
    end

    buf = IOBuffer()
    negative_base_compensation = base_digits[1] == '-' ? 1 : 0
    for i in eachindex(base_digits)
        if base_digits[i] != '.'
            print(buf, base_digits[i])
        end
        if i == 2 + indices_to_move + negative_base_compensation
            print(buf, '.')
        end
    end

    return String(take!(buf)), string(int_power - indices_to_move)
end


text:
function StencilOpState(;
    failOp::vk.VkStencilOp = vk.VK_STENCIL_OP_KEEP,
    passOp::vk.VkStencilOp = vk.VK_STENCIL_OP_KEEP,
    depthFailOp::vk.VkStencilOp = vk.VK_STENCIL_OP_KEEP,
    compareOp::vk.VkCompareOp = vk.VK_COMPARE_OP_NEVER,
    compareMask::UInt32 = UInt32(0),
    writeMask::UInt32 = UInt32(0),
    reference::UInt32 = UInt32(0)
)

    vk.VkStencilOpState(
        failOp, #::VkStencilOp
        passOp, #::VkStencilOp
        depthFailOp, #::VkStencilOp
        compareOp, #::VkCompareOp
        compareMask, #::UInt32
        writeMask, #::UInt32
        reference #::UInt32
    )
end

struct PipelineDepthStencilStateCreateInfo
    mHandleRef::Ref{vk.VkPipelineDepthStencilStateCreateInfo}

    function PipelineDepthStencilStateCreateInfo(;
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineDepthStencilStateCreateFlags = vk.VkFlags(0),
        depthTestEnable::vk.VkBool32 = VkExt.VK_FALSE,
        depthWriteEnable::vk.VkBool32 = VkExt.VK_FALSE,
        depthCompareOp::vk.VkCompareOp = VK_COMPARE_OP_NEVER,
        depthBoundsTestEnable::vk.VkBool32 = VkExt.VK_FALSE,
        stencilTestEnable::vk.VkBool32 = VkExt.VK_FALSE,
        front::vk.VkStencilOpState = StencilOpState(),
        back::vk.VkStencilOpState = StencilOpState(),
        minDepthBounds::Cfloat = Cfloat(0.0),
        maxDepthBounds::Cfloat = Cfloat(0.0)
    )

        this = new(Ref(vk.VkPipelineDepthStencilStateCreateInfo(
            vk.VK_STRUCTURE_TYPE_PIPELINE_DEPTH_STENCIL_STATE_CREATE_INFO, #sType::VkStructureType
            pNext, #::Ptr{Cvoid}
            flags, #::VkPipelineDepthStencilStateCreateFlags
            depthTestEnable, #::VkBool32
            depthWriteEnable, #::VkBool32
            depthCompareOp, #::VkCompareOp
            depthBoundsTestEnable, #::VkBool32
            stencilTestEnable, #::VkBool32
            front, #::VkStencilOpState
            back, #::VkStencilOpState
            minDepthBounds, #::Cfloat
            maxDepthBounds #::Cfloat
        )))
    end
end

function handleRef(this::PipelineDepthStencilStateCreateInfo)::Ref{vk.VkPipelineDepthStencilStateCreateInfo}
    this.mHandleRef
end


text:
using DimensionalData, Test

using DimensionalData: val, basetypeof, slicedims, dims2indices, formatdims, hasdim, swapdim,
      @dim, reducedims, dimnum, X, Y, Z, Time, Forward

dimz = (X(), Y())

@testset "permutedims" begin
    @test permutedims((Y(1:2), X(1)), dimz) == (X(1), Y(1:2))
    @test permutedims((X(1),), dimz) == (X(1), nothing)
    @test permutedims((Y(), X()), dimz) == (X(:), Y(:))
    @test permutedims([Y(), X()], dimz) == (X(:), Y(:))
    @test permutedims((Y, X),     dimz) == (X(:), Y(:))
    @test permutedims([Y, X],     dimz) == (X(:), Y(:))
    @test permutedims(dimz, (Y(), X())) == (Y(:), X(:))
    @test permutedims(dimz, [Y(), X()]) == (Y(:), X(:))
    @test permutedims(dimz, (Y, X)    ) == (Y(:), X(:))
    @test permutedims(dimz, [Y, X]    ) == (Y(:), X(:))
end

a = [1 2 3; 4 5 6]
da = DimensionalArray(a, (X((143, 145)), Y((-38, -36))))
dimz = dims(da)

@testset "slicedims" begin
    @test slicedims(dimz, (1:2, 3)) == ((X(LinRange(143,145,2); grid=RegularGrid(span=2.0)),),
                                        (Y(-36.0; grid=RegularGrid(span=1.0)),))
    @test slicedims(dimz, (2:2, :)) == ((X(LinRange(145,145,1); grid=RegularGrid(span=2.0)),
                                         Y(LinRange(-38.0,-36.0,3); grid=RegularGrid(span=1.0))), ())
    @test slicedims((), (1:2, 3)) == ((), ())
end

@testset "dims2indices" begin
    emptyval = Colon()
    @test dims2indices(grid(dimz[1]), dimz[1], Y, Nothing) == Colon()
    @test dims2indices(dimz, (Y(),), emptyval) == (Colon(), Colon())
    @test dims2indices(dimz, (Y(1),), emptyval) == (Colon(), 1)
    # Time is just ignored if it's not in dims. Should this be an error?
    @test dims2indices(dimz, (Time(4), X(2))) == (2, Colon())
    @test dims2indices(dimz, (Y(2), X(3:7)), emptyval) == (3:7, 2)
    @test dims2indices(dimz, (X(2), Y([1, 3, 4])), emptyval) == (2, [1, 3, 4])
    @test dims2indices(da, (X(2), Y([1, 3, 4])), emptyval) == (2, [1, 3, 4])
    emptyval=()
    @test dims2indices(dimz, (Y,), emptyval) == ((), Colon())
    @test dims2indices(dimz, (Y, X), emptyval) == (Colon(), Colon())
    @test dims2indices(da, X, emptyval) == (Colon(), ())
end

@testset "dimnum" begin
    @test dimnum(da, X) == 1
    @test dimnum(da, Y()) == 2
    @test dimnum(da, (Y, X())) == (2, 1)
    @test_throws ArgumentError dimnum(da, Time) == (2, 1)
end

@testset "reducedims" begin
    @test reducedims((X(3:4), Y(1:5)), (X, Y)) == (X(3), Y(1))
end

@testset "dims" begin
    @test typeof(dims(da, X)) <: X
    @test typeof(dims(dims(da), Y)) <: Y
    @test dims(da, ()) == ()
    @test_throws ArgumentError dims(da, Time)
    x = dims(da, X)
    @test dims(x) == x
end

@testset "hasdims" begin
    @test hasdim(da, X) == true
    @test hasdim(da, Time) == false
    @test hasdim(dims(da), Y) == true
    @test hasdim(dims(da), (X, Y)) == (true, true)
    @test hasdim(dims(da), (X, Time)) == (true, false)
end

@testset "swapdim" begin
    A = swapdim(da, X(LinRange(150,152,2)))
    @test val(dims(dims(A), X())) == LinRange(150,152,2)
end


text:
# from https://github.com/FluxML/Zygote.jl/blob/master/test/gradcheck.jl

function ngradient(f, xs::AbstractArray...)
  grads = zero.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

function gradcheck(f, xs...)
    n_grads = ngradient(f, xs...)
    y_grads = Yota._grad(f, xs...)[2] |> collect
    all(isapprox.(n_grads, y_grads, rtol = 1e-5, atol = 1e-5))
end

# gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
# gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)


text:
import ..update!, ..select_action, ..reset!

using ..Buffers:Experience

using StaticArrays
using Random
using LinearAlgebra:Diagonal, UniformScaling

using Ramnet.Models.AltDiscriminators:AdaptiveFunctionalDiscriminator
using Ramnet.Optimizers
import Ramnet.Loss

struct AdaptiveBinaryPolicyLoss <: Loss.AbstractLoss end

# Policy approximator with only two admissible actions: {-1, 1}.
struct AdaptiveBinaryPolicy{OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    opt::AdaptiveOptimizer{1,AdaptiveBinaryPolicyLoss}
    f::AdaptiveFunctionalDiscriminator{1,T,E}
    rng::MersenneTwister
end

function AdaptiveBinaryPolicy(::Type{O}, n::Int, encoder::E; λ::Float64=1.0, μ::Float64=0.1, η::Float64=0.1, epochs::Int=1, partitioner::Symbol=:uniform_random, seed::Union{Nothing,Int}=nothing) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    AdaptiveBinaryPolicy{OS,T,O,E}(
        AdaptiveOptimizer{1}(AdaptiveBinaryPolicyLoss(); λ, meta_rate=μ, initial_rate=η),
        AdaptiveFunctionalDiscriminator{1}(OS, n, λ, encoder, partitioner; seed),
        isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)
    )
end

# Compute π(action | observation)
function probability(action::Int, pred)
    # @show pred
    # @show typeof(pred)
    1 / (1 + exp(sign(action) * first(pred)))
end

function Loss.grad(::AdaptiveBinaryPolicyLoss, action, pred)
    # @show pred
    # @show typeof(pred)
    sign(action) * (1 - probability(action, pred))
end

function Loss.hessian(::AdaptiveBinaryPolicyLoss, action, pred)
    p = probability(action, pred)
    return -1 * p * (1 - p)
end

# TODO: Validate action (must be ether -1 or +1)
function update!(p::AdaptiveBinaryPolicy{OS,T,O,E}, observation::O, action::Int, G) where {OS,T,O,E}
    gradient, trace_update = Optimizers.grad(p.opt, p.f, observation, action)
    rate = SizedVector{1}(learning_rate(p.opt))

    train!(p.f, observation, G * rate .* gradient, rate .* SizedVector{1}(trace_update))
    # p.opt.loss.update_count = p.opt.loss.update_count == p.opt.loss.cov_decay ? p.opt.loss.cov_decay : p.opt.loss.update_count + 1

    return nothing
end

# function update!(p::AdaptiveBinaryPolicy{OS,AS,T,O,A,E}, experiences::AbstractVector{Experience{O,A}}, critic) where {OS,AS,T,O,A,E}
#     gradients = Matrix{Float64}(undef, AS, length(experiences))
#     observations = Matrix{T}(undef, OS, length(experiences))

#     for i in eachindex(experiences)

#         observation = experiences[i].observation
#         action = experiences[i].action
#         value = experiences[i].value
#         advantage = value - critic(observation)
#         gradients[:, i] = Optimizers.grad(p.opt, p, observation, action, advantage)
#         observations[:, i] = observation
#     end

#     train!(p.f, observations, SizedMatrix{AS,length(experiences)}(learning_rate(p.opt) * gradients))
# p.update_count = p.update_count == p.cov_decay ? p.cov_decay : p.update_count + 1

#     nothing
# end

function select_action(p::AdaptiveBinaryPolicy{OS,T,O,E}, observation::O) where {OS,T,O,E}
    pred, _ = predict(p.f, observation)
    return rand(p.rng) < probability(1, pred) ? 1 : -1
end

function reset!(p::AdaptiveBinaryPolicy)
    Ramnet.reset!(p.f)

    return nothing
end


text:
@namedargs mutable struct DICparam
  p::Vector{Matrix{Float64}}
  mu::Vector{Matrix{Float64}}
  sig::Vector{Vector{Float64}}
  y::Vector{Matrix{Float64}}
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
info("Testing ", MOD.translate!)
for T in FloatingPointTypes
    A = T[1 2;
          3 4;
          5 6]

    b = T[1;
          2]

    c = T[1;
          2;
          3]

    @test_approx_eq MOD.translate!(copy(A), one(T)) (A .+ one(T))
    @test_approx_eq MOD.translate!(one(T), copy(A)) (A .+ one(T))
    @test_approx_eq MOD.translate!(copy(A), b) (A .+ b')
    @test_approx_eq MOD.translate!(c, copy(A)) (A .+ c)
end

info("Testing ", MOD.perturb!)
for T in FloatingPointTypes
    A = T[1 2 3;
          4 5 6;
          7 8 9]

    B = T[1 2;
          3 4;
          5 6]

    @test_approx_eq MOD.perturb!(copy(A), one(T)) (A + one(T)*I)
    @test_throws DimensionMismatch MOD.perturb!(B, one(T))
end

info("Testing ", MOD.regularize!)
for T in FloatingPointTypes
    A = T[1 2 3;
          4 5 6;
          7 8 9]

    B = T[1 2;
          3 4;
          5 6]

    a = trace(A)/3

    @test_approx_eq MOD.regularize!(copy(A), zero(T), a) A
    @test_approx_eq MOD.regularize!(copy(A), convert(T,0.5), a) convert(T,0.5)*(A + a*I)
    @test_approx_eq MOD.regularize!(copy(A), one(T), a) a*eye(T,3)
    @test_throws DimensionMismatch MOD.regularize!(B, one(T), one(T))
end

A  = [1 2 3;
      4 5 6;
      7 8 9]
AL = [1 2 3;
      2 5 6;
      3 6 9]
AU = [1 4 7;
      4 5 8;
      7 8 9]

info("Testing ", MOD.syml)
for T in FloatingPointTypes
    B = MOD.syml(convert(Array{T},A))
    @test eltype(B) == T
    @test_approx_eq B convert(Array{T}, AL)
end

info("Testing ", MOD.symu)
for T in FloatingPointTypes
    B = MOD.symu(convert(Array{T},A))
    @test eltype(B) == T
    @test_approx_eq B convert(Array{T}, AU)
end

info("Testing ", MOD.components_svd!)
for T in FloatingPointTypes
    X = rand(T,20,5)
    V1, D1 = MOD.components_svd!(copy(X))

    D2, V2 = svd(X)[2:3]

    @test_approx_eq V1 V2
    @test_approx_eq D1 (D2 .^ 2)
end

info("Testing ", MOD.components_eig!)
for T in FloatingPointTypes
    X = rand(T,20,5)
    S = X'X

    V1, D1 = MOD.components_eig!(copy(S))

    D2, V2 = eig(S)

    @test_approx_eq abs(V1) abs(V2[:,end:-1:1])  # Signs may differ
    @test_approx_eq abs(D1) abs(D2[end:-1:1])
end

info("Testing ", MOD.components_geig!)
for T in FloatingPointTypes
    X1 = rand(T,20,5)
    X2 = rand(T,20,5)

    S1 = X1'X1
    S2 = X2'X2

    V1, D1 = MOD.components_geig!(copy(S1), copy(S2))

    D2, V2 = eig(S1, S2)

    @test_approx_eq abs(V1) abs(V2[:,end:-1:1])
    @test_approx_eq abs(D1) abs(D2[end:-1:1])
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    precompile(Tuple{typeof(setindex!), typeof(fileinfos), FileInfo, String})
    precompile(Tuple{typeof(setindex!), typeof(module2files), Vector{String}, Symbol})
    precompile(Tuple{typeof(setindex!), typeof(watched_files), WatchList, String})
    precompile(Tuple{typeof(haskey), typeof(watched_files), String})

    precompile(Tuple{typeof(_watch_package), Base.PkgId})
    precompile(Tuple{typeof(revise_dir_queued), String})
    precompile(Tuple{typeof(run_backend), REPL.REPLBackend})
    precompile(Tuple{typeof(steal_repl_backend), REPL.REPLBackend})

    return nothing
end


text:
function average end

"""
    AvgSym <: CNumber

Symbolic number representing the average over an operator.
See also: [`average`](@ref)
"""
struct AvgSym <: CNumber end

const Average = SymbolicUtils.Term{<:AvgSym}

const sym_average = begin # Symbolic function for averages
    T = SymbolicUtils.FnType{Tuple{QNumber}, AvgSym}
    SymbolicUtils.Sym{T}(:avg)
end

# Type promotion -- average(::QNumber)::Number
SymbolicUtils.promote_symtype(::typeof(sym_average), ::Type{<:QNumber}) = AvgSym

# Direct construction of average symbolic expression
function _average(operator)
    return SymbolicUtils.Term{AvgSym}(sym_average, [operator])
end

function acts_on(s::SymbolicUtils.Symbolic)
    if SymbolicUtils.istree(s)
        f = SymbolicUtils.operation(s)
        if f === sym_average
            return acts_on(SymbolicUtils.arguments(s)[1])
        else
            aon = []
            for arg in SymbolicUtils.arguments(s)
                append!(aon, acts_on(arg))
            end
            unique!(aon)
            sort!(aon)
            return aon
        end
    else
        return Int[]
    end
end

"""
    average(::QNumber)
    average(::QNumber,order)

Compute the average of an operator. If `order` is given, the [`cumulant_expansion`](@ref)
up to that order is computed immediately.
"""
average(op::QSym) = _average(op)
function average(op::QTerm)
    f = SymbolicUtils.operation(op)
    if f===(+) || f===(-) # linearity
        args = map(average, SymbolicUtils.arguments(op))
        return f(args...)
    elseif f === (*)
        # Move constants out of average
        c = op.arg_c
        op_ = QMul(1,op.args_nc)
        return c*_average(op_)
    else
        error("Unknown function $f")
    end
end
average(x::SNuN) = x
average(x,order;kwargs...) = cumulant_expansion(average(x),order;kwargs...)

function undo_average(t)
    if SymbolicUtils.istree(t)
        f = SymbolicUtils.operation(t)
        if f === sym_average
            return SymbolicUtils.arguments(t)[1]
        else
            args = map(undo_average, SymbolicUtils.arguments(t))
            return f(args...)
        end
    else
        return t
    end
end

function undo_average(eq::Symbolics.Equation)
    lhs = undo_average(eq.lhs)
    rhs = undo_average(eq.rhs)
    return Symbolics.Equation(lhs,rhs)
end

## Cumulant expansion
"""
    cumulant_expansion(avg, order::Int)

For an [`average`](@ref) of an operator, expand it in terms
of moments up to `order` neglecting their joint cumulant.

See also: https://en.wikipedia.org/wiki/Cumulant#Joint_cumulants

Examples
=======
```
julia> avg = average(a*b)
⟨a*b⟩

julia> cumulant_expansion(avg,1)
(⟨a⟩*⟨b⟩)

julia> avg = average(a*b*c)
⟨a*b*c⟩

julia> cumulant_expansion(avg,2)
((⟨a*b⟩*⟨c⟩)+(⟨a*c⟩*⟨b⟩)+(⟨a⟩*⟨b*c⟩)+(-2*⟨a⟩*⟨b⟩*⟨c⟩))
```

Optional arguments
=================
*simplify=true: Specify whether the result should be simplified.
*kwargs...: Further keyword arguments being passed to simplification.
"""
function cumulant_expansion(x::SymbolicUtils.Symbolic,order::Integer;simplify=true,kwargs...)
    if SymbolicUtils.istree(x)
        get_order(x) <= order && return x
        f = SymbolicUtils.operation(x)
        if f===sym_average
            op = SymbolicUtils.arguments(x)[1]
            return _cumulant_expansion(op.args_nc, order)
        else
            args = SymbolicUtils.arguments(x)
            cumulants = [cumulant_expansion(arg,order;kwargs...) for arg in args]
            return f(cumulants...)
        end
    else
        return x
    end
end
function cumulant_expansion(avg::SymbolicUtils.Term{<:AvgSym},order::Vector;mix_choice=maximum,kwargs...)
    aon = acts_on(avg)
    order_ = mix_choice(order[get_i(i)] for i in aon)
    return cumulant_expansion(avg,order_;kwargs...)
end
cumulant_expansion(x::Number,order;kwargs...) = x

function cumulant_expansion(x::SymbolicUtils.Symbolic,order;mix_choice=maximum,simplify=true,kwargs...)
    if SymbolicUtils.istree(x)
        f = SymbolicUtils.operation(x)
        args = SymbolicUtils.arguments(x)
        cumulants = [cumulant_expansion(arg,order;simplify=simplify,mix_choice=mix_choice) for arg in args]
        if simplify
            return SymbolicUtils.simplify(f(cumulants...);kwargs...)
        else
            return f(cumulants...)
        end
    else
        return x
    end
end
function cumulant_expansion(de::MeanfieldEquations,order;multithread=false,mix_choice=maximum,kwargs...)
    order==de.order && return de
    eqs = de.equations
    eqs_out = Vector{Symbolics.Equation}(undef, length(eqs))
    if multithread
        Threads.@threads for i=1:length(eqs)
            cr = cumulant_expansion(eqs[i].rhs,order;mix_choice=mix_choice,kwargs...)
            eqs_out[i] = Symbolics.Equation(eqs[i].lhs, cr)
        end
    else
        for i=1:length(eqs)
            cr = cumulant_expansion(eqs[i].rhs,order;mix_choice=mix_choice,kwargs...)
            eqs_out[i] = Symbolics.Equation(eqs[i].lhs, cr)
        end
    end
    return MeanfieldEquations(eqs_out,de.operator_equations,de.states,de.operators,
                            de.hamiltonian,de.jumps,de.jumps_dagger,de.rates,de.iv,de.varmap,
                            order)
end
function cumulant_expansion(de::ScaledMeanfieldEquations,order;multithread=false,mix_choice=maximum,kwargs...)
    order==de.order && return de
    eqs = de.equations
    eqs_out = Vector{Symbolics.Equation}(undef, length(eqs))
    if multithread
        Threads.@threads for i=1:length(eqs)
            cr = cumulant_expansion(eqs[i].rhs,order;mix_choice=mix_choice,kwargs...)
            cr = substitute_redundants(cr, de.scale_aons, de.names)
            eqs_out[i] = Symbolics.Equation(eqs[i].lhs, cr)
        end
    else
        for i=1:length(eqs)
            cr = cumulant_expansion(eqs[i].rhs,order;mix_choice=mix_choice,kwargs...)
            cr = substitute_redundants(cr, de.scale_aons, de.names)
            eqs_out[i] = Symbolics.Equation(eqs[i].lhs, cr)
        end
    end

    return ScaledMeanfieldEquations(eqs_out,de.operator_equations,de.states,de.operators,
                                    de.hamiltonian,de.jumps,de.jumps_dagger,de.rates,de.iv,
                                    de.varmap,order,
                                    de.scale_aons,de.names,de.was_scaled
                                    )
end

function _cumulant_expansion(args::Vector,order::Int)
    # Get all possible partitions; partitions(args,1) corresponds to the moment of order length(args)
    parts = [partitions(args,i) for i=2:length(args)]

    args_sum = Any[]
    for p in parts
        for pj in p
            n = length(pj)
            args_prod = Any[-factorial(n-1)*(-1)^(n-1)]
            for p_=pj # Product over partition blocks
                if length(p_) > order # If the encountered moment is larger than order, apply expansion
                    push!(args_prod, _cumulant_expansion(p_, order))
                else # Else, average and add its product
                    op_ = QMul(1, p_)
                    push!(args_prod, _average(op_))
                end
            end
            # Add terms in sum
            push!(args_sum, *(args_prod...))
        end
    end
    return average(+(args_sum...))
end

"""
    cumulant(x,n=get_order(x);simplify=true,kwargs...)

Compute the `n`th cumulant of `x` (either an operator or an average).
The output is simplified when `simplify=true`. Further keyword arguments are
passed on to simplification.

Examples
========
```
julia> cumulant(a)
⟨a⟩

julia> cumulant(a*b)
(⟨a*b⟩+(-1*⟨a⟩*⟨b⟩))

julia> cumulant(a*b,1)
⟨a*b⟩

julia> cumulant(a*b,3)
0
```
"""
function cumulant(op::QMul,n::Int=get_order(op);simplify=true,kwargs...)
    order = get_order(op)
    if order < n
        return zero(op)
    else
        c = _cumulant(op.args_nc, n)
        if simplify
            return SymbolicUtils.simplify(op.arg_c*c)
        else
            return op.arg_c*c
        end
    end
end
cumulant(avg::SymbolicUtils.Term{<:AvgSym},args...;kwargs...) = cumulant(SymbolicUtils.arguments(avg)[1],args...;kwargs...)
cumulant(op::QSym,n::Int=1;kwargs...) = isone(n) ? average(op) : zero(op)

function _cumulant(args::Vector,m::Int=length(args))
    parts = [partitions(args,i) for i=1:m]
    args_sum = Any[]
    for i=1:length(parts)
        p = collect(parts[i])
        for j=1:length(p) # Terms in the sum
            n = length(p[j])
            args_prod = Any[factorial(n-1)*(-1)^(n-1)]
            for p_=p[j] # Product over partition blocks
                push!(args_prod, _average(QMul(1, p_)))
            end
            # Add terms in sum
            push!(args_sum, *(args_prod...))
        end
    end
    return average(+(args_sum...))
end

"""
    get_order(arg)

Compute the order of a given argument. This is the order used to decide whether
something should be expanded using a [`cumulant_expansion`](@ref) method.

Examples
=======
```
julia> get_order(a)
1

julia> get_order(a*b)
2

julia> get_order(1)
0
```
"""
get_order(avg::SymbolicUtils.Term{<:AvgSym}) = get_order(SymbolicUtils.arguments(avg)[1])
function get_order(t::SymbolicUtils.Symbolic)
    if SymbolicUtils.istree(t)
        return maximum(map(get_order, SymbolicUtils.arguments(t)))
    else
        return 0
    end
end
get_order(::Number) = 0
get_order(q::QMul) = length(q.args_nc)
function get_order(q::QAdd)
    order = Int[get_order(arg) for arg in SymbolicUtils.arguments(t)]
    return maximum(order)
end
get_order(::QSym) = 1


text:
function test(project::String, suite_pattern::Vector{String}; coverage::Bool=false, startup_file::Bool=false, nprocs::Int=1)
    patterns = FilenameMatch.(suite_pattern)
    suites = find_test_suites(project) do suite
        isempty(suite_pattern) && return true
        any(p->ismatch(p, suite), patterns)
    end
    suites = divide_tests(suites, nprocs)
    results = pmap(run_test_suites, suites)

    print(err)
    print(report)
    return
end

function setup_test_process(project, nprocs::Int)
    test_dir(project)
    addprocs(nprocs;exeflags=["--project"])
end

function create_shared_env(project::String)
    pkgspec = PackageSpec(path=project)
    ctx = Context(env=EnvCache(project))
    TestEnv.isinstalled!(ctx, pkgspec) || throw(TestEnv.TestEnvError("$pkg not installed 👻"))
    Pkg.instantiate(ctx)
end

function divide_tests(suites::Dict{String, String}, nprocs::Int)
    suites_procs = [Dict{String, String}() for _ in 1:nprocs]

    n, r = divrem(length(suites), nprocs)
    nsuites_procs = [n for _ in 1:nprocs]
    if !iszero(r)
        for idx in 1:r
            nsuites_procs[idx] += 1
        end
    end

    proc_idx = 1; count = 0
    for (name, path) in suites
        suites_procs[proc_idx][name] = path
        count += 1

        if count == nsuites_procs[proc_idx]
            proc_idx += 1
            count = 0
        end
    end
    return suites_procs
end

function run_test_suites(suites::Dict{String, String})
    old_stdout, old_stderr = stdout, stderr
    # see https://github.com/JuliaLang/julia/issues/43586
    stdout_buf = IOContext(Pipe(), :color=>Base.get_have_color())
    stderr_buf = IOContext(Pipe(), :color=>Base.get_have_color())
    redirect_stdout(stdout_buf)
    redirect_stderr(stderr_buf)
    stdout_reader = @async read(stdout_buf, String)
    stderr_reader = @async read(stderr_buf, String)

    local err
    try
        for (suite, path) in suites
            eval_in_workspace(suite, path)
        end
    catch err
    finally
        redirect_stdout(old_stdout)
        redirect_stderr(old_stderr)
        close(stdout_buf)
        close(stderr_buf)
    end
    return fetch(stdout_reader), fetch(stderr_reader), err
end

function eval_in_workspace(name::String, path::String)
    m = Core.eval(Base.__toplevel__, :(module $(Symbol(path)) end))
    return Core.eval(m,
        Expr(:toplevel,
            :(using Test),
            :(eval(x) = $(Expr(:core, :eval))($m, x)),
            :(include(x) = $(Expr(:top, :include))($m, x)),
            :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $m, x)),
            :(Test.@testset $name begin
                include($path)
            end)
        )
    )
end

function find_test_suites(f, project::String)
    test_suites = suites = Dict{String, String}()
    project_test_dir = test_dir(project)
    for (root, dirs, files) in walkdir(project_test_dir)
        if root != project_test_dir
            prefix = relpath(root, project_test_dir) * "/"
        else
            prefix = ""
        end
        for file in files
            startswith(file, "test_") || continue
            endswith(file, ".jl") || continue

            suite_name = prefix * file[6:end-3]
            f(suite_name) || continue
            suites[suite_name] = joinpath(root, file)
        end
    end
    return test_suites
end

test_dir(project::String, paths::String...) = joinpath(project, "test", paths...)


text:
function computeFano(arr, interval)
  arrHist = fit(Histogram, arr, 0:interval:Int(round(maximum(arr))))
  ff = (mean(arrHist.weights.^2)-mean(arrHist.weights)^2)/mean(arrHist.weights)
  return ff
end


text:
using MetidaBase
using Test, Tables, TypedTables, CSV

@testset "MetidaBase.jl" begin
    io       = IOBuffer();
    # Metida table names - auto
    mt = MetidaBase.metida_table([1,2,3], ["a", "b", "c"])
    @test names(mt) == [:x1, :x2]
    # Metida table names - defined
    mt = MetidaBase.metida_table([1,2,3], ["a", "b", "c"], names = (:a, :b))
    # Push row
    pushfirst!(mt, [0, " "])
    @test mt[1, :a] == 0
    # Set element
    mt[1, :b] = "_"
    @test mt[1, :b] == "_"
    # push Named tuple as row
    ntr = NamedTuple{(:b, :a)}(["d", 10])
    pushfirst!(mt, ntr)
    @test mt[1, :b] == "d"
    # Test
    @test mt[:, :b] == ["d", "_", "a", "b", "c"]
    # Test show
    @test_nowarn Base.show(io, mt)

    # Tst size
    @test size(mt, 1) == 5
    @test size(mt, 2) == 2

    # Tables rows method
    rows = Tables.rows(mt)

    # MetidaTable > TypedTables
    df = Table(mt)

    # Enumerate
    for (i,j) in enumerate(mt)
        @test mt[i, :a] == j[1]
        @test mt[i, :b] == j[2]
    end

    # Appent one table to another
    l1 = length(mt)
    mt2 = MetidaBase.metida_table([1,2,3], ["a", "b", "c"], names = (:a, :b))
    append!(mt, mt2)
    @test l1 + length(mt2) == length(mt)
    mt2 = MetidaBase.metida_table(["e", "f", "g"], [1,2,3], names = (:b, :a))
    append!(mt, mt2)

    # CSV compat test
    @test_nowarn CSV.write(io, mt)

    # Structures
    ############################################################################
    struct ExampleIDStruct <: MetidaBase.AbstractSubject
        #time
        #obs
        id::Dict
    end

    struct ExampleResultStruct{T} <: MetidaBase.AbstractSubjectResult{T}
        data::T
        result::Dict
    end
    exiddsv = Vector{ExampleIDStruct}(undef, 3)
    for i in 1:3
        exiddsv[i] = ExampleIDStruct(Dict(:a => 1, :b => 1))
    end
    exidds = MetidaBase.DataSet(exiddsv)
    ############################################################################
    # DATASET Tests
    # Test length
    @test length(MetidaBase.getdata(exidds)) == length(exidds)
    # getid
    MetidaBase.getid(exidds[2])[:a] = 3
    MetidaBase.getid(exidds[3])[:a] = 2
    MetidaBase.getid(exidds[2])[:b] = 2
    MetidaBase.getid(exidds[3])[:b] = 3
    MetidaBase.getid(exidds, :, :a)

    exrsdsv = Vector{ExampleResultStruct}(undef, length(exidds))
    for i in 1:length(exidds)
        exrsdsv[i] = ExampleResultStruct(exidds[i], Dict(:r1 => 3, :r2 => 4))
    end
    exrsds = MetidaBase.DataSet(exrsdsv)

    @test exrsds[:, :r1][1] == 3
    @test exrsds[1, :r1] == 3

    @test MetidaBase.getid(exidds[3], :a) == 2
    sort!(exidds, :a)
    @test MetidaBase.getid(exidds[3], :a) == 3
    MetidaBase.getid(exrsds, :, :a)

    @test_nowarn sort!(exrsds, :a)

    @test first(exrsds) == exrsds[1]

    MetidaBase.uniqueidlist(exidds, [:a])
    MetidaBase.uniqueidlist(exidds, :a)

    MetidaBase.subset(exidds, Dict(:a => 1))
    MetidaBase.subset(exrsds, Dict(:a => 1))
    MetidaBase.subset(exrsds, 1:2)

    @test_nowarn map(identity, exidds)

    filtexrsds = filter(x -> x.id[:a] == 2, exidds)
    filter!(x -> x.id[:a] == 2, exidds)

    @test length(filtexrsds) == length(exidds)
    @test filtexrsds[1].id[:a] == exidds[1].id[:a] == 2

    mt  = MetidaBase.metida_table(exrsds)
    mt  = MetidaBase.metida_table(exrsds; results = :r1, ids = :a)
    smt = MetidaBase.Tables.schema(mt)
    # TypedTables export
    @test_nowarn Table(exrsds; results = :r1, ids = [:a, :b])
    #Iterators data
    v1 = [1,2,-6,missing,NaN]
    #Iterators tests
    itr1 = MetidaBase.skipnanormissing(v1)
    for i in itr1
        @test !MetidaBase.isnanormissing(i)
    end
    eachindex(itr1)
    eltype(itr1)
    keys(itr1)
    @test length(itr1) == 3

    itr2 = MetidaBase.skipnonpositive(v1)
    for i in itr2
        @test MetidaBase.ispositive(i)
    end
    eachindex(itr2)
    eltype(itr2)
    keys(itr2)
    @test length(itr2) == 2

    #Ststutils
    MetidaBase.sdfromcv(0.4) ≈ 0.38525317015992666
    MetidaBase.varfromcv(0.4) ≈ 0.1484200051182734
    MetidaBase.cvfromvar(0.4) ≈ 0.7013021443295824
    MetidaBase.cvfromsd(0.4) ≈ 0.41654636115540644

end


text:
export எழுத்துவகை_பிரி, மாத்திரை_அளவெடு, துப்புரவு_செய், தமிழுரையா
export உயிரெழுத்துக்கள், புள்ளி, ஆய்தம்

if VERSION > v"0.7.0" #FIXME: should use DEV build that made the change
    using Unicode: graphemes
end

const உயிரெழுத்துக்கள் = ["அ", "ஆ", "இ", "ஈ", "உ", "ஊ", "எ", "ஏ", "ஐ", "ஒ", "ஓ", "ஔ"]
const வல்லினங்கள் = ["க", "ச", "ட", "த", "ப", "ற"]
const மெல்லினங்கள் = ["ங",  "ஞ", "ண", "ந", "ம", "ன"]
const இடையினங்கள் = ["ய", "ர", "ல", "வ", "ழ", "ள"]
const அகரங்கள் = [வல்லினங்கள்;  மெல்லினங்கள்;  இடையினங்கள்]
const உயிர்ப்பகுதிகள் = ['ா', 'ி', 'ீ', 'ு', 'ூ', 'ெ', 'ே', 'ை', 'ொ', 'ோ', 'ௌ']
const புள்ளி = '்'
const ஆய்தம் = "ஃ"
#TODO: add grantham letters

எழுத்துவகை_பிரி(l::Char) = எழுத்துவகை_பிரி(string(l))
function எழுத்துவகை_பிரி(எழுத்து::AbstractString) #FIXME: needs better return type idea
    எழுத்து == ஆய்தம் && return :ஆய்தம்
    எழுத்து in உயிரெழுத்துக்கள் && return :உயிர்
    எழுத்து in அகரங்கள் && return :உயிர்மெய்

    ஓரெழுத்து_உறுதிசெய்(எழுத்து)
    தமிழுரையா(எழுத்து) || throw(ArgumentError("கொடுத்த குறி $எழுத்து தமிழ் எழுத்தல்ல (E:ARG_NT)"))

    if length(எழுத்து) == 1
        return :எழுத்தல்ல #தமிழ் எண்கள், நாட்காட்டிக் குறிப்புகள், முதலியவை (இப்போதைக்கு கிரந்தங்களும்)
    end

    எழுத்துப்பகுதி = எழுத்து[nextind(எழுத்து, 1)]
    எழுத்துப்பகுதி == புள்ளி && return :மெய்
    எழுத்துப்பகுதி in உயிர்ப்பகுதிகள் && return :உயிர்மெய்
end
const classify_letter = எழுத்துவகை_பிரி

மாத்திரை_அளவெடு(l::Char) = மாத்திரை_அளவெடு(string(l))
function மாத்திரை_அளவெடு(எழுத்து::AbstractString)
    எழுத்துவகை = எழுத்துவகை_பிரி(எழுத்து)

    #could give maathirai values to ௹, ௐ, etc?
    எழுத்துவகை == :எழுத்தல்ல && throw(ArgumentError("கொடுத்த குறி $எழுத்து மாத்திரை கொண்டதல்ல (E:ARG_NL)"))

    எழுத்துவகை in (:மெய், :ஆய்தம்) && return 0

    if எழுத்துவகை == :உயிர்
        if எழுத்து in உயிரெழுத்துக்கள்[[2, 4, 6, 8, 9, 11, 12]]
            return 2
        else
            return 1
        end
    end

    #அகர உயிர்மெய்கள் ஒரு codepoint நீளம் கொண்டவை
    length(எழுத்து) == 1 && return 1

    உயிர்ப்பகுதி = எழுத்து[nextind(எழுத்து, 1)]
    if உயிர்ப்பகுதி in உயிர்ப்பகுதிகள்[[1, 3, 5, 7, 8, 10, 11]]
        return 2
    else
        return 1
    end
end
const get_vowel_length = மாத்திரை_அளவெடு

ஓரெழுத்து_உறுதிசெய்(l::Char) = true
function ஓரெழுத்து_உறுதிசெய்(எழுத்து::AbstractString)
    if length(graphemes(எழுத்து)) > 1
        throw(ArgumentError("ஒரு எழுத்தை மட்டும் கொடுக்கவும் (கொடுத்தது: $எழுத்து) (E:ARG_GM)"))
    end
    true
end
const assert_single_grapheme = ஓரெழுத்து_உறுதிசெய்

function தமிழுரையா(எழுத்து::AbstractString)
    #TODO: add options to allow space, etc.
    ismatch(r"^\p{Tamil}+$", எழுத்து)
end
const is_tamil_text = தமிழுரையா

"""
துப்புரவு_செய்(உரை) -> தமிழுரை
தமிழ் எழுத்துக்கள், வெற்றிடங்கள் தவிர அனைத்தும் நீக்கு
"""
function துப்புரவு_செய்(உரை::AbstractString, எச்சரி=false)
    #TODO: add options to leave punctutations, digits, etc. in
    if VERSION > v"0.7.0" #FIXME: should use DEV build that made the change
        தமிழுரை = replace(உரை, r"[^\p{Tamil}\s]" => "")
    else
        தமிழுரை = replace(உரை, r"[^\p{Tamil}\s]", "")
    end
    நீக்கப்பட்டன = length(உரை) - length(தமிழுரை)
    if எச்சரி && (நீக்கப்பட்டன > 0)
        warn("தமிழல்லாத $(நீக்கப்பட்டன) எழுத்துக்கள் உரையிலிருந்து நீக்கப்பட்டன")
    end
    return தமிழுரை
end
const cleanup = துப்புரவு_செய்


text:
#
# indexing.jl --
#
# Array dimensions, axes, and indexing.
#

"""
    ArraySize

is the union of types eligible to define array size.  Calling
`[to_size](@ref)(dims)` on any argument `dims` such that `isa(dims,ArraySize)`
is true yields an array size in canonical form, that is an instance of
`Dims{N}` which is an alias for an `N`-tuple of `Int`.

"""
const ArraySize = Union{Integer,Tuple{Vararg{Integer}}}

"""
    to_size(dims)

converts `dims` to an instance of `Dims{N}` which is an alias for an `N`-tuple
of `Int`.  Argument `dims` can be a scalar integer or a tuple of integers.
Argument `dims` is returned if already of the correct type.  This method may
also be called as:

    to_size(dim1, dim2, ...)

to let `to_size` deal with a variable number of arguments.

The union [`ArraySize`](@ref) matches the types of acceptable argument(s) for
`to_size(arg)`: scalar integers and tuples of integers.

This method is intended for fast conversion, call `check_size(dims)` to verify
that all dimensions in `dims` are nonnegative.

"""
to_size(dims::Dims) = dims
to_size(dims::Integer...) = to_size(dims)
to_size(dims::Tuple{Vararg{Integer}}) = to_int(dims)

"""
    ArrayAxis

is the super-type of an array axis, it is an alias to `AbstractUnitRange{Int}`.

"""
const ArrayAxis = AbstractUnitRange{Int}

"""
    MaybeArrayAxis

is the union of types of arguments eligible to define an array axis.  This
alias is identical to [`IndexRange`](@ref) but has a different purpose.

Calling [`to_axis(x::MaybeArrayAxis)`](@ref) is guaranteed to yield an instance
of `ArrayAxis`.

"""
const MaybeArrayAxis = Union{Integer,AbstractUnitRange{<:Integer}}

"""
    IndexRange

is the union of types that can define a range of indices along a single
dimension of an array, that is integers and integer-valued unit ranges.  This
alias is identical to [`MaybeArrayAxis`](@ref) but has a different purpose.

"""
const IndexRange = Union{Integer,AbstractUnitRange{<:Integer}}

"""
    to_axis(ind)

yields `ind` converted into an array axis.  Argument `ind` can be a dimension
length or an integer valued unit range.  The type of the result inherits is a
sub-type of [`ArrayAxis`](@ref) which is an alias to `AbstractUnitRange{Int}`.

"""
to_axis(ind::ArrayAxis) = ind
to_axis(ind::AbstractUnitRange{<:Integer}) = to_type(ArrayAxis, ind)
to_axis(ind::Int) = Base.OneTo(ind)
to_axis(ind::Integer) = to_axis(to_type(Int, ind))

"""
    ArrayAxes{N}

is the type of array axes, that is an `N`-tuple of [`ArrayAxis`](@ref).

"""
const ArrayAxes{N} = NTuple{N,ArrayAxis}

"""
    MaybeArrayAxes{N}

is the type of arguments eligible to specify array axes, that is an `N`-tuple
of `MaybeArrayAxis`.

Calling `to_axes(x::MaybeArrayAxes{N})` is guaranteed to yield an instance
of `ArrayAxes{N}`.

"""
const MaybeArrayAxes{N} = NTuple{N,MaybeArrayAxis}

"""
    to_axes(inds)

yields `inds` converted in a proper list of array axes, that is a tuple of
[`ArrayAxis`](@ref) instances.  This method may also be called as:

    to_axes(ind1, ind2, ...)

to let `to_axes` deal with a variable number of arguments.

"""
to_axes(inds::Tuple{Vararg{ArrayAxis}}) = inds
to_axes(inds::MaybeArrayAxis...) = to_axes(inds)
to_axes(inds::Tuple{Vararg{MaybeArrayAxis}}) = map(to_axis, inds)

# Types for the sign of an offset.
const PlusMinus = Union{typeof(+),typeof(-)}

"""
    standard_size(A) -> size(A)

yields the list of dimensions of `A`, that is `size(A)`, throwing an
`ArgmentError` exception if `A` does not have standard 1-based indices.

See also [`has_standard_indexing`](@ref), [`same_standard_size`](@ref).

"""
function standard_size(A::AbstractArray)
    has_standard_indexing(A) || throw_non_standard_indexing()
    return size(A)
end

@noinline throw_non_standard_indexing() =
    throw(ArgumentError("arrays have non-standard indices"))

"""
    same_standard_size(A, B...) -> size(A)

checks whether arrays `A`, `B`, etc., all have standard indexing and the same
size which is returned.  If array sizes are not all identical, a
`DimensionMismatch` exception is thrown.  If arrays have non-standard indexing
(that is indices not starting at index one), an `ArgumentError` exception is
thrown.

See also [`standard_size`](@ref), [`has_standard_indexing`](@ref).

"""
same_standard_size(A::AbstractArray) = standard_size(A)
@inline function same_standard_size(A::AbstractArray, B::AbstractArray...)
    siz = standard_size(A)
    all_match(siz, standard_size, B...) || throw_not_same_size()
    return siz
end

"""
    same_size(A, B...) -> size(A)

checks whether arrays `A`, `B`, etc., all have the same size which is returned.
A `DimensionMismatch` exception is thrown if array sizes are not all identical.

See also [`same_standard_size`](@ref), [`same_axes`](@ref).

"""
same_size(A::AbstractArray) = size(A)
@inline function same_size(A::AbstractArray, B::AbstractArray...)
    siz = size(A)
    all_match(siz, size, B...) || throw_not_same_size()
    return siz
end

@noinline throw_not_same_size() =
    throw(DimensionMismatch("arrays must have same dimensions"))

"""
    same_axes(A, B...) -> axes(A)

checks whether arrays `A`, `B`, etc., have the same axes and returns them.  A
`DimensionMismatch` exception is thrown if axes are not all identical.

See also [`same_size`](@ref), [`all_indices`](@ref).

"""
same_axes(A::AbstractArray) = axes(A)
@inline function same_axes(A::AbstractArray, B::AbstractArray...)
    inds = axes(A)
    all_match(inds, axes, B...) || throw_not_same_axes()
    return inds
end

@noinline throw_not_same_axes() =
    throw(DimensionMismatch("arrays must have same indices"))

"""
     check_size(siz) -> len

checks the validity of the array size `siz` and yields the corresponding number
of elements (throwing an `ArgumentError` exception if this is not the case).
To be a valid array size, the values of `siz` must all be nonnegative.

See also [`to_size`](@ref).

"""
check_size(siz::Integer) =
    siz ≥ 0 ? to_int(siz) : throw(ArgumentError("invalid dimension"))

function check_size(siz::NTuple{N,Integer}) where {N}
    len = 1
    @inbounds for i in 1:N
        dim = to_int(siz[i])
        dim ≥ 0 || throw(ArgumentError("invalid dimension(s)"))
        len *= dim
    end
    return len
end

"""
    nonnegative(x) -> bool

yields whether `x` is greater or equal zero.

"""
nonnegative(x) = x ≥ zero(x)

"""
    has_standard_indexing(A)

return `true` if the indices of `A` start with 1 along all axes.  Can be called
with multiple arguments:

    has_standard_indexing(A, B, ...)

is equivalent to:

    has_standard_indexing(A) && has_standard_indexing(B) && ...

Opposite of `Base.has_offset_axes` which is not available in version of Julia
older than 0.7.

"""
has_standard_indexing(arg) = allof(x -> first(x) == 1, axes(arg)...)
has_standard_indexing(args...) = allof(has_standard_indexing, args...)
has_standard_indexing(::Array) = true
has_standard_indexing(::Colon) = true
has_standard_indexing(::Number) = true
has_standard_indexing(::Tuple) = true

"""
    cartesian_indices(A)
    cartesian_indices((n1, n2, ...))
    cartesian_indices((i1:j1, i2:j2, ...))
    cartesian_indices(CartesianIndex(i1, i2, ...), CartesianIndex(j1, j2, ...))
    cartesian_indices(R)

all yield an instance of `CartesianIndices` suitable for multi-dimensional
indexing of respectively: all the indices of array `A`, a multi-dimensional
array of dimensions `(n1,n2,...)`, a multi-dimensional region whose first and
last indices are `(i1,i2,...)` and `(j1,j2,...)` or a Cartesian region defined
by `R`, an instance of `CartesianIndices`.

"""
cartesian_indices(A::AbstractArray) = cartesian_indices(axes(A))
cartesian_indices(R::CartesianIndices) = R
@inline function cartesian_indices(start::CartesianIndex{N},
                                   stop::CartesianIndex{N}) where {N}
    CartesianIndices(map((i,j) -> i:j, start.I, stop.I))
end
cartesian_indices(inds::MaybeArrayAxis...) = cartesian_indices(inds)
cartesian_indices(inds::MaybeArrayAxes) = CartesianIndices(to_axes(inds))

# The following, would yield an `AbstractUnitRange` if a single argument is
# provided that is an integer or a range (not a tuple).  This lead to
# ambiguities so it is has been disabled.  The rules are that (i) an array
# argument (including a range) yields the indices for this array, (ii) a tuple
# of dimensions, a tuple of unit ranges or a pair of `CartesianIndices` yields
# the indices of the Cartesian region defined by these arguments.
#
#cartesian_indices(dim::Int) = Base.OneTo(dim)
#cartesian_indices(dim::Integer) = Base.OneTo(to_int(dim))
#cartesian_indices(rng::AbstractUnitRange{Int}) = rng
#cartesian_indices(rng::AbstractUnitRange{<:Integer}) = convert(UnitRange{Int}, rng)

"""
    axis_limits(I) = (i0,i1)

yields the limits `i0` and `i1` of index range `I` as a 2-tuple of `Int`'s and
such that `i0:i1` represents the same indices as `I` (although not in the same
order if `step(I) < 0`).  If `step(I)` is not equal to ±1, an `ArgumentError`
exception is thrown.

"""
axis_limits(I::AbstractUnitRange{<:Integer}) =
    (to_int(first(I)), to_int(last(I)))
axis_limits(I::AbstractRange{<:Integer}) =
    ((i0, i1, s) = (to_int(first(I)), to_int(last(I)), step(I));
     (s == +1 ? (i0,i1) :
      s == -1 ? (i1,i0) : throw_invalid_range_step()))

@noinline throw_invalid_range_step() =
    throw(ArgumentError("expecting a range with a step equal to ±1"))

"""

Assuming `A` and `B` are arrays with `N` dimensions:

    common_indices(A, B) -> inds

yields the set of all the indices that are valid for both `A` and `B`.  The
result is similar to `axes(A)` or `axes(B)`, that is an `N`-tuple of integer
valued unit ranges.

An offset `k` with a sign may be specified:

    common_indices(A, B, ±, k)

to obtain the set of all indices `i` such that `A[i]` and `B[i ± k]` are valid
and where here and above `±` is either `+` or `-`.  Offset `k` can be a tuple
of integers or a Cartesian index.

Arguments `A` and `B` may be both tuples of indices or index ranges or both
scalar or index range which specify the size or the axes of the arrays to be
indexed.  This is used in the following example, where we want to do `A[i] =
B[i]*C[i + k]` given the offset `k` and for all valid indices `i`:

    I = common_indices(same_axes(A, B), axes(C), +, k)
    @inbounds @simd for i in CartesianIndices(I)
       A[i] = B[i]*C[i + k]
    end

Note that `same_axes(A,B)` is called to get the axes of `A` and `B` while
asserting that they are the same, as a result no bound checking is necessary
and the loop can be optimized for vectorization.

""" common_indices

@inline function common_indices(A::AbstractArray{<:Any,N},
                                B::AbstractArray{<:Any,N}) where {N}
    common_indices(axes(A), axes(B))
end
@inline function common_indices(A::AbstractArray{<:Any,N},
                                B::AbstractArray{<:Any,N},
                                op::PlusMinus,
                                K::CartesianIndex{N}) where {N}
    common_indices(axes(A), axes(B), op, K)
end
@inline function common_indices(A::AbstractArray{<:Any,N},
                                B::AbstractArray{<:Any,N},
                                op::PlusMinus,
                                K::NTuple{N,Integer}) where {N}
    common_indices(axes(A), axes(B), op, K)
end

@inline function common_indices(I::NTuple{N,IndexRange},
                                J::NTuple{N,IndexRange}) where {N}
    map(common_indices, I, J)
end

@inline function common_indices(I::NTuple{N,IndexRange},
                                J::NTuple{N,IndexRange},
                                op::PlusMinus,
                                K::CartesianIndex{N}) where {N}
    common_indices(I, J, op, K.I)
end

@inline function common_indices(I::NTuple{N,IndexRange},
                                J::NTuple{N,IndexRange},
                                ::typeof(+),
                                K::NTuple{N,Integer}) where {N}
    map(_common_indices_plus, I, J, K)
end

@inline function common_indices(I::NTuple{N,IndexRange},
                                J::NTuple{N,IndexRange},
                                ::typeof(-),
                                K::NTuple{N,Integer}) where {N}
    map(_common_indices_minus, I, J, K)
end

@inline common_indices(i::IndexRange, j::IndexRange) =
    max(to_int(first(i)),
        to_int(first(j))):min(to_int(last(i)),
                              to_int(last(j)))

@inline common_indices(i::IndexRange, j::IndexRange, ::typeof(+), k::Integer) =
    _common_indices_plus(i, j, k)

@inline common_indices(i::IndexRange, j::IndexRange, ::typeof(-), k::Integer) =
    _common_indices_minus(i, j, k)

@inline _common_indices_plus(i::IndexRange, j::IndexRange, k::Integer) =
    max(to_int(first(i)),
        to_int(first(j) - k)):min(to_int(last(i)),
                                  to_int(last(j) - k))

@inline _common_indices_minus(i::IndexRange, j::IndexRange, k::Integer) =
    max(to_int(first(i)),
        to_int(first(j) + k)):min(to_int(last(i)),
                                  to_int(last(j) + k))

"""
    split_interval(I, J, +/-, k) -> Ia, Ib, Ic

given unit ranges `I` and `J` and offset `±k`, yields 3 unit ranges, such that
`Ia ∪ Ib ∪ Ic = I` and:

- `∀ i ∈ Ia`, `i ± k < first(J)`;

- `∀ i ∈ Ib`, `i ± k ∈ J`;

- `∀ i ∈ Ic`, `i ± k > last(J)`.

Unit ranges may be replaced by their first and last values:

    split_interval(first(I), last(I), first(J), last(J), +/-, k)

yields the same result as above.

""" split_interval

@inline function split_interval(I::AbstractUnitRange{<:Integer},
                                J::AbstractUnitRange{<:Integer},
                                pm::PlusMinus, k::Integer)
    split_interval(first(I), last(I), first(J), last(J), pm, k)
end

@inline function split_interval(imin::Integer, imax::Integer,
                                jmin::Integer, jmax::Integer,
                                pm::PlusMinus, k::Integer)
    split_interval(to_int(imin), to_int(imax),
                   to_int(jmin), to_int(jmax), pm, to_int(k))
end

@inline function split_interval(imin::Int, imax::Int,
                                jmin::Int, jmax::Int, ::typeof(-), k::Int)
    split_interval(imin, imax, jmin, jmax, +, -k)
end

@inline function split_interval(imin::Int, imax::Int,
                                jmin::Int, jmax::Int, ::typeof(+), k::Int)
    #
    # Intervals for different conditions: (a) below lower bound, (b) within
    # bounds or (c) beyond upper bound.
    #
    # (a)  i ∈ imin:imax  and  j = i + k < jmin
    #
    Ia = imin:min(imax,jmin-k-1)
    #
    # (b)  i ∈ imin:imax  and  jmin ≤ j = i + k ≤ jmax
    #
    Ib = max(imin,jmin-k):min(imax,jmax-k)
    #
    # (c)  i ∈ imin:imax  and  j = i + k > jmax
    #
    Ic = max(imin,jmax-k+1):imax
    return Ia, Ib, Ic
end

"""
    all_indices(A...)

yields an iterable object for visiting each index of array(s) `A` in an
efficient manner. For array types that have opted into fast linear indexing
(like `Array`), this is simply the range `1:length(A)`. For other array types,
return a specialized Cartesian range to efficiently index into the array(s) with
indices specified for every dimension.

If more than one `AbstractArray` argument are supplied, `all_indices` will
create an iterable object that is fast for all arguments (a `UnitRange` if all
inputs have fast linear indexing, a `CartesianIndices` otherwise).  A
`DimensionMismatch` exception is thrown if the arrays have different axes so
that it is always safe to use `@inbounds` in front of a loop like:

   for i in all_indices(A, B, C, D)
       A[i] = B[i]*C[i] + D[i]
   end

when `A`, `B` etc. are all (abstract) arrays.

This method is similar to `eachindex` except that a `DimensionMismatch`
exception is thrown if arrays have different axes.  For linearly indexed
arrays, `eachindex` only checks that they have the same linear index range
(that is the same number of elements, not the same shape).

""" all_indices

@inline all_indices(A::AbstractArray) =
    eachindex(IndexStyle(A), A)
@inline all_indices(A::AbstractArray, B::AbstractArray) =
    all_indices(IndexStyle(A, B), A, B)
@inline all_indices(A::AbstractArray, B::AbstractArray...) =
    all_indices(IndexStyle(A, B...), A, B...)

@inline function all_indices(::IndexLinear, A::AbstractArray,
                             B::AbstractArray...)
    all_match(axes(A), axes, B...) ||
        throw_indices_mismatch(IndexLinear(), A, B...)
    return eachindex(IndexLinear(), A)
end

@inline function all_indices(::IndexCartesian, A::AbstractArray,
                             B::AbstractArray...)
    inds = axes(A)
    all_match(inds, axes, B...) ||
        throw_indices_mismatch(IndexCartesian(), A, B...)
    # The following is the same as `eachindex(IndexCartesian(),A)` which yields
    # `CartesianIndices(axes(A))`.
    return CartesianIndices(inds)
end

@noinline function throw_indices_mismatch(::IndexLinear,
                                          A::AbstractArray...)
    io = IOBuffer()
    write(io, "all arguments of `all_indices` must have the same shape, got ")
    m = length(A)
    for i in 1:m
        write(io, (i == 1 ? "(" : i == m ? " and (" : ", )"))
        dims = size(A[i])
        n = length(dims)
        for j in 1:n
            j > 1 && write(io, ",")
            print(io, dims[j])
        end
        write(io, (n == 1 ? ",)" : ")"))
    end
    throw(DimensionMismatch(String(take!(io))))
end

@noinline function throw_indices_mismatch(::IndexCartesian,
                                          A::AbstractArray...)
    io = IOBuffer()
    write(io, "all arguments of `all_indices` must have the same axes, got ")
    m = length(A)
    for i in 1:m
        write(io, (i == 1 ? "(" : i == m ? " and (" : ", )"))
        inds = axes(A[i])
        n = length(inds)
        for j in 1:n
            j > 1 && write(io, ",")
            print(io, first(inds[j]), ":", last(inds[j]))
        end
        write(io, (n == 1 ? ",)" : ")"))
    end
    throw(DimensionMismatch(String(take!(io))))
end


text:
function lab_14(r::Robot)
    for i ∈ (Nord, Sud, West, Ost)
        go_to_border_come_back_and_return_distance!(r, i; markers = true, go_around_barriers = true)
    end
    putmarker!(r)
end


text:
using Pkg
Pkg.add("NeuralNetDiffEq")
Pkg.add("Flux")
Pkg.add("DifferentialEquations")
Pkg.add("LinearAlgebra")


text:
"""The module provides utilities to convert numbers into and from byte representations,
functions to check whether the system is little-endian or big-endian, and functions to
create bitmasks.
"""
module Utils
"""
    to_bytes(num::Number) -> Vector{UInt8}

Creates the byte representation of the number `num`.

# Arguments
- `num::Number`: the number from which we retrieve the bytes.

# Returns
- `Vector{UInt8}`: the bytes representation of the number `num`

# Examples
```jldoctest
using CANalyze.Utils
bytes = Utils.to_bytes(UInt16(0xAAFF))

# output
2-element Vector{UInt8}:
 0xff
 0xaa
```
"""
function to_bytes(num::Number)::Vector{UInt8}
    return reinterpret(UInt8, [num])
end

"""
    from_bytes(type::Type{T}, array::AbstractArray{UInt8}) where {T <: Number} -> T

Creates a value of type `T` constituted by the byte-array `array`. If the `array` length
is smaller than the size of `T`, `array` is filled with enough zeros.

# Arguments
- `type::Type{T}`: the type to which the byte-array is transformed
- `array::AbstractArray{UInt8}`: the byte array

# Returns
- `T`: the value constructed from the byte sequence

# Examples
```jldoctest
using CANalyze.Utils
bytes = Utils.from_bytes(UInt16, UInt8[0xFF, 0xAA])

# output
0xaaff
```
"""
function from_bytes(type::Type{T}, array::AbstractArray{UInt8})::T where {T <: Number}
    if length(array) < sizeof(T)
        for i=1:(sizeof(T) - length(array))
            push!(array, UInt8(0))
        end
    end

    values = reinterpret(type, array)
    return values[1]
end

"""
    is_little_endian() -> Bool

Returns whether the system has little-endian byte-order

# Returns
- `Bool`: The system has little-endian byte-order
"""
function is_little_endian()::Bool
    x::UInt16 = 0x0001
    lst = reinterpret(UInt8, [x])

    if lst[1] == 0x01
        return true
    else
        return false
    end
end

"""
    is_big_endian() -> Bool

Returns whether the system has big-endian byte-order

# Returns
- `Bool`: The system has big-endian byte-order
"""
function is_big_endian()::Bool
    return !is_little_endian()
end

"""
    mask(::Type{T}, length::UInt8, shift::UInt8) where {T <: Integer} -> T

Creates a mask of type `T` with `length` number of bits and right-shifted by `shift`
number of bits.

# Arguments
- `Type{T}`: the type of the mask
- `length::UInt8`: the number of bits
- `shift::UInt8`: the right-shift

# Returns
- `T`: the mask defined by `length` and `shift`
"""
function mask(::Type{T}, length::UInt8, shift::UInt8)::T where {T <: Integer}
    ret::T = mask(T, length)
    ret <<= shift
    return ret
end

"""
    mask(::Type{T}, length::Integer, shift::Integer) where {T <: Integer} -> T

Creates a mask of type `T` with `length` number of bits and right-shifted by `shift`
number of bits.

# Arguments
- `Type{T}`: the type of the mask
- `length::Integer`: the number of bits
- `shift::Integer`: the right-shift

# Returns
- `T`: the mask defined by `length` and `shift`

# Examples
```jldoctest
using CANalyze.Utils
m = Utils.mask(UInt64, 32, 16)

# output
0x0000ffffffff0000
```
"""
function mask(::Type{T}, length::Integer, shift::Integer)::T where {T <: Integer}
    l = convert(UInt8, length)
    s = convert(UInt8, shift)
    return mask(T, l, s)
end

"""
    mask(::Type{T}, length::UInt8) where {T <: Integer} -> T

Creates a mask of type `T` with `length` number of bits.

# Arguments
- `Type{T}`: the type of the mask
- `length::UInt8`: the number of bits

# Returns
- `T`: the mask defined by `length`
"""
function mask(::Type{T}, length::UInt8)::T where {T <: Integer}
    ret = zero(T)
    if length > 0
        for i in 1:(length-1)
            ret += 1
            ret <<= 1
        end
        ret += 1
    end
    return ret
end

"""
    mask(::Type{T}, length::Integer) where {T <: Integer} -> T

Creates a mask of type `T` with `length` number of bits.

# Arguments
- `Type{T}`: the type of the mask
- `length::Integer`: the number of bits

# Returns
- `T`: the mask defined by `length`

# Examples
```jldoctest
using CANalyze.Utils
m = Utils.mask(UInt64, 32)

# output
0x00000000ffffffff
```
"""
function mask(::Type{T}, length::Integer)::T where {T <: Integer}
    l = convert(UInt8, length)
    return mask(T, l)
end

"""
    mask(::Type{T}) where {T <: Integer} -> T

Creates a full mask of type `T` with `8sizeof(T)` bits.

# Arguments
- `Type{T}`: the type of the mask

# Returns
- `T`: the full mask

# Examples
```jldoctest
using CANalyze.Utils
m = Utils.mask(UInt64)

# output
0xffffffffffffffff
```
"""
function mask(::Type{T})::T where {T <: Integer}
    return full_mask(T)
end

"""
    full_mask(::Type{T}) where {T <: Integer} -> T

Creates a full mask of type `T` with `8sizeof(T)` bits.

# Arguments
- `Type{T}`: the type of the mask

# Returns
- `T`: the full mask

# Examples
```jldoctest
using CANalyze.Utils
m = Utils.full_mask(Int8)

# output
-1
```
"""
function full_mask(::Type{T})::T where {T <: Integer}
    ret::T = zero(T)
    for i in 0:(8sizeof(T) - 2)
        ret += 1
        ret <<= 1
    end
    ret += 1
    return ret
end

"""
    zero_mask(::Type{T}) where {T <: Integer} -> T

Creates a zero mask of type `T` where every bit is unset.

# Arguments
- `Type{T}`: the type of the mask

# Returns
- `T`: the zero mask

# Examples
```jldoctest
using CANalyze.Utils
m = Utils.zero_mask(UInt8)

# output
0x00
```
"""
function zero_mask(::Type{T})::T where {T <: Integer}
    return zero(T)
end

"""
    bit_mask(::Type{T}, bits::Set{UInt16}) where {T <: Integer} -> T

Creates a bit mask of type `T` where every bit inside `bits` is set.

# Arguments
- `Type{T}`: the type of the mask
- `bits::Set{UInt16}`: the set of bits we want to set

# Returns
- `T`: the mask

# Examples
```jldoctest
using CANalyze.Utils
m = Utils.bit_mask(UInt8, Set{UInt16}([0,1,2,3,4,5,6,7]))

# output
0xff
```
"""
function bit_mask(::Type{T}, bits::Set{UInt16})::T where {T <: Integer}
    result = zero(T)
    for bit in bits
        result |= mask(T, 1, bit)
    end
    return T(result)
end

function bit_mask(::Type{T}, bits::Integer...)::T where {T}
    bits = Set{UInt16}(bits)
    return bit_mask(T, bits)
end

function bit_mask(::Type{S}, bits::AbstractArray{T,N})::S where {S,T,N}
    bits = Set{UInt16}(bits)
    return bit_mask(S, bits)
end

export to_bytes, from_bytes
export is_little_endian, is_big_endian
export mask, zero_mask, full_mask, bit_mask
end


text:
"Create an equivalent smooth circuit from the given circuit."
function smooth(root::Δ)::Δ
    new_root = smooth(root[end])
    node2dag(new_root)
end

function smooth(root::ΔNode)::ΔNode
    lit_nodes = literal_nodes(root)
    f_con(n) = (n, BitSet())
    f_lit(n) = (n, BitSet(variable(n)))
    f_a(n, call) = begin
        parent_scope = BitSet()
        smooth_children = Vector{ΔNode}(undef, num_children(n))
        map!(smooth_children, children(n)) do child
            (smooth_child, scope) = call(child)
            union!(parent_scope, scope)
            smooth_child
        end
        return (conjoin_like(n, smooth_children), parent_scope)
    end
    f_o(n, call) = begin
        parent_scope = mapreduce(c -> call(c)[2], union, children(n))
        smooth_children = Vector{ΔNode}(undef, num_children(n))
        map!(smooth_children, children(n)) do child
            (smooth_child, scope) = call(child)
            smooth_node(smooth_child, setdiff(parent_scope, scope), lit_nodes)
        end
        return (disjoin_like(n, smooth_children), parent_scope)
    end
    foldup(root, f_con, f_lit, f_a, f_o, Tuple{ΔNode, BitSet})[1]
end

"Return a smooth version of the node where the `missing_scope` variables are added to the scope, using literals from `lit_nodes`"
function smooth_node(node::ΔNode, missing_scope, lit_nodes)
    if isempty(missing_scope)
        return node
    else
        get_lit(l) = get!(() -> literal_like(node, l), lit_nodes, l)
        ors = map(collect(missing_scope)) do v
            lit = var2lit(Var(v))
            lit_node = get_lit(lit)
            not_lit_node = get_lit(-lit)
            disjoin_like(node, lit_node, not_lit_node)
        end
        return conjoin_like(node, node, ors...)
    end
end

"""
Forget variables from the circuit.
Warning: this may or may not destroy the determinism property.
"""
function forget(is_forgotten::Function, circuit::Δ)
    forgotten = Dict{ΔNode,ΔNode}()
    (_, true_node) = constant_nodes(circuit) # reuse constants when possible
    if isnothing(true_node)
        true_node = true_like(circuit[end])
    end
    forget_node(n::ΔNode) = forget_node(GateType(n),n)
    forget_node(::ConstantGate, n::ΔNode) = n
    forget_node(::LiteralGate, n::ΔNode) =
        is_forgotten(variable(n)) ? true_node : n
    function forget_node(::⋀Gate, n::ΔNode)
        forgotten_children = map(c -> forgotten[c], children(n))
        conjoin_like(n, forgotten_children...)
    end
    function forget_node(::⋁Gate, n::ΔNode)
        forgotten_children = map(c -> forgotten[c], children(n))
        disjoin_like(n, forgotten_children...)
    end
    for node in circuit
        forgotten[node] = forget_node(node)
    end
    node2dag(forgotten[circuit[end]])
end

function forget2(is_forgotten::Function, circuit::Δ)::Δ
    new_root = forget2(is_forgotten, circuit[end])
    node2dag(new_root)
end

function forget2(is_forgotten::Function, root::ΔNode)::ΔNode
    (_, true_node) = constant_nodes2(root) # reuse constants when possible
    if isnothing(true_node)
        true_node = true_like(root)
    end
    f_con(n) = n
    f_lit(n) = is_forgotten(variable(n)) ? true_node : n
    f_a(n, cn) = conjoin_like(n, cn)
    f_o(n, cn) = disjoin_like(n, cn)
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, ΔNode)
end

"Remove all constant leafs from the circuit"
function propagate_constants(circuit::Δ)
    proped = Dict{ΔNode,ΔNode}()
    propagate(n::ΔNode) = propagate(GateType(n),n)
    propagate(::LeafGate, n::ΔNode) = n
    function propagate(::⋀Gate, n::ΔNode)
        proped_children = map(c -> proped[c], children(n))
        if any(c -> is_false(c), proped_children)
            false_like(n)
        else
            proped_children = filter(c -> !is_true(c), proped_children)
            conjoin_like(n, proped_children...)
        end
    end
    function propagate(::⋁Gate, n::ΔNode)
        proped_children = map(c -> proped[c], children(n))
        if any(c -> is_true(c), proped_children)
            true_like(n)
        else
            proped_children = filter(c -> !is_false(c), proped_children)
            disjoin_like(n, proped_children...)
        end
    end
    for node in circuit
        proped[node] = propagate(node)
    end
    node2dag(proped[circuit[end]])
end

@inline normalize(n::ΔNode, old_n::ΔNode, kept) = ()

"Return the circuit conditioned on given constrains"
function condition(circuit::Δ, lit::Lit)::Δ
    new_root = condition(circuit[end], lit)
    node2dag(new_root)
end

function condition(root::ΔNode, lit::Lit)::ΔNode
    literals = literal_nodes(root)
    (false_node, ) = constant_nodes2(root) # reuse constants when possible
    if isnothing(false_node)
        false_node = false_like(root)
    end

    if !haskey(literals, -lit)
        root
    elseif haskey(literals, -lit) && !haskey(literals, lit)
        false_node
    else
        f_con(n) = (n, true)
        f_lit(n) = begin
            if variable(n) == lit2var(lit) && literal(n) != lit
                (false_node, false)
            else
                (n, true)
            end
        end
        f_a(n, cv) = begin
            kept = last.(cv)
            new_children = convert(Vector{ΔNode}, first.(cv)[kept])
            if all(kept)
                (conjoin_like(n, new_children), true)
            else
                (nothing, false)
            end
        end
        f_o(n, cv) = begin
            kept = last.(cv)
            new_children = convert(Vector{ΔNode}, first.(cv)[kept])
            if any(kept) && (length(new_children) == 1) && isliteralgate(new_children[1])
                (new_children[1], true)
            elseif any(kept)
                new_n = disjoin_like(n, new_children)
                normalize(new_n, n, kept)
                (new_n, true)
            else
                (nothing, false)
            end
        end
        foldup_aggregate(root, f_con, f_lit, f_a, f_o, Tuple{Union{Nothing,ΔNode}, Bool})[1]
    end
end

import Base.split

"Return the circuit after spliting on edge `edge` and variable `var`"
function split(circuit::Union{Δ, ΔNode}, edge::Tuple{ΔNode, ΔNode}, var::Var; depth=0, sanity_check=true)
    or = first(edge)
    and = last(edge)

    # sanity check
    if sanity_check
        @assert depth >= 0
        @assert GateType(or) isa ⋁Gate && GateType(and) isa ⋀Gate && and in children(or)
        if circuit isa Δ
            @assert or in circuit && and in circuit
        else
            # TODO
        end
        literals = literal_nodes(and)
        @assert haskey(literals, var2lit(var)) && haskey(literals, - var2lit(var))
    end

    # split
    new_children1 = map(children(and)) do c
        condition(c, var2lit(var))
    end

    new_children2 = map(children(and)) do c
        condition(c, - var2lit(var))
    end

    new_and1 = conjoin_like(and, new_children1)
    new_and2 = conjoin_like(and, new_children2)
    new_or = disjoin_like(or, [[new_and1, new_and2]; filter(c -> c != and, children(or))])
    new_or = copy(new_or, depth)

    replace_node(circuit, or, new_or), new_or
end

"Clone the `and` node and redirect one of its parents to the new copy"
function clone(circuit::Union{Δ, ΔNode}, or1::ΔNode, or2::ΔNode, and::ΔNode; depth=1)
    # sanity check
    @assert depth >= 1
    @assert GateType(or1) isa ⋁Gate && GateType(or2) isa ⋁Gate && GateType(and) isa ⋀Gate
    @assert and in children(or1) && and in children(or2)
    if circuit isa Δ
        @assert or1 in circuit && or2 in circuit && and in circuit
    else
        # TODO
    end

    # clone
    new_and = copy(and, depth)
    new_or2 = disjoin_like(or2, [[new_and]; filter(c -> c != and, children(or2))])

    replace_node(circuit, or2, new_or2)
end

function replace_node(circuit::Δ, old::ΔNode, new::ΔNode)::Δ
    node2dag(replace_node(circuit[end], old, new))
end

function replace_node(root::ΔNode, old::ΔNode, new::ΔNode)::ΔNode
    @assert GateType(old) == GateType(new)
    f_con(n) = old == n ? new : n
    f_lit = f_con
    f_a(n, cns) = old == n ? new : conjoin_like(n, cns)
    f_o(n, cns) = old == n ? new : begin
            new_n = disjoin_like(n, cns)
            normalize(new_n, n, trues(length(cns)))
            new_n
        end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, ΔNode)
end


text:
@enum htsFormatCategory unknown_category sequence_data variant_data index_file region_list category_maximum=32767
@enum htsExactFormat unknown_format binary_format text_format sam bam bai cram crai vcf bcf csi gzi tbi bed format_maximum=32767
@enum htsCompression no_compression gzip bgzf custom compression_maximum=32767

immutable Record
    tid::Int32 # chromosome ID defined by bam_hdr_t
    pos::Int32 # 0-based leftmost coordinate
    bin::UInt16 # bin calculated by bam_reg2bin
    qual::UInt8 # mapping quality
    l_qname::UInt8 #length of the query name
    flag::UInt16 # bitwise flag
    n_cigar::UInt16 # number of CIGAR operations
    l_qseq::Int32 # length of the query sequence (read)
    mtid::Int32 # chrosome ID of next read in template
    mpos::Int32 # 0-based of next read in template, defined by bam_hdr_t
    isize::Int32 # ??? template length ???

    l_data::Int32 #current length of data
    m_data::Int32 # maximum length of data
    data::Ptr{Cuchar} # qname-cigar-seq-qual-aux how to load them all?
    id::UInt64   # BAM_ID ???
end

type KString
    l::Csize_t #8
    m::Csize_t #8
    s::Ptr{Cchar} #8
    function KString(l,m,s)
        new(l,m,s)
    end
    function KString(str::ASCIIString)
        s = convert(Ptr{Cchar}, pointer(str.data))
        l = length(str)
        m = l
        KString(l,m,s)
    end
end

bitstype 32 BINS

type HTSFile
    bins::UInt32
    lineno::Int64

    l::Csize_t
    m::Csize_t
    s::Ptr{Cchar}

    fn::Ptr{Cchar}
    fn_aux::Ptr{Cchar}
    bgzf::Ptr{Void}

    category::htsFormatCategory
    format::htsExactFormat

    major::Cshort
    minor::Cshort

    compression_level::Cshort
    specific::Ptr{Void}
end

#typealias Header Void
#=
 @abstract Structure for the alignment header.
 @field n_targets   number of reference sequences
 @field l_text      length of the plain text in the header
 @field target_len  lengths of the reference sequences
 @field target_name names of the reference sequences
 @field text        plain text
 @field sdict       header dictionary
=#
type Header
    n_targets::Int32
    ignore_sam_err::Int32
    l_text::UInt32
    target_len::Ptr{UInt32}
    cigar_tab::Ptr{Int8}
    target_name::Ptr{Ptr{Cchar}}
    text::Ptr{Cchar}
    sdict::Ptr{Void}
end

function strptr(phdr::Ptr{Header})
    hdr = unsafe_load(phdr)
    text = convert(Ptr{UInt8}, hdr.text)
    l_text = hdr.l_text
    strptr(text,1,hdr.l_text)
end

function strptr(pkstr::Ptr{KString})
    kstr = unsafe_load(pkstr)
    p = convert(Ptr{UInt8},kstr.s)
    strptr(p,1,kstr.l)
end
function show(io::IO,ks::KString)
    data = Array{UInt8,1}(ks.l)
    for i = 1:ks.l
        data[i] = unsafe_load(ks.s,i)
    end
    print(io,ASCIIString(data))
end

function bam_hdr_init()
    ccall((:bam_hdr_init,libhts),Ptr{Header},())
end
@doc """ input: BGZF *fp
         output: bam_hdr_t*
""" ->
function bam_hdr_read(bam_hdl::Ptr{Void})
    ccall((:bam_hdr_read,libhts),Ptr{Header},(Ptr{Void},),bam_hdl)
end
@doc """ input:  BGZF *fp
                 const bam_hdr_t *h
         output: int
""" ->
function bam_hdr_write(fp,h::Ptr{Header})
    ccall((:bam_hdr_write,libhts),Cint,(Ptr{Void},Ptr{Header}),fp,h)
end
@doc """ void bam_hdr_destroy(bam_hdr_t *h)
""" ->
function bam_hdr_destroy(h)
    ccall((:bam_hdr_destroy,libhts),Void,(Ptr{Header},),h)
end
@doc """ int bam_name2id(bam_hdr_t *h,const char*ref)
""" ->
function bam_name2id(h::Ptr{Void},ref::Ptr{Cchar})
    ccall((:bam_name2id,libhts),Cint,(Ptr{Void},Ptr{Cchar}),h,ref)
end
@doc """ bam_hdr_t* bam_hdr_dup(const bam_hdr_t *h0)
""" ->
function bam_hdr_dup(h0)
    ccall((:bam_hdr_dup,libhts),Ptr{Header},(Ptr{Header},),h)
end
@doc """ bam1_t* bam_init1(void)
""" ->
function bam_init1()
    ccall((:bam_init1,libhts),Ptr{Record},())
end
@doc """ void bam_destroy1(bam1_t *b)
""" ->
function bam_destroy1(b::Ptr{Record})
    ccall((:bam_destroy1,libhts),Void,(Ptr{Record},),b)
end
@doc """ int bam_read1(BGZF *fp, bam1_t *b)
""" ->
function bam_read!(fp,b)
    ccall((:bam_read1,libhts),Cint,(Ptr{Void},Ptr{Record}),fp,b)
end
@doc """ int bam_write1(BGZF *fp, const bam1_t *b)
""" ->
function bam_write1(fp,b)
    ccall((:bam_write1,libhts),Cint,(Ptr{Void},Ptr{Record}),fp,b)
end
@doc """ bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc)
""" ->
function bam_copy(bdst,bsrc)
    ccall((:bam_copy1,libhts),Ptr{Record},(Ptr{Record},Ptr{Record}),bdst,bsrc)
end
@doc """ bam1_t *bam_dup1(const bam1_t *bsrc)
""" ->
function bam_dup(bsrc)
    ccall((:bam_dup1,libhts),Ptr{Record},(Ptr{Record},),bsrc)
end
@doc """ int bam_cigar2qlen(int n_cigar, const uint32_t *cigar)
""" ->
function bam_cigar2qlen(n_cigar::Cint,cigar::Ptr{UInt32})
    ccall((:bam_cigar2qlen,libhts),Ptr{Record},(Cint,Ptr{UInt32}),n_cigar,cigar)
end
@doc """ int bam_cigar2rlen(int n_cigar, const uint32_t *cigar)
""" ->
function bam_cigar2rlen(n_cigar::Cint,cigar::Ptr{UInt32})
    ccall((:bam_cigar2rlen,libhts),Ptr{Record},(Cint,Ptr{UInt32}),n_cigar,cigar)
end
@doc """ int32_t bam_endpos(const bam1_t *b)
""" ->
function bam_endpos(b)
    ccall((:bam_endpos,libhts),Int32,(Ptr{Record},),b)
end
@doc """ int bam_str2flag(const char *str)
""" ->
function bam_str2flag(str::Ptr{Cchar})
    ccall((:bam_str2flag,libhts),Cint,(Ptr{Cchar},),str)
end
@doc """ char *bam_flag2str(int flag)
""" ->
function bam_flag2str(flag::Cint)
    ccall((:bam_flag2str,libhts),Ptr{Cchar},(Cint,),flag)
end

##############################################
###      BAM/CRAM indexing                 ###
##############################################
#=
    // These BAM iterator functions work only on BAM files.  To work with either
    // BAM or CRAM files use the sam_index_load() & sam_itr_*() functions.
    #define bam_itr_destroy(iter) hts_itr_destroy(iter)
    #define bam_itr_queryi(idx, tid, beg, end) sam_itr_queryi(idx, tid, beg, end)
    #define bam_itr_querys(idx, hdr, region) sam_itr_querys(idx, hdr, region)
    #define bam_itr_next(htsfp, itr, r) hts_itr_next((htsfp)->fp.bgzf, (itr), (r), 0)
=#
@doc """ #define bam_itr_destroy(iter) hts_itr_destroy(iter)
         void hts_itr_destroy(hts_itr_t *iter)
""" ->
function bam_itr_destroy(iter)
    ccall((:hts_itr_destroy,libhts),Void,(Ptr{Void},),iter)
end

@doc """ #define bam_itr_queryi(idx, tid, beg, end) sam_itr_queryi(idx, tid, beg, end)
         hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end)
""" ->
function bam_itr_queryi(idx::Ptr{Void}, tid::Cint, beg::Cint, ed::Cint)
    ccall((:sam_itr_queryi,libhts), Ptr{Void},(Ptr{Void},Ptr{Cint},Ptr{Cint},Ptr{Cint}),(idx, tid, beg, ed))
end
@doc """ #define bam_itr_querys(idx, hdr, region) sam_itr_querys(idx, hdr, region)
         hts_itr_t *sam_itr_querys(const hts_idx_t *idx, bam_hdr_t *hdr, const char *region)
""" ->
function bam_itr_querys(idx, hdr, region::Ptr{Cchar})
    ccall((:sam_itr_querys,libhts),Ptr{Void},(Ptr{Void},Ptr{Void},Ptr{Cchar}),idx, hdr, region)
end

@doc """ #define bam_itr_next(htsfp, itr, r) hts_itr_next((htsfp)->fp.bgzf, (itr), (r), 0)
         int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data)
""" ->
function bam_itr_next(bgzf, itr, r)
    ccall((:hts_itr_next,libhts),Cint,(Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),bgzf, itr, r, C_NULL)
end

@doc """
// build .csi or .bai BAM index file.  Does not work with CRAM.
// It is recommended to use the sam_index_* functions below instead.
#define bam_index_load(fn) hts_idx_load((fn), HTS_FMT_BAI)
hts_idx_t *hts_idx_load(const char *fn, int fmt);
""" ->
function bam_index_load(fn)
    ccall((:hts_idx_load,libhts),Ptr{Void},(Ptr{Cchar},Cint),fn,1)
end

@doc """
// build .csi or .bai BAM index file.  Does not work with CRAM.
// It is recommended to use the sam_index_* functions below instead.
#define bam_index_build(fn, min_shift) (sam_index_build((fn), (min_shift)))
int sam_index_build(const char *fn, int min_shift)
""" ->
function bam_index_build(fn, min_shift)
    ccall((:sam_index_build,libhts),Cint,(Ptr{Cchar},Cint),fn,min_shift)
end

@doc """
/// Load a BAM (.csi or .bai) or CRAM (.crai) index file
/** @param fp  File handle of the data file whose index is being opened
    @param fn  BAM/CRAM/etc filename to search alongside for the index file
    @return  The index, or NULL if an error occurred.
*/
hts_idx_t *sam_index_load(htsFile *fp, const char *fn);
""" ->
function sam_index_load(fp, fn::Ptr{Cchar})
    ccall((:sam_index_load,libhts),Ptr{Void},(Ptr{Void},Ptr{Cchar}),fp,fn)
end

@doc """
/// Load a specific BAM (.csi or .bai) or CRAM (.crai) index file
/** @param fp     File handle of the data file whose index is being opened
    @param fn     BAM/CRAM/etc data file filename
    @param fnidx  Index filename, or NULL to search alongside @a fn
    @return  The index, or NULL if an error occurred.
*/
hts_idx_t *sam_index_load2(htsFile *fp, const char *fn, const char *fnidx);
""" ->
function sam_index_load2(fp, fn::Ptr{Cchar}, fnidx::Ptr{Cchar});
    ccall((:sam_index_load2,libhts), Ptr{Void},(Ptr{Void},Ptr{Cchar},Ptr{Cchar}),fp,fn,fnidx)
end
@doc """
/// Generate and save an index file
/** @param fn        Input BAM/etc filename, to which .csi/etc will be added
    @param min_shift Positive to generate CSI, or 0 to generate BAI
    @return  0 if successful, or negative if an error occurred (usually -1; or
             -2: opening fn failed; -3: format not indexable)
*/
int sam_index_build(const char *fn, int min_shift);
""" ->
function sam_index_build(fn::Ptr{Cchar}, min_shift::Cint);
    ccall((:sam_index_build,libhts),Cint,(Ptr{Cchar},Cint),fn,  min_shift)
end
@doc """
/// Generate and save an index to a specific file
/** @param fn        Input BAM/CRAM/etc filename
    @param fnidx     Output filename, or NULL to add .bai/.csi/etc to @a fn
    @param min_shift Positive to generate CSI, or 0 to generate BAI
    @return  0 if successful, or negative if an error occurred.
*/
int sam_index_build2(const char *fn, const char *fnidx, int min_shift);
""" ->
function sam_index_build2(fn::Ptr{Cchar}, fnidx::Ptr{Cchar}, min_shift::Cint)
    ccall((:sam_index_build2,libhts),Cint,(Ptr{Cchar},Ptr{Cchar},Cint),fn, fnidx, min_shift)
end
@doc """
    #define sam_itr_destroy(iter) hts_itr_destroy(iter)
    void hts_itr_destroy(hts_itr_t *iter)
""" ->
function sam_itr_destroy(iter)
    ccall((:hts_itr_destroy,libhts),Void,(Ptr{Void},),iter)
end

@doc """    hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end);
""" ->
function sam_itr_queryi(idx, tid, beg, ed);
    ccall((:sam_itr_queryi,libhts),Ptr{Void},(Ptr{Void},Cint,Cint,Cint), idx, tid, beg, ed)
end
@doc """
    hts_itr_t *sam_itr_querys(const hts_idx_t *idx, bam_hdr_t *hdr, const char *region);
""" ->
function sam_itr_querys(idx, hdr::Ptr{Header}, region::Ptr{Cchar})
    ccall((:sam_itr_querys,libhts),Ptr{Void},(Ptr{Void},Ptr{Header},Ptr{Cchar}),idx, hdr, region)
end
@doc """
    #define sam_itr_next(htsfp, itr, r) hts_itr_next((htsfp)->fp.bgzf, (itr), (r), (htsfp))
    int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data)
""" ->
function sam_itr_next(bgzf, itr, r,handle)
    ccall((:hts_itr_next,libhts),Cint,(Ptr{Void},Ptr{Void},Ptr{Record},Ptr{Void}),bgzf,itr,r,handle)
end

#=
    /***************
     *** SAM I/O ***
     ***************/
=#
#=
@doc """
    #define sam_open(fn, mode) (hts_open((fn), (mode)))
    htsFile *hts_open(const char *fn, const char *mode)
""" ->
sam_open(fn::AbstractString,mode::AbstractString) = hts_open(fn,mode)
=#
@doc """
    #define sam_open_format(fn, mode, fmt) (hts_open_format((fn), (mode), (fmt)))
    htsFile *hts_open_format(const char *fn, const char *mode, const htsFormat *fmt);
""" ->
function sam_open_format(fn::Ptr{Cchar}, mode::Ptr{Cchar}, fmt)
    ccall((:hts_open_format,libhts),Ptr{Void},(Ptr{Cchar},Ptr{Cchar},Ptr{Void}),fn, mode, fmt)
end
@doc """
    #define sam_close(fp) hts_close(fp)
    int hts_close(htsFile *fp)
""" ->
function sam_close(fp)
    ret = ccall((:hts_close,libhts),Cint,(Ptr{Void},),fp)
    @assert ret == 0 || info("non-zeros status in the sam_close")
end
@doc """
    int sam_open_mode(char *mode, const char *fn, const char *format);
""" ->
function sam_open_mode(mode::Ptr{Cchar}, fn::Ptr{Cchar}, format::Ptr{Cchar});
    ccall((:sam_open_mode,libhts),Cint,(Ptr{Cchar},Ptr{Cchar},Ptr{Cchar}), mode, fn, format)
end
@doc """
    // A version of sam_open_mode that can handle ,key=value options.
    // The format string is allocated and returned, to be freed by the caller.
    // Prefix should be "r" or "w",
    char *sam_open_mode_opts(const char *fn,
                             const char *mode,
                             const char *format);
""" ->
function sam_open_mode_opts(fn::Ptr{Cchar},mode::Ptr{Cchar},format::Ptr{Cchar})
    ccall((:sam_open_mode_opts,libhts),Ptr{Cchar},(Ptr{Cchar},Ptr{Cchar},Ptr{Cchar}),fn,mode,format)
end

@doc """
    bam_hdr_t *sam_hdr_parse(int l_text, const char *text);
    translate a string to a bam_hdr_t type
""" ->
function sam_hdr_parse(l_text::Cint, text::Cstring)
    phdr = ccall((:sam_hdr_parse,libhts),Ptr{Header},(Cint,Cstring),l_text,text)
end

@doc """
    bam_hdr_t *sam_hdr_read(samFile *fp);
""" ->
function sam_hdr_read(fp::Ptr{Void})
    ccall((:sam_hdr_read,libhts),Ptr{Header},(Ptr{Void},),fp)
end

@doc """
    int sam_hdr_write(samFile *fp, const bam_hdr_t *h);
""" ->
function sam_hdr_write(fp::Ptr{Void}, h::Ptr{Header})
    ccall((:sam_hdr_write,libhts),Cint,(Ptr{Void},Ptr{Header}),fp,h)
end
@doc """
    int sam_parse1(kstring_t *s, bam_hdr_t *h, bam1_t *b);
""" ->
function sam_parse1(s::Ptr{KString}, h::Ptr{Header}, b::Ptr{Record})
    ccall((:sam_parse1,libhts),Cint,(Ptr{KString},Ptr{Header},Ptr{Record}), s, h, b)
end
@doc """
        int sam_format1(const bam_hdr_t *h, const bam1_t *b, kstring_t *str);
""" ->
function sam_format1(h::Ptr{Header}, b::Ptr{Record}, str::Ptr{KString})
    ccall((:sam_format1,libhts),Cint,(Ptr{Header},Ptr{Record},Ptr{KString}), h, b, str)
end
@doc """
    int sam_read1(samFile *fp, bam_hdr_t *h, bam1_t *b);
""" ->
function sam_read1(fp, h::Ptr{Header}, b::Ptr{Record})
    ccall((:sam_read1,libhts),Cint,(Ptr{Void},Ptr{Header},Ptr{Record}), fp, h,b)
end
@doc """
    int sam_write1(samFile *fp, const bam_hdr_t *h, const bam1_t *b);
""" ->
function sam_write1(fp, h::Ptr{Header}, b::Ptr{Record})
    ccall((:sam_write1,libhts),Cint,(Ptr{Void},Ptr{Header},Ptr{Record}),fp, h, b)
end

#=
    /*************************************
     *** Manipulating auxiliary fields ***
     *************************************/
=#
@doc """
    uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]);
""" ->
function bam_aux_get(b,tag::Ptr{Cchar})
    ccall((:bam_aux_get,libhts),UInt8,(Ptr{Void},Ptr{Cchar}), b, tag)
end
@doc """
    int32_t bam_aux2i(const uint8_t *s);
""" ->
function bam_aux2i(s::Ptr{UInt8})
    ccall((:bam_aux2i,libhts),Int32,(Ptr{UInt8},),s)
end
@doc """
    double bam_aux2f(const uint8_t *s);
""" ->
function bam_aux2f(s::Ptr{UInt8})
    ccall((:bam_aux2f,libhts),Cdouble,(Ptr{UInt8},),s)
end
@doc """
    char bam_aux2A(const uint8_t *s);
""" ->
function bam_aux2A(s::Ptr{UInt8})
    ccall((:bam_aux2A,libhts),Cchar,(Ptr{UInt8},),s)
end
@doc """
    char *bam_aux2Z(const uint8_t *s);
""" ->
function bam_aux2Z(s::Ptr{UInt8})
    ccall((:bam_aux2Z,libhts),Ptr{Cchar},(Ptr{UInt8},),s)
end
@doc """
    void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data);
""" ->
function bam_aux_append(b::Ptr{Record}, tag::Ptr{Cchar}, tp::Cchar, len::Cint, data::Ptr{UInt8})
    warn("Not sure how to deal with char tag[2]")
    ccall((:bam_aux_append,libhts),Void,(Ptr{Record},Ptr{Char},Cchar,Cint,UInt8), b, tag, tp, len, data)
end

@doc """
    int bam_aux_del(bam1_t *b, uint8_t *s);
""" ->
function bam_aux_del(b::Ptr{Record}, s::Ptr{UInt8})
    ccall((:bam_aux_del,libhts),Cint,(Ptr{Record},Ptr{UInt8}),b, s)
end

function bam_get_cigar(rec::Record)
    n_cigar = rec.n_cigar
    p = convert(Ptr{UInt32}, (rec.data + rec.l_qname))
    s = ""
    for i in 1:n_cigar
        tmp = unsafe_load(p,i)
        l4 = tmp & 0b1111
        h28 = tmp >> 4
        s = string(s,h28,cigarstring[l4+1])
    end
    s
end
function bam_seqi(s,i)
    (unsafe_load(s,(i-1)>>1+1) >> ((~(i-1)&1)<<2)) & 0xf
end
function bam_get_seq(rec::Record)
    pseq = rec.data + rec.n_cigar<<2 + rec.l_qname
    s = ""
    for i in 1:rec.l_qseq
        s =string(s, seq_nt16_str[bam_seqi(pseq, i)+1])
    end
    s
end
function bam_get_qual(rec::Record)
    len_qual = (rec.l_qseq+1)>>1
    p = rec.data + rec.n_cigar<<2 + rec.l_qname + (rec.l_qseq+1)>>1
    data = Array{UInt8,1}(len_qual)
    for i = 1:len_qual
        data[i] = unsafe_load(p,i) + 33
    end
    ASCIIString(data)
end
function bam_get_l_aux(rec::Record)
    p = rec.l_data - rec.n_cigar<<2 - rec.l_qname - rec.l_qseq - (rec.l_qseq + 1)>>1
end
function bam_get_aux(rec::Record)
    p = rec.data + rec.n_cigar<<2 + rec.l_qname + (rec.l_qseq+1)>>1 + rec.l_qseq
    len_aux = bam_get_l_aux(rec)
    data = Array{UInt8,1}(len_aux)
    for i = 1:len_aux
        data[i] = unsafe_load(p,i)
    end
    ASCIIString(data)
end

function hts_open{T<:AbstractString}(fn::T,mode::T)
    fp = pointer(fn.data)
    mdp = pointer(mode.data)
    ret = ccall((:hts_open,libhts),Ptr{Void},(Ptr{Cchar},Ptr{Cchar}),fp,mdp)
    if ret == C_NULL
        error("hts_open return C_NULL")
    end
    ret
end

function sam_read!(bam_hdl::Ptr{Void},bam_hdr_hdl::Ptr{Header},b::Ptr{Record})
    ccall((:sam_read1,libhts),Cint,(Ptr{Void},Ptr{Header},Ptr{Record}),bam_hdl,bam_hdr_hdl,b)
end

function sam_format!(h::Ptr{Header},b::Ptr{Record},pstr::Ptr{KString})
    ccall((:sam_format1,libhts),Cint,(Ptr{Header},Ptr{Record},Ptr{KString}),h,b,pstr)
end
function sam_parse!(pks,h,b)
    ccall((:sam_parse1,libhts),Cint,(Ptr{KString},Ptr{Void},Ptr{Record}),pks,h,b)
end


text:
# Computational Ghost Imaging

export ghostimage, simulate_ghostimage

function ghostimage(cs)
    N   = length(cs.measurements)
    Σs  = sum(cs.measurements)
    ΣsI = sum(cs.patterns .* cs.measurements, dims=1)
    ΣI  = sum(cs.patterns, dims=1)

    renormalize((ΣsI / N) - (Σs * ΣI / (N*N)))
end

simulate_ghostimage(img, m) = ghostimage(measure(img, m))


text:
@doc raw"""
    gu2w(g, u) -> w, un

Returns `w, un` (as a `NamedTuple`), from `g, u`, where:

```math
\mathbf{w} = g \frac{\mathbf{u}}{\|\mathbf{u}\|}
```

and `un` are the norms ``\|\mathbf{u}\|``.
"""
function gu2w(g::AbstractArray, u::AbstractArray)
    un = sqrt.(sum!(similar(g), abs2.(u)))
    return (w = g .* u ./ un, un = un)
end

@doc raw"""
    w2gu(w, un) -> g, u

Returns `g, u` (as a `NamedTuple`), such that

```math
\mathbf{w} = g \frac{\mathbf{u}}{\|\mathbf{u}\|}
```

where the norms ``\|\mathbf{u}\|`` are given by `un`.
"""
function w2gu(w::AbstractArray, un::AbstractArray)
    g = sqrt.(sum!(similar(un), abs2.(w)))
    u = w ./ g .* un
    return (g = g, u = u)
end

@doc raw"""
    ∂wnorm(∂w, g, u)

Given the gradients `∂w` of a function `f(w)` with respect to `w`,
returns the gradients `∂g, ∂u` of `f` with respect to the re-parameterization:

```math
\mathbf{w} = g \frac{\mathbf{u}}{\|\mathbf{u}\|}
```
"""
function ∂wnorm(∂w::AbstractArray, g::AbstractArray, u::AbstractArray)
    @assert size(∂w) == size(u)
    # see Eqs. (3) and (4) of Salimans & Kingma 2016
    un = sqrt.(sum!(similar(g), abs2.(u)))
    ∂g = sum!(similar(g), ∂w .* u) ./ un
    ∂u = (∂w - ∂g .* u ./ un) .* g ./ un
    return (g = ∂g, u = ∂u)
end