text:
__precompile__()
module Fluids

import Base: dot ,sqrt , +, -, *, /, ^, mod, mean, std, size, IndexStyle, getindex, setindex!, eltype, transpose

import Derivatives: dx, dy, dz

using ReadGlobal, MappedArrays, StaticArrays

export Grid, Scalar, ConstScalar, FVector, Tensor, SymTensor, dx, dy, dz, ∇, ∇² , AntiSymTensor, ⊗

struct Grid
  xl::Float64
  yl::Float64
  zl::Float64
end

# =============================== Scalars ==========================================

abstract type AbstractScalar{T,N} <: AbstractArray{T,N} end

^(s::AbstractScalar,n::T) where T<:Union{AbstractScalar,Number} = Scalar(s.grid, s .^ n)
^(s::AbstractScalar,n::Integer) = Scalar(s.grid, s .^ n)
^(n::Number,s::AbstractScalar) = Scalar(s.grid, n .^ s)
+(s::AbstractScalar,n::T) where T<:Union{Number,AbstractScalar} = Scalar(s.grid, s .+ n)
+(s::AbstractScalar,n::T) where T<:Union{Number} = Scalar(s.grid, s .+ n)
+(n::Number,s::AbstractScalar) = Scalar(s.grid, s .+ n)
-(s::AbstractScalar,n::T) where T<:Union{Number,AbstractScalar} = Scalar(s.grid, s .- n)
-(s::AbstractScalar,n::T) where T<:Union{Number} = Scalar(s.grid, s .- n)
-(n::Number,s::AbstractScalar) = Scalar(s.grid, n .- s)
/(s::AbstractScalar,n::T) where T<:Union{Number,AbstractScalar} = Scalar(s.grid, s ./ n)
/(s::AbstractScalar,n::T) where T<:Union{Number} = Scalar(s.grid, s ./ n)
/(n::Number,s::AbstractScalar) = Scalar(s.grid, n ./ s)
*(s::AbstractScalar,n::T) where T<:Union{Number,AbstractScalar} = Scalar(s.grid, s .* n)
*(s::AbstractScalar,n::T) where T<:Union{Number} = Scalar(s.grid, s .* n)
*(n::Number,s::AbstractScalar) = Scalar(s.grid, n .* s)
sqrt(s::AbstractScalar) = Scalar(s.grid, sqrt.(s))
-(s::AbstractScalar) = Scalar(s.grid, -(s))
mean(s::AbstractScalar) = mean(s.data)
std(s::AbstractScalar) = std(s.data)

std(n::T) where T<:Number = zero(T)

struct Scalar{T<:Number,N} <: AbstractScalar{T,N}
  grid::Grid
  data::Array{T,N}
end

Scalar(filename::String) = begin
  nx,ny,nz,xl,yl,zl = getdimsize()
  Scalar(Grid(xl,yl,zl),readfield(filename,nx,ny,nz))
end

size(S::Scalar) = size(S.data)
IndexStyle(::Type{T}) where {T<:Scalar} = Base.IndexLinear()
getindex(S::Scalar,I...) = S.data[I...]
setindex!(S::Scalar,v,I) =  setindex!(S.c,v,I)
eltype(S::Scalar{T,N}) where {T,N} = T

dx(field::Scalar,n::Integer=1) = Scalar(field.grid,dx(field.data,field.grid.xl,n))
dy(field::Scalar,n::Integer=1) = Scalar(field.grid,dy(field.data,field.grid.yl,n))
dz(field::Scalar,n::Integer=1) = Scalar(field.grid,dz(field.data,field.grid.zl,n))

struct ConstScalar{T<:Number,N} <: AbstractScalar{T,N}
  grid::Grid
  size::NTuple{N,Int}
  data::T
end

size(S::ConstScalar) = S.size
IndexStyle(::Type{T}) where {T<:ConstScalar} = Base.IndexCartesian()
getindex(S::ConstScalar,I...) = S.data
eltype(S::ConstScalar{T,N}) where {T,N} = T

^(s::ConstScalar,n::Number) = ConstScalar(s.grid,s.size, s.data^n)
^(n::Number,s::ConstScalar) = ConstScalar(s.grid,s.size, n^s.data)
+(s::ConstScalar,n::Number) = ConstScalar(s.grid,s.size, s.data+n)
+(n::Number,s::ConstScalar) = ConstScalar(s.grid, s.data + n)
-(s::ConstScalar,n::Number) = ConstScalar(s.grid,s.size, s.data - n)
-(n::Number,s::ConstScalar) = ConstScalar(s.grid,s.size, n - s.data)
/(s::ConstScalar,n::Number) = ConstScalar(s.grid,s.size, s.data / n)
/(n::Number,s::ConstScalar) = ConstScalar(s.grid,s.size, n / s.data)
*(s::ConstScalar,n::Number) = ConstScalar(s.grid, s.data * n)
*(n::Number,s::ConstScalar) = ConstScalar(s.grid,s.size, n * s.data)
sqrt(s::ConstScalar) = ConstScalar(s.grid,s.size, sqrt(s.data))
-(s::ConstScalar) = ConstScalar(s.grid,s.size, -s.data)

# =============================== Vectors ==========================================

struct FVector{T<:AbstractFloat,N}
  grid::Grid
  x::Array{T,N}
  y::Array{T,N}
  z::Array{T,N}
end

FVector(file1::String,file2::String,file3::String) = begin
  nx,ny,nz,xl,yl,zl = getdimsize()
  FVector(Grid(xl,yl,zl),readfield(file1,nx,ny,nz),readfield(file2,nx,ny,nz),readfield(file3,nx,ny,nz))
end

getindex(u::FVector,I...) = SVector{3}(u.x[I...],u.y[I...],u.z[I...])

dx(field::FVector,n::Integer=1) = FVector(field.grid,
                          dx(field.x,field.grid.xl,n),
                          dx(field.y,field.grid.xl,n),
                          dx(field.z,field.grid.xl,n))
dy(field::FVector,n::Integer=1) = FVector(field.grid,
                          dy(field.x,field.grid.yl,n),
                          dy(field.y,field.grid.yl,n),
                          dy(field.z,field.grid.yl,n))
dz(field::FVector,n::Integer=1) = FVector(field.grid,
                          dz(field.x,field.grid.zl,n),
                          dz(field.y,field.grid.zl,n),
                          dz(field.z,field.grid.zl,n))

+(v::FVector,u::FVector) = FVector(v.grid,v.x .+ u.x, v.y .+ u.y, v.z .+ u.z)
-(v::FVector,u::FVector) = FVector(v.grid,v.x .- u.x, v.y .- u.y, v.z .- u.z)
-(v::FVector) = FVector(v.grid,-v.x, -v.y, -v.z)
⊗(v::FVector,u::FVector) = Tensor(v.grid,
v.x .* u.x, v.x .* u.y, v.x .* u.z,
v.y .* u.x, v.y .* u.y, v.y .* u.z,
v.z .* u.x, v.z .* u.y, v.z .* u.z)

mean(v::FVector) = SVector{3}(mean(v.x),mean(v.y),mean(v.z))
std(v::FVector;mean=nothing) = mean === nothing ? SVector{3}(std(v.x),std(v.y),std(v.z)) : SVector{3}(std(v.x,mean=mean[1]),std(v.y,mean=mean[2]),std(v.z,mean=mean[3]))

# =============================== Tensors ==========================================

abstract type AbstractTensor{T,N} end

getindex(S::AbstractTensor,I...) = @SMatrix [S.xx[I...] S.xy[I...] S.xz[I...];
                                          S.yx[I...] S.yy[I...] S.yz[I...];
                                          S.zx[I...] S.zy[I...] S.zz[I...]]

-(A::AbstractTensor) = Tensor(A.grid,
-A.xx , -A.xy , -A.xz ,
-A.yx , -A.yy , -A.yz ,
-A.zx , -A.zy , -A.zz )

+(A::AbstractTensor,B::AbstractTensor) = Tensor(A.grid,
A.xx .+ B.xx,A.xy .+ B.xy,A.xz .+ B.xz,
A.yx .+ B.yx,A.yy .+ B.yy,A.yz .+ B.yz,
A.zx .+ B.zx,A.zy .+ B.zy,A.zz .+ B.zz)

-(A::AbstractTensor,B::AbstractTensor) = Tensor(A.grid,
A.xx .- B.xx, A.xy .- B.xy, A.xz .- B.xz,
A.yx .- B.yx, A.yy .- B.yy, A.yz .- B.yz,
A.zx .- B.zx, A.zy .- B.zy, A.zz .- B.zz)

dot(A::AbstractTensor,B::AbstractTensor) = Tensor(A.grid,
A.xx .* B.xx .+ A.xy .* B.yx + A.xz .* B.zx,
A.xx .* B.xy .+ A.xy .* B.yy + A.xz .* B.zy,
A.xx .* B.xz .+ A.xy .* B.yz + A.xz .* B.zz,
A.yx .* B.xx .+ A.yy .* B.yx + A.yz .* B.zx,
A.yx .* B.xy .+ A.yy .* B.yy + A.yz .* B.zy,
A.yx .* B.xz .+ A.yy .* B.yz + A.yz .* B.zz,
A.zx .* B.xx .+ A.zy .* B.yx + A.zz .* B.zx,
A.zx .* B.xy .+ A.zy .* B.yy + A.zz .* B.zy,
A.zx .* B.xz .+ A.zy .* B.yz + A.zz .* B.zz)

function ^(A::AbstractTensor,n::Integer)
  if n <= 0
    error("Only positive integers are accepted")
  elseif n == 1
    return A + zero(eltype(A))
  elseif n == 2
    return A ⋅ A
  else
    return A ⋅ (A^(n-1))
  end
end

mean(t::AbstractTensor) = SMatrix{3,3}(mean(t.xx),mean(t.xy),mean(t.xz),mean(t.yx),mean(t.yy),mean(t.yz),mean(t.zx),mean(t.zy),mean(t.zz))
std(t::AbstractTensor;mean=nothing) = mean === nothing ? SMatrix{3,3}(std(t.xx),std(t.xy),std(t.xz),std(t.yx),std(t.yy),std(t.yz),std(t.zx),std(t.zy),std(t.zz)) : SMatrix{3,3}(std(t.xx,mean=mean[1,1]),std(t.xy,mean=mean[1,2]),std(t.xz,mean=mean[1,3]),std(t.yx,mean=mean[2,1]),std(t.yy,mean=mean[2,2]),std(t.yz,mean=mean[2,3]),std(t.zx,mean=mean[3,1]),std(t.zy,mean=mean[3,2]),std(t.zz,mean=mean[3,3]))

struct Tensor{T<:AbstractFloat,N} <: AbstractTensor{T,N}
  grid::Grid
  xx::Array{T,N}
  xy::Array{T,N}
  xz::Array{T,N}
  yx::Array{T,N}
  yy::Array{T,N}
  yz::Array{T,N}
  zx::Array{T,N}
  zy::Array{T,N}
  zz::Array{T,N}
end

struct SymTensor{T<:AbstractFloat,N} <: AbstractTensor{T,N}
  grid::Grid
  xx::Array{T,N}
  yy::Array{T,N}
  zz::Array{T,N}
  xy::Array{T,N}
  xz::Array{T,N}
  yz::Array{T,N}
  yx::Array{T,N}
  zx::Array{T,N}
  zy::Array{T,N}

  function SymTensor(grid::Grid, xx::Array{T,N}, yy::Array{T,N}, zz::Array{T,N}, xy::Array{T,N}, xz::Array{T,N}, yz::Array{T,N}) where {T<:AbstractFloat,N}
    return new{T,N}(grid,xx,yy,zz,xy,xz,yz,xy,xz,yz)
  end
end

function SymTensor(f1::String,f2::String,f3::String,f4::String,f5::String,f6::String)
  nx,ny,nz,xl,yl,zl = getdimsize()
  SymTensor(Grid(xl,yl,zl),readfield(f1,nx,ny,nz),readfield(f2,nx,ny,nz),readfield(f3,nx,ny,nz),readfield(f4,nx,ny,nz),readfield(f5,nx,ny,nz),readfield(f6,nx,ny,nz))
end

SymTensor(f1::String) = SymTensor(f1,replace(f1,"11","22"),replace(f1,"11","33"),replace(f1,"11","12"),replace(f1,"11","13"),replace(f1,"11","23"))

struct ZeroArray{T,N} <: AbstractArray{T,N}
  size::NTuple{N,Int}
end
getindex(S::ZeroArray{T,N},I...) where {T,N} = zero(T)
eltype(S::ZeroArray{T,N}) where {T,N} = T
size(S::ZeroArray) = S.size

struct AntiSymTensor{T<:AbstractFloat,N} <: AbstractTensor{T,N}
  grid::Grid
  xy::Array{T,N}
  xz::Array{T,N}
  yz::Array{T,N}
  xx::ZeroArray{T,N}
  yy::ZeroArray{T,N}
  zz::ZeroArray{T,N}
  yx::MappedArray{T,N,Array{T,N},typeof(-),typeof(-)}
  zx::MappedArray{T,N,Array{T,N},typeof(-),typeof(-)}
  zy::MappedArray{T,N,Array{T,N},typeof(-),typeof(-)}

  function AntiSymTensor{T,N}(grid::Grid,xy::Array{T,N}, xz::Array{T,N}, yz::Array{T,N}) where {T,N}
    return new{T,N}(grid,xy,xz,yz,ZeroArray{T,N}(size(xy)),ZeroArray{T,N}(size(xy)),ZeroArray{T,N}(size(xy)),mappedarray((-,-),xy),mappedarray((-,-),xz),mappedarray((-,-),yz))
  end
end

function AntiSymTensor(f1::String,f2::String,f3::String)
  nx,ny,nz,xl,yl,zl = getdimsize()
  u12 = readfield(f1,nx,ny,nz)
  AntiSymTensor{eltype(u12),3}(Grid(xl,yl,zl),u12,readfield(f2,nx,ny,nz),readfield(f3,nx,ny,nz))
end

function AntiSymTensor(f1::String)
  nx,ny,nz,xl,yl,zl = getdimsize()
  u12 = readfield(f1,nx,ny,nz)
  AntiSymTensor{eltype(u12),3}(Grid(xl,yl,zl),u12,readfield(replace(f1,"12","13"),nx,ny,nz),readfield(replace(f1,"12","23"),nx,ny,nz))
end

# =============================== Mixed Functions ==========================================

dot(t::AbstractTensor,u::FVector{T,N}) where {T,N} = FVector(u.grid,
t.xx .* u.x .+ t.xy .* u.y + t.xz .* u.z,
t.yx .* u.x .+ t.yy .* u.y + t.yz .* u.z,
t.zx .* u.x .+ t.zy .* u.y + t.zz .* u.z)

dot(u::FVector{T,N},t::AbstractTensor) where {T,N} = FVector(u.grid,
t.xx .* u.x .+ t.yx .* u.y + t.zx .* u.z,
t.xy .* u.x .+ t.yy .* u.y + t.zy .* u.z,
t.xz .* u.x .+ t.yz .* u.y + t.zz .* u.z)

function ∇(field::Scalar)
   FVector(field.grid,dx(field).data,dy(field).data,dz(field).data)
end
function ∇(field::FVector)
  Tensor(field.grid,
  dx(field.x,field.grid.xl),dx(field.y,field.grid.xl),dx(field.z,field.grid.xl),
  dy(field.x,field.grid.yl),dy(field.y,field.grid.yl),dy(field.z,field.grid.yl),
  dz(field.x,field.grid.zl),dz(field.y,field.grid.zl),dz(field.z,field.grid.zl),)
end

function ∇²(field::Scalar)
  Scalar(field.grid,
  dx(field.data,field.grid.xl,2) .+ dy(field.data,field.grid.yl,2) .+ dz(field.data,field.grid.zl,2))
end

function ∇²(field::FVector)
  FVector(field.grid,
  dx(field.x,field.grid.xl,2) .+ dy(field.x,field.grid.yl,2) .+ dz(field.x,field.grid.zl,2),
  dx(field.y,field.grid.xl,2) .+ dy(field.y,field.grid.yl,2) .+ dz(field.y,field.grid.zl,2),
  dx(field.z,field.grid.xl,2) .+ dy(field.z,field.grid.yl,2) .+ dz(field.z,field.grid.zl,2))
end

dot(v::FVector,u::FVector) = Scalar(v.grid, v.x .* u.x .+ v.y .* u.y .+ v.z .* u.z)

+(v::SymTensor,u::AntiSymTensor) = Tensor(v.grid,   v.xx     , v.xy .+ u.xy, v.xz .+ u.xz,
                                                 v.xy .- u.xy,     v.yy    , v.yz .+ u.yz,
                                                 v.xz .- u.xz, v.yz .- u.yz,     v.zz)
+(u::AntiSymTensor, v::SymTensor) = v+u

*(s::AbstractScalar,v::FVector) = FVector(v.grid,s .* v.x, s .* v.y, s .* v.z)

*(u::FVector, s::AbstractScalar) = *(s,u)

/(v::FVector,s::AbstractScalar) = FVector(v.grid, v.x ./ s, v.y ./ s, v.z ./ s)

dot(f::typeof(∇),u::FVector) = Scalar(u.grid, dx(u.x,u.grid.xl) .+ dy(u.y,u.grid.yl) .+ dz(u.z,u.grid.zl))

dot(f::typeof(∇),u::AbstractTensor) = FVector(u.grid,
dx(u.xx,u.grid.xl) .+ dy(u.yx,u.grid.yl) .+ dz(u.zx,u.grid.zl),
dx(u.xy,u.grid.xl) .+ dy(u.yy,u.grid.yl) .+ dz(u.zy,u.grid.zl),
dx(u.xz,u.grid.xl) .+ dy(u.yz,u.grid.yl) .+ dz(u.zz,u.grid.zl))

Base.cross(f::typeof(∇),u::FVector) = FVector(u.grid,
dy(u.z,u.grid.yl) .- dz(u.y,u.grid.zl),
dz(u.x,u.grid.zl) .- dx(u.z,u.grid.xl),
dx(u.y,u.grid.xl) .- dy(u.x,u.grid.yl))

mod(v::FVector) = Scalar(v.grid, sqrt.(v.x .^2 .+ v.y .^2 .+ v.z .^2))
mod(v::AbstractTensor) = Scalar(v.grid, @.(sqrt(v.xx^2 + v.xy^2 + v.xz^2 + v.yx^2 + v.yy^2 + v.yz^2 + v.zx^2 + v.zy^2 + v.zz^2)))

end # module


text:
# ---
# title: 128. Longest Consecutive Sequence
# id: problem128
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Union Find
# link: <https://leetcode.com/problems/longest-consecutive-sequence/description/>
# hidden: true
# ---
#
# Given an unsorted array of integers `nums`, return _the length of the longest
# consecutive elements sequence._
#
# **Follow up:**  Could you implement the `O(n)` solution?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [100,4,200,1,3,2]
#     Output: 4
#     Explanation: The longest consecutive elements sequence is [1, 2, 3, 4]. Therefore its length is 4.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [0,3,7,2,5,8,4,6,0,1]
#     Output: 9
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 104`
#   * `-109 <= nums[i] <= 109`
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
using StatsBase
using Base.Threads

import Base.sort
import Base.isless
import Base.==

type Node
    item_ids::Array{Int16,1}
    transacts::BitArray
end

function ==(a::Node, b::Node)
    res = a.item_ids == b.item_ids && a.transacts == b.transacts
    res
end

function ==(x::Array{Node,1}, y::Array{Node,1})
    n = length(x)
    res = length(y) == n
    i = 1
    while res && i < n
        println(i)
        res = x[i] == y[i]
        i += 1
    end
    res
end

"""
    isless(x::Array{Int,1}, y::Array{Int,1})
This `isless()` method is defined to handle two Array{Int,1} so that we can
use `sortperm()` function on arrays of arrays.
"""
function isless(x::Array{Int,1}, y::Array{Int,1})
    res = true
    if length(x) ≠ length(y)
        error("isless() cannot be used to compare vectors of differing lengths")
    else
        for i = 1:length(x)
            if x[i] > y[i]
                res = false
                break
            end
        end
    end
    res
end

function sort(x::Array{Array{Int,1},1})
    perm = sortperm(x)
    x[perm]
end

@code_warntype sortperm([[3, 13, 21], [3, 12, 14]])

function prefix(nd::Node, k::Int)
    res = view(nd.item_ids, 1:k-2)
    res
end

groceries = ["asparagus", "broccoli", "carrots", "cauliflower", "celery",
             "corn", "cucumbers", "lettuce", "mushrooms", "onions",
             "peppers", "potatos", "spinach", "zucchini", "tomatoes",
             "apples", "avocados", "bananas", "berries", "cherries",
             "grapefruit", "grapes", "kiwis", "lemons", "melon",
             "oranges", "peaches", "nectarines", "pears", "plums",
             "butter", "milk", "sour cream", "whipped cream", "yogurt",
             "bacon", "beef", "chicken", "ground beef", "turkey",
             "crab", "lobster", "oysters", "salmon", "shrimp",
             "tilapia", "tuna", "flour", "sugar", "yeast",
             "cookies", "crackers", "nuts", "oatmeal", "popcorn",
             "pretzels", "cosmetics", "floss", "mouthwash", "toothpaste",
             "lime", "almonds", "cashews", "ketchup", "mustard"]

transactions = [sample(groceries, 12, replace = false) for x in 1:1_000_000];

function get_unique_items{M}(T::Array{Array{M, 1}, 1})
    dict = Dict{M, Int}()

    for t in T
        for i in t
            dict[i] = 1
        end
    end
    return [x for x in keys(dict)]
end

function occurrence(T::Array{Array{String, 1}, 1})
    n = length(T)
    uniq_items = get_unique_items(T)
    sort!(uniq_items)
    p = length(uniq_items)
    res = BitArray(n, p)

    for j = 1:p
        for i = 1:n
            res[i, j] = uniq_items[j] ∈ T[i]
        end
    end
    res
end

t = [["a", "b"],
     ["b", "c", "d"],
     ["a", "c"],
     ["e", "b"],
     ["a", "c", "d"],
     ["a", "e"],
     ["a", "b", "c"],
     ["c", "b", "e"]]

@code_warntype get_unique_items(t)
@code_warntype occurrence(t)

occ1 = occurrence(t)

function is_frequent(indcs::Array{Int,1}, idx::Int, occ::BitArray{2}, minsupp::Int)
    # NOTE: This function's number of allocations is constant with
    # respect to the size of its input vectors (this is good!).
    n = size(occ, 1)
    bit_col = trues(n)
    for j in indcs
        bit_col &= view(occ, :, j)
    end
    cnt = sum(bit_col & view(occ, :, idx))

    res = cnt ≥ minsupp
    res
end

transactions = [sample(groceries, 4, replace = false) for x in 1:100_000];
occ2 = occurrence(transactions);
@code_warntype is_frequent([1], 2, occ1, 2)
@time is_frequent([1], 2, occ2, 2)

function bitwise_and!(x::BitArray, y::BitArray)
    # NOTE: This function's number of allocations is constant with respect to the size
    # of inputs. But it's slower than `x & y`, which doesn't modify in place and uses more mem.
    for i = 1:length(x)
        @inbounds x[i] &= y[i]
    end
end

n = 100_000_000
@code_warntype bitwise_and!(bitrand(n), bitrand(n))
@time bitwise_and!(bitrand(n), bitrand(n))

function inplace_bitwise_and!(res::BitArray, x::BitArray, y::BitArray)
    for i = 1:length(res)
        res[i] = x[i] & y[i]
    end
end

@code_warntype inplace_bitwise_and!(falses(n), bitrand(n), bitrand(n))
@time inplace_bitwise_and!(falses(n), bitrand(n), bitrand(n))

function transactions_to_nodes(T::Array{Array{String,1},1})
    n = length(T)
    uniq_items = get_unique_items(T)
    sort!(uniq_items)
    p = length(uniq_items)
    occur = falses(n, p)

    for j = 1:p
        @simd for i = 1:n
            @inbounds occur[i, j] = uniq_items[j] ∈ T[i]
        end
    end
    nodes = Array{Node,1}(p)
    for j = 1:p
        @inbounds nodes[j] = Node([Int16(j)], occur[:, j])
    end

    return nodes
end

# Get size of largest transaction
function max_transaction(T::Array{Array{String,1},1})
    res = 0
    n = length(T)
    for i = 1:n
        len = length(T[i])
        res = (len > res) ? len : res
    end
    res
end

t = [["a", "b"],
     ["b", "c", "d"],
     ["a", "c"],
     ["e", "b"],
     ["a", "c", "d"],
     ["a", "e"],
     ["a", "b", "c"],
     ["c", "b", "e"]]

@code_warntype transactions_to_nodes(t)

@time transactions_to_nodes(t)

n = 100_000
transactions = [sample(groceries, 20, replace = false) for x in 1:n];
@time w1 = transactions_to_nodes(transactions);

function merge_nodes(node1, node2, k)
    ids = Array{Int,1}(k)
    ids[1:k-1] = deepcopy(node1.item_ids[1:(k-1)])
    if k == 2
        ids[k] = node2.item_ids[1]
    elseif k > 2
        ids[k] = node2.item_ids[k-1]
    end
    transacts = node1.transacts & node2.transacts
    nd = Node(ids, transacts)
    return nd
end

function gen_next_layer(prev::Array{Node,1}, minsupp = 1)
    if length(prev) == 0
        n = 0
    else
        k = length(prev[1].item_ids) + 1
        n = length(prev)
    end
    nodes = Array{Node,1}(0)             # next layer of nodes

    for i = 1:(n-1)
        for j = (i+1):n
            if k == 2 || prefix(prev[i], k) == prefix(prev[j], k)
                nd = merge_nodes(prev[i], prev[j], k, n_obs)
                if sum(nd.transacts) ≥ minsupp
                    push!(nodes, nd)
                end
            end
        end
    end
    nodes
end

a1 = transactions_to_nodes(t)
@code_warntype gen_next_layer(a1)

n = 100_000
transactions = [sample(groceries, 12, replace = false) for x in 1:n]
@time w1 = transactions_to_nodes(transactions);
@time gen_next_layer(w1);
a2 = gen_next_layer(a1)
a3 = gen_next_layer(a2)
a4 = gen_next_layer(a3)

function frequent(T::Array{Array{String,1},1}, minsupp = 0)
    nodes = transactions_to_nodes(T)
    max_items = max_transaction(T)

    F = Array{Array{Node,1},1}(max_items)
    F[1] = nodes
    k = 2

    while k <= max_items
        # println(k)
        F[k] = gen_next_layer(F[k-1], minsupp)
        k += 1
    end
    F
end

n = 100_000
t = [sample(groceries, 50, replace = false) for _ in 1:n];

@code_warntype frequent(t, 1)
@time f = frequent(t, round(Int, n*0.2));


text:
"""
`JuLIP.FIO` : provides some basis file IO. All sub-modules and derived
packages should use these interface functions so that the file formats
can be changed later. This submodule provides

- load_dict
- save_dict
- decode_dict
"""
module FIO

using JSON

export load_dict, save_dict, save_json, load_json,
       decode_dict, read_dict, write_dict

#######################################################################
#                     Conversions to and from Dict
#######################################################################

# eventually we want to be able to serialise all JuLIP types to Dict
# and back and those Dicts may only contain elementary data types
# this will then allow us to load them back via decode_dict without
# having to know the type in the code

"""
`decode_dict(D::Dict) -> ?`

Looks for a key `__id__` in `D` whose value is used to dynamically dispatch
the decoding to
```julia
convert(Val(Symbol(D["__id__"])))
```
That is, a user defined type must implement this `convert(::Val, ::Dict)`
utility function. The typical code would be
```julia
module MyModule
    struct MyStructA
        a::Float64
    end
    Dict(A::MyStructA) = Dict( "__id__" -> "MyModule_MyStructA",
                               "a" -> a )
    MyStructA(D::Dict) = A(D["a"])
    Base.convert(::Val{:MyModule_MyStructA})
end
```
The user is responsible for choosing an id that is sufficiently unique.

The purpose of this function is to enable the loading of more complex JSON files
and automatically converting the parsed Dict into the appropriate composite
types. It is intentional that a significant burden is put on the user code
here. This will maximise flexibiliy, e.g., by introducing version numbers,
and being able to read old version in the future.
"""
function decode_dict(D::Dict)
    if !haskey(D, "__id__")
        error("JuLIP.IO.decode_dict: `D` has no key `__id__`")
    end
    return read_dict(Val(Symbol(D["__id__"])), D)
end

read_dict(D::Dict) = decode_dict(D)

read_dict(v::Val, D::Dict) = convert(v, D)

#######################################################################
#                     JSON
#######################################################################

function load_dict(fname::AbstractString)
    return JSON.parsefile(fname)
end

function save_dict(fname::AbstractString, D::Dict; indent=2)
    f = open(fname, "w")
    JSON.print(f, D, indent)
    close(f)
    return nothing
end

save_json = save_dict
load_json = load_dict

## Some useful utility functions

write_dict(A::Matrix{T}) where {T <: Number} =
    Dict("__id__" => "JuLIP_Matrix",
         "T"      => string(T),
         "nrows"  => size(A, 1),
         "ncols"  => size(A, 2),
         "data"   => A[:])

function read_dict(::Val{:JuLIP_Matrix}, D::Dict)
   T = Meta.eval(Meta.parse(D["T"]))
   return resize(T.(D["data"]), D["ncols"], D["nrows"])
end

##
end


text:
using CLIMA
using MPI
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.VTK: writemesh
using Logging
using LinearAlgebra
using Random
using StaticArrays
using CLIMA.DGmethods:
    DGModel,
    Vars,
    vars_state_conservative,
    number_state_conservative,
    init_ode_state
using CLIMA.ColumnwiseLUSolver: banded_matrix, banded_matrix_vector_product!
using CLIMA.DGmethods.NumericalFluxes:
    RusanovNumericalFlux,
    CentralNumericalFluxSecondOrder,
    CentralNumericalFluxGradient
using CLIMA.MPIStateArrays: MPIStateArray, euclidean_distance

using Test

include("../DGmethods/advection_diffusion/advection_diffusion_model.jl")

struct Pseudo1D{n, α, β, μ, δ} <: AdvectionDiffusionProblem end

function init_velocity_diffusion!(
    ::Pseudo1D{n, α, β},
    aux::Vars,
    geom::LocalGeometry,
) where {n, α, β}
    # Direction of flow is n with magnitude α
    aux.u = α * n

    # diffusion of strength β in the n direction
    aux.D = β * n * n'
end

function initial_condition!(
    ::Pseudo1D{n, α, β, μ, δ},
    state,
    aux,
    x,
    t,
) where {n, α, β, μ, δ}
    ξn = dot(n, x)
    # ξT = SVector(x) - ξn * n
    state.ρ = exp(-(ξn - μ - α * t)^2 / (4 * β * (δ + t))) / sqrt(1 + t / δ)
end

let
    # boiler plate MPI stuff
    CLIMA.init()
    ArrayType = CLIMA.array_type()

    mpicomm = MPI.COMM_WORLD
    Random.seed!(777 + MPI.Comm_rank(mpicomm))

    # Mesh generation parameters
    N = 4
    Nq = N + 1
    Neh = 10
    Nev = 4

    @testset "$(@__FILE__) DGModel matrix" begin
        for FT in (Float64, Float32)
            for dim in (2, 3)
                for single_column in (false, true)
                    # Setup the topology
                    if dim == 2
                        brickrange = (
                            range(FT(0); length = Neh + 1, stop = 1),
                            range(FT(1); length = Nev + 1, stop = 2),
                        )
                    elseif dim == 3
                        brickrange = (
                            range(FT(0); length = Neh + 1, stop = 1),
                            range(FT(0); length = Neh + 1, stop = 1),
                            range(FT(1); length = Nev + 1, stop = 2),
                        )
                    end
                    topl = StackedBrickTopology(mpicomm, brickrange)

                    # Warp mesh
                    function warpfun(ξ1, ξ2, ξ3)
                        # single column currently requires no geometry warping

                        # Even if the warping is in only the horizontal, the way we
                        # compute metrics causes problems for the single column approach
                        # (possibly need to not use curl-invariant computation)
                        if !single_column
                            ξ1 = ξ1 + sin(2π * ξ1 * ξ2) / 10
                            ξ2 = ξ2 + sin(2π * ξ1) / 5
                            if dim == 3
                                ξ3 = ξ3 + sin(8π * ξ1 * ξ2) / 10
                            end
                        end
                        (ξ1, ξ2, ξ3)
                    end

                    # create the actual grid
                    grid = DiscontinuousSpectralElementGrid(
                        topl,
                        FloatType = FT,
                        DeviceArray = ArrayType,
                        polynomialorder = N,
                        meshwarp = warpfun,
                    )
                    d = dim == 2 ? FT[1, 10, 0] : FT[1, 1, 10]
                    n = SVector{3, FT}(d ./ norm(d))

                    α = FT(1)
                    β = FT(1 // 100)
                    μ = FT(-1 // 2)
                    δ = FT(1 // 10)
                    model = AdvectionDiffusion{dim}(Pseudo1D{n, α, β, μ, δ}())

                    # the nonlinear model is needed so we can grab the state_auxiliary below
                    dg = DGModel(
                        model,
                        grid,
                        RusanovNumericalFlux(),
                        CentralNumericalFluxSecondOrder(),
                        CentralNumericalFluxGradient(),
                    )

                    vdg = DGModel(
                        model,
                        grid,
                        RusanovNumericalFlux(),
                        CentralNumericalFluxSecondOrder(),
                        CentralNumericalFluxGradient();
                        direction = VerticalDirection(),
                        state_auxiliary = dg.state_auxiliary,
                    )

                    A_banded = banded_matrix(
                        vdg,
                        MPIStateArray(dg),
                        MPIStateArray(dg);
                        single_column = single_column,
                    )

                    Q = init_ode_state(dg, FT(0))
                    dQ1 = MPIStateArray(dg)
                    dQ2 = MPIStateArray(dg)

                    vdg(dQ1, Q, nothing, 0; increment = false)
                    Q.data .= dQ1.data

                    vdg(dQ1, Q, nothing, 0; increment = false)
                    banded_matrix_vector_product!(vdg, A_banded, dQ2, Q)
                    @test all(isapprox.(
                        Array(dQ1.realdata),
                        Array(dQ2.realdata),
                        atol = 100 * eps(FT),
                    ))

                    α = FT(1 // 10)
                    function op!(LQ, Q)
                        vdg(LQ, Q, nothing, 0; increment = false)
                        @. LQ = Q + α * LQ
                    end

                    A_banded = banded_matrix(
                        op!,
                        vdg,
                        MPIStateArray(dg),
                        MPIStateArray(dg);
                        single_column = single_column,
                    )

                    Q = init_ode_state(dg, FT(0))
                    dQ1 = MPIStateArray(vdg)
                    dQ2 = MPIStateArray(vdg)

                    op!(dQ1, Q)
                    Q.data .= dQ1.data

                    op!(dQ1, Q)
                    banded_matrix_vector_product!(vdg, A_banded, dQ2, Q)
                    @test all(isapprox.(
                        Array(dQ1.realdata),
                        Array(dQ2.realdata),
                        atol = 100 * eps(FT),
                    ))
                end
            end
        end
    end
end

nothing


text:
using Gen
using GenDirectionalStats
using Test
using LinearAlgebra: norm
import Rotations
using SpecialFunctions: besseli

@testset "hopf fibration" begin

    # identity rotation

    (direction, rot2) = to_direction_and_plane_rotation(Rot3(1.0, 0.0, 0.0, 0.0))
    @test isapprox(Rotations.rotation_angle(rot2), 0.0, atol=1e-5)
    @test isapprox(direction, UnitVector3(0.0, 0.0, 1.0))

    rot3 = from_direction_and_plane_rotation(UnitVector3(0.0, 0.0, 1.0), Rotations.RotMatrix{2,Float64}(0.0))
    @test isapprox(rot3, Rot3(1.0, 0.0, 0.0, 0.0))

    # check that mapping between (UnitVector3, PlaneRotation) parametrization
    # and Rotation3D is an involution

    for i in 1:20
        rot = uniform_rot3()
        (direction, plane_rotation) = to_direction_and_plane_rotation(rot)
        new_rot = from_direction_and_plane_rotation(direction, plane_rotation)
        @test isapprox(new_rot, rot)
    end

   for i in 1:20
        direction = uniform_3d_direction()
        plane_rotation = uniform_rot2()
        rot = from_direction_and_plane_rotation(direction, plane_rotation)
        (new_direction, new_plane_rotation) = to_direction_and_plane_rotation(rot)
        @test isapprox(new_direction, direction)
        @test isapprox(new_plane_rotation, plane_rotation)
    end

end

function rand_axis()
    x, y, z = rand(3)
    axis = [x, y, z]
    axis / norm(axis)
end

@gen function foo()
    # uniform distribution, so all moves should be accepted
    rot ~ uniform_rot3()
end

@testset "uniform_angle_fixed_axis_mh" begin

    # smoke test (these should all accept)
    trace = simulate(foo, ())
    for i=1:10
        @test uniform_angle_fixed_axis_mh(trace, :rot, rand_axis();
            check=true, egocentric=false)[2]
        @test uniform_angle_fixed_axis_mh(trace, :rot, rand_axis();
            check=true, egocentric=true)[2]
    end
end

@testset "flip_around_fixed_axis_mh" begin

    # smoke test
    for i=1:10
        trace = simulate(foo, ())
        trace, acc = flip_around_fixed_axis_mh(trace, :rot, rand_axis(), check=true)
        @test acc
    end

    # actual test (since this move is deterministic)
    # test a flip around the y-axis
    rot = Rot3(Rotations.AngleAxis(pi/4, 1.0, 0.0, 0.0))
    trace, = generate(foo, (), choicemap((:rot, rot)))
    expected = Rot3(
        [-1.0 0.0 0.0;
        0.0 cos(pi/4) sin(pi/4);
        0.0 sin(pi/4) -cos(pi/4)])

    # allocentric
    axis = [0., cos(pi/4), sin(pi/4)]
    new_tr, = flip_around_fixed_axis_mh(trace, :rot, axis; check=true, egocentric=false)
    @test isapprox(new_tr[:rot], expected)

    # egocentric
    axis = [0., 1., 0.]
    new_tr, = flip_around_fixed_axis_mh(trace, :rot, axis; check=true, egocentric=true)
    @test isapprox(new_tr[:rot], expected)
end

@testset "small_angle_fixed_axis_mh" begin

    # smoke test (these should all accept)
    trace = simulate(foo, ())
    for i=1:10
        @test small_angle_fixed_axis_mh(trace, :rot, rand_axis(), 0.1;
            check=true, egocentric=false)[2]
        @test small_angle_fixed_axis_mh(trace, :rot, rand_axis(), 0.1;
            check=true, egocentric=true)[2]
    end

end

@testset "small_angle_random_axis_mh" begin

    # smoke test (these should all accept)
    trace = simulate(foo, ())
    for i=1:10
        @test small_angle_random_axis_mh(trace, :rot, 0.1;
            check=true, egocentric=false)[2]
        @test small_angle_random_axis_mh(trace, :rot, 0.1;
            check=true, egocentric=true)[2]
    end
end

@testset "log besseli 1" begin

    actual = GenDirectionalStats.log_besseli(1, 300.0; n_terms=1000)
    expected = log(besseli(1, 300.0))
    @test isapprox(actual, expected)

    function test_logbesseli_2(x)
        approx = Inf
        prev_approx = -Inf
        n_terms = 1

        while abs(approx - prev_approx) > 0.01
            prev_approx = approx
            approx = GenDirectionalStats._log_besseli(1, x; n_terms=n_terms)
            n_terms += 1
        end
        return n_terms
    end

    x = 10000
    n_terms = test_logbesseli_2(x)
    @show n_terms
    @show n_terms / x
    @time GenDirectionalStats.log_besseli(1, 10)
end


text:
#=

  Some examples of the regular global constraint in Julia ConstraintSolver.jl

  Here's a description of the transition tables.
  Note that input can never be 0 since it is used for no state.

  Example: regexp 123*21

     Transition = [
       # inputs: 1  2  3
                [2, 0, 0], # transitions from state 1: -2-> state 2  [start]
                [0, 3, 0], # transitions from state 2: -3-> state 3 #
                [0, 4, 3], # transitions from state 3: -3-> state 3, -2-> state 4
                [5, 0, 0], # transitions from state 4: -1-> state 5  [final]
                [0, 0, 0]  # transitions from state 5: (none)
                ],
     InitialState = 1,
     AcceptingStates = [5]

  (This is the first example, :1, below.)

  The states are numbered from 1..Transition.length.
  The initial state is in InitialState, i.e. the state Transition[InitialState].
  The columns are the inputs (the column index), e.g. the first state ([2, 0, 0])
  means that it only accept a "1" (first column), and then go to state 2 ([0, 3, 0],
  which only accept a "2" (second column) and go to state 3 ([0, 4, 3]) with the
  following meaning:

      - if "1": not accepted
      - if "2": goto state 4
      - if "3": goto state 3 (loop)

  And the DFA now go through the states. The accepting state in this
  example is only state 5 (AcceptingStates must be a list).

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# the regular constraints is in constraints_utils.jl
#
function regular_test(problem,n=4,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    transition       = problem[:transition]
    n_states          = problem[:n_states]
    input_max        = problem[:input_max]
    initial_state    = problem[:initial_state]
    accepting_states = problem[:accepting_states]

    @variable(model, 1 <= x[1:n] <= 3, Int)

    a = regular(model,x,n_states,input_max,transition,initial_state,accepting_states)

    # @objective(model,Max,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                a_val = convert.(Integer,JuMP.value.(a; result=sol))
                println("x:$x_val")
                println("a:$a_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

regular_problems = Dict(

#
# regexp 123*21
#
:1 => Dict(
   :info => "regex: 123*21",
   :transition => resize_matrix([
                [2, 0, 0], # transitions from state 1: -2-> state 2
                [0, 3, 0], # transitions from state 2: -3-> state 3 #
                [0, 4, 3], # transitions from state 3: -3-> state 3, -2-> state 4
                [5, 0, 0], # transitions from state 4: -1-> state 5
                [0, 0, 0]  # transitions from state 5: (none)
                ]),
   :n_states => 5,
   :input_max => 3,
   :initial_state => 1,
   :accepting_states => [5],
   :lens => 4:10, # The lengths to check ('n' in the model)
   :all_solutions => true, # Find all solutions?
   ),

# This example is from
# "Handbook of Constraint programming", page 180
#
#    ^(c+|a+b+a+){n}$
#    ^(3+|1+2+1+){n}$
#
# Or:  aa*bb* + cc*
#
# "It accepts the strings 'aaabaa' and 'cc', but not 'aacbba'."
#
:2 => Dict(
  :info => "It accepts the strings 'aaabaa' and 'cc', but not 'aacbba'.",
  :transition => resize_matrix([
      # inputs
      # 1 2 3
       [2,0,5], # 1 start
       [2,3,0], # 2 loop + next
       [4,3,0], # 3 next + loop
       [4,0,0], # 4 loop + end state
       [0,0,5]  # 5 loop and end state
     ]),
  :n_states => 5,
  :input_max => 3,
  :accepting_states => [4,5],
  :initial_state => 1,
  :lens => [4,5],
  :all_solutions => true,
 ),

#
# Example1 from Choco http://www.emn.fr/x-info/choco-solver/doku.php?id=regular
#
# First example
#        t.add(new Transition(0, 1, 1));
#        t.add(new Transition(1, 1, 2));
#        t.add(new Transition(2, 1, 3));
#
#        t.add(new Transition(3, 3, 0));
#        t.add(new Transition(0, 3, 0));
#
:3  => Dict(
  :info => "From Choco I",
  :transition => resize_matrix([
    [2, 0, 1], # 1
    [3, 0, 0], # 2
    [0, 4, 0], # 3
    [0, 0, 1]  # 4
   ]),
  :n_states => 4,
  :input_max => 3,
  :initial_state => 1,
  :accepting_states => [4],
  :lens => [10],
  :all_solutions => true,
),

#
# Example2 from Choco http://www.emn.fr/x-info/choco-solver/doku.php?id=regular
#
# This is an all_equal constraint. The initial_state decides which
# values it will be

:4 => Dict(
  :info => "From Choco 2. All equal constraint. initial_state (here 3) decides the value.",
  :transition => resize_matrix([
     [1, 0, 0, 0], # 1 (start)
     [0, 2, 0, 0], # 2
     [0, 0, 3, 0], # 3
     [0, 0, 0, 4]  # 4
   ]),
  :n_states => 4,
  :input_max => 4,
  :initial_state => 3,
  :accepting_states => 1:4,
  :lens => 3:6,
  :all_solutions => true,

),

#
# Example3 from Choco http://www.emn.fr/x-info/choco-solver/doku.php?id=regular
#
# This is stretchPath, not regular
#  lgt.add(new int[]{2, 2, 2}); // stretches of value 2 are exactly of size 2
#  lgt.add(new int[]{0, 2, 2}); // stretches of value 0 are exactly of size 2
#  lgt.add(new int[]{1, 2, 3}); // stretches of value 1 are at least of size 2 and at most 3
#
# This should be interpreted as:
#  0{2}1{2,3}2{2}
# for all permutations of 0,1,and 2.
# Or rather:
#  1{2}2{2,3}3{2}
# for all permutations of 1,2, and 3 since 0 is not allowed as an input value...
#
:5 => Dict(
  :info => "stretch_path: 1{2}2{2,3}3{2} for all permutations of 1,2, and 3",
  :transition => resize_matrix([
     [2, 4, 7], # 1 (start)
     [3, 0, 0], # 2
     [0, 4, 7], # 3 (end)
     [0, 5, 0], # 4
     [2, 6, 7], # 5 (end)
     [2, 0, 7], # 6 (end)
     [0, 0, 8], # 7
     [2, 4, 0]  # 8 (end)
   ]),
  :n_states => 8,
  :input_max => 3,
  :initial_state => 1,
  :accepting_states => [3,5,6,8],
  :lens => 4:6,
  :all_solutions => true,
),

#
# All different.
# This works in principle (here is n=3), but the automaton will
# be forbiddingly large for larger arrays.
#
:6 => Dict(
  :info => "all_different (size 3)",
  :transition => resize_matrix([
     #  1   2   3
      [ 2,  7, 12], # 1 (start)
      [ 0,  3,  5], # 2
      [ 0,  0,  4], # 3
      [ 0,  0,  0], # 4
      [ 0,  6,  0], # 5
      [ 0,  0,  0], # 6
      [ 8,  0, 10], # 7
      [ 0,  0,  9], # 8
      [ 0,  0,  0], # 9
      [11,  0,  0], # 10
      [ 0,  0,  0], # 11
      [13, 15,  0], # 12
      [ 0, 14,  0], # 13
      [ 0,  0,  0], # 14
      [ 16, 0,  0], # 15
      [ 0,  0,  0]  # 16
   ]),
  :n_states => 16,
  :input_max => 3,
  :initial_state => 1,
  :accepting_states => [4,6,9,11,14,16],
  :lens => [3],
  :all_solutions => true,
 ),

#
# "Need Regular Expression for Finite Automata: Even number of 1s and Even number of 0s"
# http://stackoverflow.com/questions/17420332/need-regular-expression-for-finite-automata-even-number-of-1s-and-even-number-o
#
#     Q1 <- 1 -> Q2
#     ^          ^
#     |          |
#     2          2
#     |          |
#     v          v
#     Q4 <- 1 -> Q3
#
:7 => Dict(
  :info => "Even number of 1s and even number of 0s",
  :transition => resize_matrix([
     #  1   2
      [ 2,  4], # state 1 Q1 start, final
      [ 1,  3], # state 2 Q2
      [ 4,  2], # state 3 Q3
      [ 3,  1]  # state 4 Q4
   ]),
  :n_states => 4,
  :input_max => 2,
  :initial_state => 1,
  :accepting_states => [1],
  :lens => 1:6, # must be an even number
  :all_solutions=> true
  ),

)

# Test the problems
for p in sort(collect(keys(regular_problems)))
    println("\n\nPROBLEM $p")
    problem = regular_problems[p]
    local all_solutions = get(problem,:all_solutions,false)
    lens = get(problem,:lens, 5)
    info = get(problem,:info, "<no info")
    println("Info: $info")
    for n in lens
        println("\nn:$n")
        @time regular_test(regular_problems[p],n,true,all_solutions)
    end
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using LocalHillResponses
using Test

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


text:
module Generative

using Flux, Distributions
using Flux: glorot_uniform, @functor
using Random
import Random: GLOBAL_RNG
import Distributions: logpdf

include("vae.jl")
include("data_loader/data_loader.jl")

end # module

max_stars_repo_name: HamletWantToCode/Generative

max_issues_repo_name: HamletWantToCode/Generative

max_forks_repo_name: HamletWantToCode/Generative


text:
#####
##### Centered second-order advection scheme
#####

struct CenteredSecondOrder <: AbstractAdvectionScheme end

boundary_buffer(::CenteredSecondOrder) = 0

const C2 = CenteredSecondOrder

@inline momentum_flux_uu(i, j, k, grid, ::C2, U, u) = ℑxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, U) * ℑxᶜᵃᵃ(i, j, k, grid, u)
@inline momentum_flux_uv(i, j, k, grid, ::C2, V, u) = ℑxᶠᵃᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, V) * ℑyᵃᶠᵃ(i, j, k, grid, u)
@inline momentum_flux_uw(i, j, k, grid, ::C2, W, u) = ℑxᶠᵃᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, W) * ℑzᵃᵃᶠ(i, j, k, grid, u)

@inline momentum_flux_vu(i, j, k, grid, ::C2, U, v) = ℑyᵃᶠᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, U) * ℑxᶠᵃᵃ(i, j, k, grid, v)
@inline momentum_flux_vv(i, j, k, grid, ::C2, V, v) = ℑyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, V) * ℑyᵃᶜᵃ(i, j, k, grid, v)
@inline momentum_flux_vw(i, j, k, grid, ::C2, W, v) = ℑyᵃᶠᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, W) * ℑzᵃᵃᶠ(i, j, k, grid, v)

@inline momentum_flux_wu(i, j, k, grid, ::C2, U, w) = ℑzᵃᵃᶠ(i, j, k, grid, Ax_ψᵃᵃᶠ, U) * ℑxᶠᵃᵃ(i, j, k, grid, w)
@inline momentum_flux_wv(i, j, k, grid, ::C2, V, w) = ℑzᵃᵃᶠ(i, j, k, grid, Ay_ψᵃᵃᶠ, V) * ℑyᵃᶠᵃ(i, j, k, grid, w)
@inline momentum_flux_ww(i, j, k, grid, ::C2, W, w) = ℑzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, W) * ℑzᵃᵃᶜ(i, j, k, grid, w)

# Calculate the flux of a tracer quantity c through the faces of a cell.
# In this case, the fluxes are given by u*Ax*c̄ˣ, v*Ay*c̄ʸ, and w*Az*c̄ᶻ.
@inline advective_tracer_flux_x(i, j, k, grid, ::C2, U, c) = Ax_uᶠᶜᶜ(i, j, k, grid, U) * ℑxᶠᵃᵃ(i, j, k, grid, c)
@inline advective_tracer_flux_y(i, j, k, grid, ::C2, V, c) = Ay_vᶜᶠᶜ(i, j, k, grid, V) * ℑyᵃᶠᵃ(i, j, k, grid, c)
@inline advective_tracer_flux_z(i, j, k, grid, ::C2, W, c) = Az_wᶜᶜᵃ(i, j, k, grid, W) * ℑzᵃᵃᶠ(i, j, k, grid, c)


text:
export references

"""
    references(dist::SemiMetric, db; <kwargs>) -> SubDatabase

Computes a set of references, a sample of `db`, using the [`KCenters`](@ref) specification, it is a wrapper to `kcenters` function in the `KCenters.jl` package.
The set of references are computed taking into account the specified `r` but also the distance function `dist`.

# Arguments

- `dist`: a distance function
- `db`: the database to be sampled

# Keyword arguments
- `k::Int`: the number of centers to compute, defaults to `sqrt(|db|)`
- `sample::Real`: indicates the sampling size before computing the set of `k` references, defaults to `log(|db|) k`; `sample=0` means for no sampling.
- `maxiters::Int`: number of iterationso  of the Lloyd algorithm that should be applied on the initial computation of centers, that is, `maxiters > 0` applies `maxiters` iterations of the algorithm.
- `tol::Float64`: change tolerance to stop the Lloyd algorithm (error changes smaller than `tol` among iterations will stop the algorithm)
- `initial`: initial centers or a strategy to compute initial centers, symbols `:rand`, `:fft`, and `:dnet`.
There are several interactions between initial values and parameter interactions (described in `KCenters` object), for instance,
the `maxiters > 0` apply the Lloyd's algorithm to the initial computation of references.

- if `initial=:rand`:
  - `maxiters = 0` will retrieve a simple random sampling
  - `maxiters > 0' achieve kmeans-centroids, `maxiters` should be set appropiately for the the dataset
- if `initial=:dnet`:
  - `maxiters = 0` computes a pure density-net
  - `maxiters > 0` will compute a kmeans centroids but with an initialization based on the dnet
- if `initial=:fft`:
  - `maxiters = 0` computes `k` centers with the farthest first traversal algorithm
  - `maxiters > 0` will use the FFT based kcenters as initial points for the Lloyd algorithm

Note 1: `maxiters > 0` needs to compute centroids and these centroids should be _defined_
for the specific data model, and also be of use in the specific metric distance and error function.

Note 2: The error function is defined as the mean of distances of all objects in `db` to its associated nearest centers in each iteration.

Note 3: The computation of references on large databases can be prohibitive, in these cases please consider to work on a sample of `db`
"""
function references(
      dist::SemiMetric, db;
      k=ceil(Int, sqrt(length(db))),
      sample=ceil(Int, log2(length(db)) * k),
      maxiters=0,
      tol=0.001, initial=:rand)
    0 < k < length(db) || throw(ArgumentError("invalid relation between k and n, must follow 0 < k < n"))
    C = if sample > 0
      s = unique(rand(eachindex(db), sample))
      kcenters(dist, SubDatabase(db, s), k; initial, maxiters, tol)
    else
      kcenters(dist, db, k; initial, maxiters, tol)
    end

    C.centers[C.dmax .> 0.0]  # centers covering more than itself
end

max_stars_repo_name: sadit/NeighborhoodApproximationIndex

max_issues_repo_name: sadit/NeighborhoodApproximationIndex

max_forks_repo_name: sadit/NeighborhoodApproximationIndex


text:
using MultiResolutionIterators

animal_info = [
    [["Turtles", "are", "reptiles", "."],
     ["They", "have", "shells", "."],
     ["They", "live", "in", "the", "water", "."]],
    [["Cats", "are", "mammals", "."],
     ["They", "live", "on", "the", "internet", "."]]
    ]

# Get rid of document boundaries
flatten_levels(animal_info, 1) |> full_consolidate

# Get rid of sentence boundaries, so documents made up of words
flatten_levels(animal_info, 2) |> full_consolidate

# Get rid of document and sentence boundries
flatten_levels(animal_info, 1:2) |> full_consolidate

# Get rid of all boundaries, just a stream of characters
flatten_levels(animal_info, ALL_LEVELS) |> full_consolidate

# Get rid of word boundaries so each document is a a stream of characters
flatten_levels(animal_info, [1,3]) |> full_consolidate

# Join all words using spaces, keep other structure
join_levels(animal_info, Dict(3=>" ")) |> full_consolidate

### Using Named levels

# Declare an indexer.
struct AnimalTextIndexer end;

# Overload `levelname_map` this so it knows the name mapping
MultiResolutionIterators.levelname_map(::AnimalTextIndexer) = [
    :documents=>1,
    :sentences=>2,
    :words=>3, :tokens=>3, # can have multiple aliases for same level
    :characters=>4 # As characters themselves are not iterable this name/level has little effect
]
indexer = AnimalTextIndexer();

# Merge all sentences
flatten_levels(animal_info, lvls(indexer, :sentences)) |> full_consolidate

# I.e keep documents, and words
flatten_levels(animal_info, (!lvls)(indexer, :words, :documents)) |> full_consolidate

# Merge everything **except** words
flatten_levels(animal_info, (!lvls)(indexer, :words)) |> full_consolidate

# Merge everything **except** words and sentences
flatten_levels(animal_info, (!lvls)(indexer, :words, :sentences)) |> full_consolidate

# i.e. merge documents
flatten_levels(animal_info, lvls(indexer, :documents)) |> full_consolidate

# # Join all words using spaces, join all sentences with new lines, all documents with double new lines
join_levels(animal_info,
 lvls(indexer,Dict(:words=>" ", :sentences=>"\n", :documents=>"\n---\n"))) |>
 full_consolidate |> print


text:
using ShaleDrillingLikelihood: vw_revenue, vw_cost, vw_extend, vw_scrap,
    DrillingRevenueUnconstrained, DrillingRevenueConstrained,
    ConstrainedIdx, UnconstrainedFmConstrainedIdx

export updateThetaUnconstrained!, ThetaConstrained, ThetaRho

ThetaRho() = 0.726029
AlphaPsi() = 0x1.55b10dcc51d5fp-2
AlphaG() = 0x1.075c51cf96449p-1
AlphaT() = STARTING_α_t
Alpha0() = -2.8
AlphaD() = -0.2
BetaPsi() = 0.11805850128182346
AlphaTFE() = [ 0.20665908809883377, 0.1661000339389365, 0.14007105659535182,
    0.14457016546750073, 0.19677616668888495, 0.26422850311680984,
    0.5206302285265718, 0.635563439673686,]
Sigma_eta() = 0x1.8c7275fb0e3d2p-4
Sigma_u() = 0x1.2ca92ef5843d8p-2
Gamma0() = -0x1.c98f8c57dd39bp+3

@deprecate SigmaSq_eta() Sigma_eta()
@deprecate SigmaSq_u() Sigma_u()

Theta(m, args...) = throw(error("Default Theta not assigned for $m"))
Theta(m::AbstractDynamicDrillModel, args...; kwargs...) = Theta(reward(m), args...; kwargs...)
function Theta(m::DrillReward, args...; kwargs...)
    c = Theta(cost(m), args...; kwargs...)
    e = Theta(extend(m), args...; kwargs...)
    s = Theta(scrap(m), args...; kwargs...)
    r = Theta(revenue(m), args...; kwargs...)
    return vcat(c,e,s,r)
end

function ThetaConstrained(m::DrillReward{<:DrillingRevenueUnconstrained}, theta)
    rwrd_c = ConstrainedProblem(m, theta)
    theta_rev = vw_revenue(m, theta)
    idxit = Base.Iterators.flatten(ConstrainedIdx(revenue(m)))
    r = [x for (i,x) in enumerate(theta_rev) if i ∉ idxit]
    c = vw_cost(m,theta)
    e = vw_extend(m,theta)
    s = vw_scrap(m,theta)
    return vcat(c,e,s,r)
end

function updateThetaUnconstrained!(m::DrillReward{<:DrillingRevenueConstrained}, thetau, thetac)
    rwrd_u = UnconstrainedProblem(m)
    _nparm(m) == length(thetac) ||
        throw(DimensionMismatch("_nparm(m) = $(_nparm(m)) != length(thetac) = $(length(thetac))"))
    _nparm(rwrd_u) == length(thetau) ||
        throw(DimensionMismatch("_nparm(rwrd_u) = $(_nparm(rwrd_u)) != length(thetau) = $(length(thetau))"))
    vw_cost(rwrd_u, thetau) .= vw_cost(m, thetac)
    vw_extend(rwrd_u, thetau) .= vw_extend(m, thetac)
    vw_scrap(rwrd_u, thetau) .= vw_scrap(m, thetac)

    revidx = UnconstrainedFmConstrainedIdx(revenue(m))
    thetarev_u = vw_revenue(rwrd_u, thetau)
    thetarev_u[revidx] .= vw_revenue(m, thetac)
    return thetau
end

function updateThetaUnconstrained!(m::DrillReward{<:DrillingRevenueUnconstrained}, thetau, thetac)
    rwrd_c = ConstrainedProblem(m)
    return updateThetaUnconstrained!(rwrd_c, thetau, thetac)
end

# Royalty
function Theta(m::RoyaltyModel; θρ=ThetaRho(), kwargs...)
    #            θρ  ψ        g    x_r      κ₁   κ₂
    return vcat(θρ, 1.0,    0.1, -0.3,    -0.6, 0.6)
end

# Production
function Theta(m::ProductionModel; αψ=AlphaPsi(), αg=AlphaG(), kwargs...)
    #           αψ, αg, γx   σ2η, σ2u  # NOTE: η is iwt, u is iw
    return vcat(αψ, αg, 0.2, 0.3, 0.4)
end

# Test Drill
function Theta(m::TestDrillModel; θρ=ThetaRho(), αψ=AlphaPsi(), kwargs...)
    return vcat(αψ, 2.0, -2.0, -0.75, θρ)
end

# from datasets
# ---------

function Theta(d::DataRoyalty; θρ=ThetaRho(), αψ=AlphaPsi(), kwargs...)
    βψ = BetaPsi()
    if num_choices(d) == 6 && _num_x(d) == 4
        betas = [ 0.5966055028536498, 1.1822487035901503, -1.7017414563283797, 0.14048454874430671]
        kappas = [3.8373126232392716, 4.172804362649826, 5.0264099395765385, 5.925934015361181, 6.501397017005945]
        return vcat(θρ, βψ, betas, kappas)
    end
    throw(error("no starting values available"))
end

function Theta(d::DataProduce; αψ=AlphaPsi(), αg=AlphaG(), αt=AlphaT(),
         αTFE=AlphaTFE(), ση=Sigma_eta(), σu=Sigma_u(), γ0=Gamma0(), αD=AlphaD(), kwargs...)
    if _num_x(d) == 2
        return vcat(αψ, αg, γ0,     ση, σu)
    elseif _num_x(d) == 3
        return vcat(αψ, αg, αt, γ0, ση, σu)
    elseif _num_x(d) == 4
        return vcat(αψ, αg, αt, αD, γ0, ση, σu)
    elseif _num_x(d) == 10
        return vcat(αψ, αg, αTFE, γ0, ση, σu)
    end
    throw(error("no starting values available"))
end

# revenue
# ---------
function Theta(m::DrillingRevenue{Constrained}; θρ=ThetaRho(), kwargs...)
    return vcat(Alpha0(), θρ)
end

function Theta(m::DrillingRevenue{Unconstrained, NoTrend};
    θρ=ThetaRho(), αψ=AlphaPsi(), αg=AlphaG(), α0=Alpha0(), kwargs...)
    return vcat(α0, αg, αψ, θρ)
end

function Theta(m::DrillingRevenue{Unconstrained, TimeTrend};
    θρ=ThetaRho(), αψ=AlphaPsi(), αg=AlphaG(), αT=AlphaT(), α0=Alpha0(), kwargs...)
    return vcat(α0, αg, αψ, αT, θρ)
end

function Theta(m::DrillingRevenue{Unconstrained, TimeTrendDgt0};
    θρ=ThetaRho(), αψ=AlphaPsi(), αg=AlphaG(), αT=AlphaT(), α0=Alpha0(), αD=AlphaD(), kwargs...)
    return vcat(α0, αg, αψ, αT, αD, θρ)
end

function Theta(m::DrillingRevenue{Unconstrained, TimeFE};
    θρ=ThetaRho(), αψ=AlphaPsi(), αg=AlphaG(), αT=AlphaT(), αTFE=AlphaTFE(), α0=Alpha0(), kwargs...)
    _nparm(tech(m)) == length(αTFE) || throw(DimensionMismatch())
    return vcat(α0, αg, αψ, αTFE, θρ)
end

# cost
# ---------

Theta(m::DrillingCost_constant, args...; kwargs...) = vcat(-5.5)

Theta(m::DrillingCost_dgt1, args...; kwargs...) = vcat(-5.5, 1.0)

function Theta(m::DrillingCost_TimeFE, args...; kwargs...)
    c2plus = 1.407
    n = length(ShaleDrillingLikelihood.yearrange(m))
    if n == 5
        timefe = [-9.48651914568503, -6.198466654405615, -4.859543515359247, -4.391926980826075, -4.464122201335934,]
    else
        throw(error("no starting values for length(yearrange($m)) == 5"))
    end
    return vcat(timefe, c2plus)
end

function Theta(m::DrillingCost_TimeFE_rigrate, args...; kwargs...)
    c2plus = 1.407
    crig = -1.0
    n = length(ShaleDrillingLikelihood.yearrange(m))
    if n == 5
        timefe = [-9.48651914568503, -6.198466654405615, -4.859543515359247, -4.391926980826075, -4.464122201335934,]
    else
        throw(error("no starting values for length(yearrange($m)) == 5"))
    end
    return vcat(timefe, c2plus, crig)
end

function Theta(m::DrillingCost_TimeFE_costdiffs, args...; kwargs...)
    n = length(ShaleDrillingLikelihood.yearrange(m))
    if n == 5
        timefe = [-10.548923591552184, -7.034018996118814, -5.758999044446409, -5.539143772743451, -5.579481358288621,]
    else
        throw(error("no starting values for length(yearrange($m)) == 5"))
    end
    return vcat(timefe, 0.017248457193998108, -1.155685021020903, 1.6270414642692912,)
end

function Theta(m::DrillingCost_DoubleTimeFE, args...; kwargs...)
    n = length(ShaleDrillingLikelihood.yearrange(m))
    d1D1 = -0.765358
    if n == 5
        x = [-10.755, -7.973, -6.463, -6.114, -5.986,]
        timefe = vcat(x, x.+d1D1)
    else
        throw(error("no starting values for length(yearrange($m)) == 5"))
    end
    return vcat(timefe, 0.230703, -d1D1+2.01314)
end

# extension
# ---------
Theta(m::ExtensionCost_Zero    , args...; kwargs...) = zeros(0)
Theta(m::ExtensionCost_Constant, args...; kwargs...) = vcat(-1.0)

Theta(m::ScrapValue_Zero    , args...; kwargs...) = zeros(0)
Theta(m::ScrapValue_Constant, args...; kwargs...) = vcat(0.1)
Theta(m::ScrapValue_Price   , args...; kwargs...) = vcat(0.1, 0.1)

CoefLinks(r) = (zeros(Int,0), zeros(Int,0))

CoefLinks(r::DrillReward{<:DrillingRevenue{Unconstrained, NoTrend}}) =
    [(idx_produce_ψ, idx_drill_ψ,), (idx_produce_g, idx_drill_g)]

CoefLinks(r::DrillReward{<:DrillingRevenue{Unconstrained, TimeTrend}}) =
    [
        (idx_produce_ψ, idx_drill_ψ,),
        (idx_produce_g, idx_drill_g),
        (idx_produce_t, idx_drill_t),
    ]

CoefLinks(r::DrillReward{<:DrillingRevenue{Unconstrained, TimeTrendDgt0}}) =
    [
        (idx_produce_ψ, idx_drill_ψ,),
        (idx_produce_g, idx_drill_g),
        (idx_produce_t, idx_drill_t),
        (idx_produce_D, idx_drill_D),
    ]

function CoefLinks(r::DrillReward{<:DrillingRevenue{Unconstrained, TimeFE}})
    idx_d_t = idx_drill_t(r)
    idx_p_t = idx_produce_t() .+ (idx_d_t .- first(idx_d_t))
    cfd = vcat(idx_drill_ψ(r), idx_drill_g(r), idx_d_t)
    cfp = vcat(idx_produce_ψ(), idx_produce_g(), idx_p_t)
    return (cfp, cfd,)
end

# CoefLinks(r::DrillReward) = CoefLinks(revenue(r))
CoefLinks(m::DynamicDrillModel) = CoefLinks(reward(m))
CoefLinks(d::DataDrill) = CoefLinks(_model(d))


text:
struct UniqueValues{S, T1<:AbstractArray, T2<:AbstractArray{S}}
    values::T1
    unique::T2
    value2index::Dict{S, Int}
end

function UniqueValues(col, s = unique(sort(col)))
    value2index = Dict(zip(s, 1:length(s)))
    UniqueValues(col, s, value2index)
end

function (cs::UniqueValues)(s, el)
    scale = to_scale(s)
    @assert scale !== nothing
    @assert typeof(scale) !== typeof(s)
    cs(scale, el)
end

(cs::UniqueValues)(scale::Function, el) = scale(el)

function (cs::UniqueValues)(scale::AbstractArray, el)
    scale[(cs.value2index[el] - 1) % length(scale) + 1]
end

struct Group
    columns::NamedTuple
    f::FunctionOrAnalysis
end

Group(c::NamedTuple) = Group(c, tuple)

Group(v, f::FunctionOrAnalysis = tuple) = Group((color = v,), f)
Group(f::FunctionOrAnalysis = tuple; kwargs...) = Group(values(kwargs), f)

struct ByColumn end
const bycolumn = ByColumn()

Base.isless(::ByColumn, ::ByColumn) = false

combine_gog(f1, f2) = (args...) -> f1(to_tuple(f2(args...))...)
combine_gog(f1, f2::typeof(tuple)) = f1
combine_gog(f1::typeof(tuple), f2) = f2
combine_gog(f1::typeof(tuple), f2::typeof(tuple)) = tuple

Base.merge(g1::Group, g2::Group) = Group(merge(g1.columns, g2.columns), combine_gog(g1.f, g2.f))
Base.merge(f::FunctionOrAnalysis, g::Group) = merge(Group(f), g)
Base.merge(g::Group, f::FunctionOrAnalysis) = merge(g, Group(f))

Base.:*(g1::Group, g2::Group) = merge(g1, g2)

width(v::Tuple) = length(v)
width(v::NamedTuple) = maximum(width, v)
width(v::AbstractVector) = 1
width(v::AbstractArray) = mapreduce(length, *, axes(v)[2:end])

column_length(v::Union{Tuple, NamedTuple}) = column_length(v[1])
column_length(v::AbstractVector) = length(v)
column_length(v::AbstractArray) = length(axes(v)[1])

extract_view(v::Union{Tuple, NamedTuple}, idxs) = map(x -> extract_view(x, idxs), v)
extract_view(v::AbstractVector, idxs) = view(v, idxs)
extract_view(v::AbstractArray, idxs) = view(v, idxs, axes(v)[2:end]...)

extract_view(v::Tuple, idxs, n) = extract_view(v[n], idxs)
extract_view(v::AbstractVector, idxs, n) = view(v, idxs)
function extract_view(v::AbstractArray, idxs, n)
    ax = axes(v)[2:end]
    c = CartesianIndices(ax)[n]
    view(v, idxs, Tuple(c)...)
end

extract_view(v::NamedTuple, idxs, n) = map(t -> extract_view(t, idxs, n), v)

concatenate(v::Union{Tuple, NamedTuple}...) = map(concatenate, v...)
concatenate(v::AbstractArray...) = vcat(v...)


text:
using Pkg
Pkg.activate("../../..")
using InteractionNetworkModels, Distributions, BenchmarkTools
using Distances, StructuredDistances
# The Model(s)

model_mode = Hollywood(-3.0, Poisson(7), 10)
S = sample(model_mode, 10)
S = [[1,1,1,1], [2,2,2,2], [3,3,3,3,3,3]]
V = 1:20

d_lcs = MatchingDistance(FastLCS(100))
d_lcs_t = ThreadedMatchingDistance(LCS())

K_inner = DimensionRange(1, 10)
K_outer = DimensionRange(1, 50)

model = SIM(S, 5.0, d_lcs, V, K_inner, K_outer)
model_t = SIM(S, 5.0, d_lcs_t, V, K_inner, K_outer)

# model_f = SIM(S, 4.0, d_f, V, 50, 50)
mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=5, β=0.6, ν_td=3,
    len_dist=TrGeometric(0.8, K_inner.l, K_inner.u),
    lag=1,
    K=200,
    burn_in=1000
)

@time out=mcmc_sampler(
    model,
    lag=20,
    init=model.mode,
    burn_in=0,
    desired_samples=200
)

@time out=mcmc_sampler(
    model_t,
    lag=20,
    init=model.mode,
    burn_in=0,
    desired_samples=200
)


text:
### In this file we show a step-by-step example for generating a DeBurijnGraph
###  using paired end reads generated by using the Pseudoseq package

using BioSequenceGraphs
using BioSequences
using Pseudoseq

## First we generate data using a reference file in fasta format.
## Below is a small example for illustrative purposes.

Pseudoseq.sequence("ecoli-ref.fasta", "pe-reads.fastq"; ng = 5000, flen = 700, cov = 1, paired = true, rdlen = 20, err = 0.001)

## This generated a file named "pe-reads.fastq" with coverage 2, read length 250 and error rate 1

 r = FASTQ.Reader(open("pe-reads.fastq", "r"))

read_set = Set{LongDNASeq}()

for record in r
    push!(read_set,FASTQ.sequence(record))
end

## generate a set of unique kmers from the set of reads of length 10.
## In our case this generated 742343  kmers
    kmer_set = BioSequenceGraphs.new_extract_canonical_kmers(read_set,10)

## At the end we have generated a dbg with Base.length(BioSequenceGraphs.nodes(dbg))=  742343
<<<<<<< HEAD
=======
start = time(); BioSequenceGraphs.new_deBruijn_Constructor(kmer_set); diff  = time()-start
## Finished  in 605.3275558948517

## First  we test the simple path finder on this dbg
>>>>>>> e837ae91cbfe8ae18f58a8c8082260f9f12361e3

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
module AptMoo

export moo, moomoo, moomoomoo

function moo()
    print(raw"""
                 (__)
                 (oo)
           /------\/
          / |    ||
         *  /\---/\
            ~~   ~~
..."Have you mooed today?"...
""")
end

function moomoo()
    print(raw"""
                 (__)
         _______~(..)~
           ,----\(oo)
          /|____|,'
         * /"\ /\
           ~ ~ ~ ~
..."Have you mooed today?"...
""")
end

function moomoomoo()
    print(raw"""
                     \_/
   m00h  (__)       -(_)-
      \  ~Oo~___     / \
         (..)  |\
___________|_|_|_____________
..."Have you mooed today?"...
""")
end

end # module


text:
push!(LOAD_PATH, "src/")
using ArgParse, nanoJulia, DataFrames, CSV

# Arguments parser
function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table! s begin
		"--inputfile", "-i"
			help = "The path to sequencing summary file"
			arg_type = String
			required = true
		"--outputdir", "-o"
			help = "Folder to output QC results"
			required = true
	end

	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()
	inputfile = parsed_args["inputfile"]
    outputdir = parsed_args["outputdir"]

	if isdir(outputdir)
		println("\033[1;33m[Warning] The output folder has already existed, so it will be renamed with additional suffix \"_bak\"\033[0m")
		new_outputdir = string(outputdir, "_bak")
		mv(outputdir, new_outputdir)
	end
	mkdir(outputdir)

	length2qualityTextFile = joinpath(outputdir, "readlength_vs_readquality.tsv")

    println("\033[1;32m* Start parsing this sequencing summary file...\033[0m")
	output_table = readSeqSummary(inputfile)
	println("\033[1;32m* Writing output file...\033[0m")
	output_table |> CSV.write(length2qualityTextFile)
	totalLength = nanoJulia.totalLen(output_table.length)
	n50 = nanoJulia.calculate_N50(output_table.length, totalLength)
	println("\033[1;32m* Generating summry file and figures...\033[0m")
	generateStatSummary(output_table, totalLength, n50, outputdir)
end

main()


text:
"""
    abstract type AbstractParticle

Subtypes of `AbstractParticle` represent various particles and/or Hamiltonians, such as
particles described the Dirac equation or the non-relativistic Schrödinger equation.

The subtypes are used to determine which physical particle species is meant when dispatching
on generic functions calculating some properties of or related to the particle.
"""
abstract type AbstractParticle end
Base.Broadcast.broadcastable(x::AbstractParticle) = Ref(x)

"""
    struct DiracParticle <: AbstractParticle

Represents a relativistic particle, the free-particle dynamics of which are described by the
Dirac equation.

```math
(\\gamma^\\mu \\partial_\\mu - m) \\Psi = 0
```

Note that we set the zero of the energy to ``+mc^2``.
"""
struct DiracParticle <: AbstractParticle
    mass :: Float64
    function DiracParticle(mass)
        mass > 0 || throw(DomainError(mass, "mass must be positive"))
        new(mass)
    end
end
"""
    const DiracElectron :: DiracParticle

A relativistic Dirac fermion ([`DiracParticle`](@ref)) with the mass of an electron.
"""
const DiracElectron = DiracParticle(1.0)
"""
    const DiracMuon :: DiracParticle

A relativistic Dirac fermion ([`DiracParticle`](@ref)) with the mass of a muon.
"""
const DiracMuon = DiracParticle(muon_mass_au)

"""
    struct NRParticle <: AbstractParticle

Represents a non-relativistic particle, described by the non-relativistic Schrödinger
equation.

```math
i \\frac{\\partial\\Psi}{\\partial t} = - \\frac{\\nabla^2}{2 m^2} \\Psi
```
"""
struct NRParticle <: AbstractParticle
    mass :: Float64
    function NRParticle(mass)
        mass > 0 || throw(DomainError(mass, "mass must be positive"))
        new(mass)
    end
end
"""
    const NRElectron

A non-relativistic particle ([`NRParticle`](@ref)) with the mass of an electron.
"""
const NRElectron = NRParticle(1.0)
"""
    const NRMuon

A non-relativistic particle ([`NRParticle`](@ref)) with the mass of a muon.
"""
const NRMuon = NRParticle(muon_mass_au)

"""
    hydrogenic_energy(
        dynamics::AbstractParticle[, nucleus::AbstractNuclearModel], Z;
        qnumbers...
    ) -> Float64

Returns the energy (in atomic units) of a hydrogenic bound state for given dynamics
(particle), central potential generated by a nuclear charge distribution, and quantum numbers.

Each nuclear model is, by the definition of [`AbstractNuclearModel`](@ref), normalized to unity.
This means that the nuclear charge ``Z`` needs to be passed as a separate argument.
Also, the second argument (`nucleus`) can be omitted and it defaults to [`PointNucleus`](@ref).

## Quantum numbers

The keyword arguments (`qnumbers...`) can be used to specify all the necessarily quantum
numbers to uniquely identify the state. In some cases, some quantum numbers may be omitted.

In the non-relativistic case we need to specify `n` and `ℓ`, whereas in the relativistic case
it is `n` and `κ`.
"""
function hydrogenic_energy end
# We let the second argument default to `PointNucleus`:
hydrogenic_energy(p::AbstractParticle, Z::Number; kwargs...) = hydrogenic_energy(p, PointNucleus(), Z; kwargs...)

function hydrogenic_energy(p::NRParticle, ::PointNucleus, Z::Real; n::Integer, ℓ = 0, m = 0, ms = 0)
    n >= 1 || throw(DomainError(n, "n must be a positive integer"))
    - p.mass * Z^2 / (2 * n^2)
end

function hydrogenic_energy(p::DiracParticle, ::PointNucleus, Z::Number; n::Integer, κ::Integer, mj = 0)
    n >= 1 || throw(DomainError(n, "n must be a positive integer"))
    if Z isa Real
        Z < 0 && throw(DomainError(Z, "0 ≤ Z ≤ α⁻¹ for real Z"))
        Z > 1/α && throw(DomainError(Z, "0 ≤ Z ≤ α⁻¹ for real Z, pass as complex(Z) to evaluate at Z > α⁻¹"))
    end
    # Note: using this expression is not as numerically precise as it could probably be, in
    # terms of floating point precision. It leads to an error of about 1e-12.
    (1 / sqrt(1 + (Z*α / (n - abs(κ) + sqrt(κ^2 - (Z*α)^2)))^2) - 1) * p.mass * α^-2
end

function hydrogenic_energy(p::NRParticle, nm::UniformShellNucleus, Z::Real; n::Integer, ℓ::Integer, m = 0, ms = 0, verbose=false)
    n >= 1 || throw(DomainError(n, "n must be a positive integer"))
    f = E -> f_andrae_nonrel_topslice(; Z=Z, R=nm.R, r=nm.R, ℓ=ℓ, E = E)
    # Pick the initial energy range to be between the point energies for n and n+1:
    emin, emax = hydrogenic_energy(p, Z; n = n), hydrogenic_energy(p, Z; n = n + 1)
    # In the expressions, we need to evaluate sqrt(2*(E+Z/R)), so E ≥ -Z/R for f to even
    # make sense. So we must constrain the bounding box when Z becomes large enough.
    emin = max(emin, -Z/nm.R)
    find_zero(f, (emin, emax), verbose = verbose)
end

# See docs/src/implementation/hydrogenic.jl for documentation.
function f_andrae_nonrel_topslice(; Z, R, ℓ, r, E)
    # Inner variables
    α = sqrt(2 * (E + Z/R))
    xin = α*r
    J1 = sf_bessel_jl(ℓ + 1, xin)
    J2 = sf_bessel_jl(ℓ, xin)
    # Outer variables
    βnℓ = sqrt(-2*E)
    νnℓ = Z/βnℓ
    b = 2*(ℓ+1)
    a = b/2 - νnℓ
    xout = 2*βnℓ*r
    U1 = try
        sf_hyperg_U(a+1, b+1, xout)
    catch e
        @error "sf_hyperg_U(a+1, b+1, xout)" a b xout r
        rethrow(e)
    end
    U2 = try
        sf_hyperg_U(a, b, xout)
    catch e
        @error "sf_hyperg_U(a, b, xout)" a b xout r
        rethrow(e)
    end
    # Multiply both sides and subtract
    βnℓ*J2*U2 + 2*βnℓ*a*J2*U1 - α*J1*U2
end

function hydrogenic_energy(p::DiracParticle, nm::UniformShellNucleus, Z::Real; n::Integer, κ::Integer, mj = 0, verbose=false)
    n >= 1 || throw(DomainError(n, "n must be a positive integer"))
    # We'll do root finding in a bracketing box, so we first need to decide the bracketing
    # interval. As we know that the energies always slightly shift up, we choose the
    # interval to be between the n and n+1 point energies. In most cases, this will
    # guarantee that we will have exactly one zero within the interval. Where this will fail
    # is when the FNC shift is so large that energy is greater than the next PNC energy.
    emin = hydrogenic_energy(p, Z; κ = κ, n = n)
    emax = hydrogenic_energy(p, Z; κ = κ, n = n + 1)
    # We shift the bracketing interval a bit because of the numerical errors in
    # hydrogenic_energy. With 64-bit floats, it's accuracy is only about 1e-12, so we are
    # conservative here and shift everything by 1e-10. This is still many orders of magnitude
    # smaller than the energy difference between any two consective n values (for low enough
    # n values), so this is pretty safe thing to do.
    #
    # However, the result is that we are guaranteed to have a sane bracketing interval,
    # whereas without this shift we sometimes place the start or end of the interval _after_
    # the corresponding zero, as opposed to having it before it.
    emin -= 1e-10
    emax -= 1e-10
    # In the expressions, we need to evaluate α_nκ, which require E ≥ -Z/R for f to not be
    # complex. So we constrain the bounding box when Z becomes large enough.
    emin = max(emin, -Z/nm.R)
    f = E -> f_andrae_dirac_topslice(; Z = Z, R = nm.R, κ = κ, E = E)
    find_zero(f, (emin, emax), verbose = verbose)
end

# See docs/src/implementation/hydrogenic.jl for documentation.
function f_andrae_dirac_topslice(; Z, R, κ, E)
    # Note: E = 0 is defined to be mc², so the energy will always be negative for bound
    # states
    c = α^-1 # since we are in atomic units
    ℓ = κ2ℓ(κ)
    # Following are from equation (282)
    #
    # Wnκ = E + c^2
    # C⁺, C⁻ = sqrt(c^2 + Wnκ), sqrt(c^2 - Wnκ)
    # βnκ = C⁺ * C⁻ / c
    #
    # However, as the expressions in the paper are subject to numerical roundoff errors that
    # limit the resolution of E to about 1e-12 (for 64-bit floats), we slightly rearrange
    # the equations for βnκ and also νnκ below.
    Ec⁻² = E/c^2
    βnκ = sqrt(-E * (2 + Ec⁻²))
    x_out = 2*βnκ*R # x at r = R for outer part
    A = Z/βnκ + κ
    γ = sqrt(κ^2 - (Z/c)^2)
    b = 2γ +  1
    # We also modify νnκ
    #
    # νnκ = 0.5 + Z*Wnκ/(βnκ * c^2)
    #
    # to reduce roundoff errors.
    νnκ = 0.5 + (Z/βnκ) * (Ec⁻² + 1)
    a = b/2 - νnκ
    # From 283
    U00 = sf_hyperg_U(a, b, x_out)
    U10 = sf_hyperg_U(a + 1, b, x_out)
    U11 = sf_hyperg_U(a + 1, b + 1, x_out)
    U21 = sf_hyperg_U(a + 2, b + 1, x_out)
    Ud = U00 + A*U10 # denomimator in (283)
    # From 273 & 274
    αnκ = sqrt(2*(E + Z/R) + ((E + Z/R)/c)^2)
    x_in = αnκ*R # x at r = R for inner part
    J0 = sf_bessel_jl(ℓ, x_in)
    J1 = sf_bessel_jl(ℓ + 1, x_in)
    # Combining all of them now:
    ((ℓ + 1 - γ)/R + βnκ)*Ud*J0 - αnκ*J1*Ud + 2*βnκ*(a*U11 + A*(a + 1)*U21)*J0
end

# From AtomicLevels.jl (MIT license)
function κ2ℓ(κ::Integer)
    κ == zero(κ) && throw(ArgumentError("κ can not be zero"))
    (κ < 0) ? -(κ + 1) : κ
end


text:
# TODO:
# - type tags?
# - eltype of Ptr and Array is only `first Element`
# - recursive types

import LLVM: refcheck
LLVM.@checked struct TypeTree
    ref::API.CTypeTreeRef
end
Base.unsafe_convert(::Type{API.CTypeTreeRef}, tt::TypeTree) = tt.ref
LLVM.dispose(tt::TypeTree) = API.EnzymeFreeTypeTree(tt)

TypeTree() = TypeTree(API.EnzymeNewTypeTree())
TypeTree(CT, ctx) = TypeTree(API.EnzymeNewTypeTreeCT(CT, ctx))
function TypeTree(CT, idx, ctx)
    tt = TypeTree(CT, ctx)
    only!(tt, idx)
    return tt
end
Base.copy(tt::TypeTree) = TypeTree(API.EnzymeNewTypeTreeTR(tt))
Base.copy!(dst::TypeTree, src::TypeTree) = API.EnzymeSetTypeTree(dst, src)

# function Base.show(io::IO, ::MIME"text/plain", tt::TypeTree)
#     print(io, "TypeTree: ")
#     for data in tt.data
#         println(io)
#         print(io, "- ")
#         show(io, data)
#     end
# end

function only!(tt::TypeTree, offset::Integer)
    API.EnzymeTypeTreeOnlyEq(tt, offset)
end

function shift!(tt::TypeTree, dl, offset, maxSize, addOffset)
    API.EnzymeTypeTreeShiftIndiciesEq(tt, dl, offset, maxSize, addOffset)
end

function merge!(dst::TypeTree, src::TypeTree; consume=true)
    API.EnzymeMergeTypeTree(dst, src)
    LLVM.dispose(src)
    return nothing
end

function typetree(::Type{T}, ctx, dl) where T <: Integer
    tt = TypeTree()
    for i in 1:sizeof(T)
        merge!(tt, TypeTree(API.DT_Integer, i-1, ctx))
    end
    return tt
end

function typetree(::Type{Float16}, ctx, dl)
    return TypeTree(API.DT_Half, 0, ctx)
end

function typetree(::Type{Float32}, ctx, dl)
    return TypeTree(API.DT_Float, 0, ctx)
end

function typetree(::Type{Float64}, ctx, dl)
    return TypeTree(API.DT_Double, 0, ctx)
end

function typetree(::Type{<:Ptr{T}}, ctx, dl) where T
    tt = typetree(T, ctx, dl)
    merge!(tt, TypeTree(API.DT_Pointer, ctx))
    only!(tt, 0)
    return tt
end

function typetree(::Type{<:Array{T}}, ctx, dl) where T
    offset = 0

    tt = typetree(T, ctx, dl)
    merge!(tt, TypeTree(API.DT_Pointer, ctx))
    only!(tt, offset)

    offset += sizeof(Ptr{Cvoid})

    sizeofstruct = offset + 2 + 2 + 4 + 2 * sizeof(Csize_t)
    if true # STORE_ARRAY_LEN
        sizeofstruct += sizeof(Csize_t)
    end

    for i in offset:sizeofstruct
        merge!(tt, TypeTree(API.DT_Integer, i, ctx))
    end
    return tt
end

function typetree(@nospecialize(T), ctx, dl)
    if fieldcount(T) == 0
        error("$T is unknown leaf")
    end

    tt = TypeTree()
    for f in 1:fieldcount(T)
        offset  = fieldoffset(T, f)
        subT    = fieldtype(T, f)
        subtree = typetree(subT, ctx, dl)

        # Allocated inline so adjust first path
        if subT.isinlinealloc
            shift!(subtree, dl, 0, sizeof(subT), offset)
        else
            merge!(subtree, TypeTree(API.DT_Pointer, ctx))
            only!(subtree, offset)
        end

        merge!(tt, subtree)
    end
    return tt
end

struct FnTypeInfo
    rTT::TypeTree
    argTTs::Vector{TypeTree}
    known_values::Vector{API.IntList}
end
Base.cconvert(::Type{API.CFnTypeInfo}, fnti::FnTypeInfo) = fnti
function Base.unsafe_convert(::Type{API.CFnTypeInfo}, fnti::FnTypeInfo)
    args_kv = Base.unsafe_convert(Ptr{API.IntList}, Base.cconvert(Ptr{API.IntList}, fnti.known_values))
    rTT = Base.unsafe_convert(API.CTypeTreeRef, Base.cconvert(API.CTypeTreeRef, fnti.rTT))

    tts = API.CTypeTreeRef[]
    for tt in fnti.argTTs
        raw_tt = Base.unsafe_convert(API.CTypeTreeRef, Base.cconvert(API.CTypeTreeRef, tt))
        push!(tts, raw_tt)
    end
    argTTs = Base.unsafe_convert(Ptr{API.CTypeTreeRef}, Base.cconvert(Ptr{API.CTypeTreeRef}, tts))
    return API.CFnTypeInfo(argTTs, rTT, args_kv)
end


text:
"""
"""
mesh

function mesh(bc, st, n, args...)
    coor = _mesh.(bc, st, n, args...)
    n̅ = @. n + 2width(st) * has(bc)
    eye = Ones{Bool}.(n̅)
    x̅, x̅ₘ, dx̅ = _mesh(coor, eye)
    n̅, x̅, x̅ₘ, dx̅
end

function _mesh(bc::Dirichlet, st, n::Int, func=identity)
    w = width(st)
    func.(range(-w, n+w-1, n+2w) ./ (n-1)),
    func.(range(1-2w, 2(n+w)-1, n+2w) ./ (2n-2)),
    func.(range(1-w, n+w, n+2w) ./ (n-1)) .-
    func.(range(-w, n+w-1, n+2w) ./ (n-1))
end

_mesh(::Periodic, ::Any, n::Int, func=identity) =
    func.(range(0, n-1, n) ./ n),
    func.(range(1, 2n-1, n) ./ 2n),
    func.(range(1, n, n) ./ n) .-
    func.(range(0, n-1, n) ./ n)

_mesh(coor::NTuple{1}, eye::NTuple{1}) =
    (coor[1][1],),
    (coor[1][2],),
    (coor[1][3],)

_mesh(coor::NTuple{2}, eye::NTuple{2}) =
    (ApplyArray(kron, eye[2], coor[1][1]),
     ApplyArray(kron, coor[2][1], eye[1])),
    (ApplyArray(kron, eye[2], coor[1][2]),
     ApplyArray(kron, coor[2][2], eye[1])),
    (ApplyArray(kron, eye[2], coor[1][3]),
     ApplyArray(kron, coor[2][3], eye[1]))

_mesh(coor::NTuple{3}, eye::NTuple{3}) =
    (ApplyArray(kron, eye[3], eye[2], coor[1][1]),
     ApplyArray(kron, eye[3], coor[2][1], eye[1]),
     ApplyArray(kron, coor[3][1], eye[2], eye[1])),
    (ApplyArray(kron, eye[3], eye[2], coor[1][2]),
     ApplyArray(kron, eye[3], coor[2][2], eye[1]),
     ApplyArray(kron, coor[3][2], eye[2], eye[1])),
    (ApplyArray(kron, eye[3], eye[2], coor[1][3]),
     ApplyArray(kron, eye[3], coor[2][3], eye[1]),
     ApplyArray(kron, coor[3][3], eye[2], eye[1]))