text:
# interface for PetscMat that also maps cleanly onto the Julia Matrix interface
# the Julia sparse matrix interface is ill-defined, so this may take some
# liberties

#=
function PetscMat(mglobal::Integer, nglobal::Integer, comm::MPI_Comm; mlocal=PETSC_DECIDE, nlocal=PETSC_DECIDE)

  return mat
end
=#

"""
  Constructor

  **Inputs**

   * mglobal: first dimension global size (or PETSC_DECIDE)
   * nglobal: second dimension global size (or PETSC_DECIDE)
   * format: matrix format
   * comm: MPI communicator

  **Keyword Arguments**

   * mlocal: first dimension local size (or PETSC_DECIDE)
   * nlocal: second dimension local size (or PETSC_DECIDE)
"""
function PetscMat(mglobal::Integer, nglobal::Integer, format, comm::MPI_Comm;
                  mlocal=PETSC_DECIDE, nlocal=PETSC_DECIDE)

  mat = PetscMat(comm)
  MatSetSizes(mat, mlocal, nlocal, mglobal, nglobal)
  MatSetType(mat, format)

  return mat
end

#=
function PetscDestroy(vec::PetscMat)
  if (vec.pobj != 0)
    err = ccall( (:MatDestroy,  libpetsclocation), PetscErrorCode, (Ptr{Ptr{Void}},), &vec.pobj);
  end
  vec.pobj = 0

#    println("Petsc Mat Destroy called")
#    sleep(5)
  return 0
end
=#

"""
  PetscDestroy for AbstractMatrix.  No-op
"""
function PetscDestroy(mat::AbstractMatrix)
end

# 1-based indexing for both Petsc and regular matrices
#----------------------------------------------------------------------------
"""
  1-based indexing for both regular and Pets matrices.
  Note that Petsc treats arrays at being row-major, so it is recommened
  to set MatSetOption(mat, MAT_ROW_ORIENTED, PETSC_FALSE) before using
  this function.

  **Inputs**

   * flag: INSERT_VALUES or ADD_VALUES.  Note that the first one
           result in non-deterministic behavior in parallel (in the general
           case)

  * vals: the values, must be length(idxm) x length(idxn)

  **Inputs/Outputs**

   * mat: a matrix, can be a Petsc matrix or a julia matrix
   * idxm: the row numbers
   * idxn: the column numbers

  Note that idxm and idxn are listed as input/outputs because they may be
  modified by this function, however when the function returns they
  will have the same values as on entry.  This is needed to accomodate the
  fact that Petsc uses 1 based indexing internally.

  This function is optimized for PetscMat and SparseMatrixCSC

  Aliasing restriction: idxm and idxn cannot alias
"""
function set_values1!(mat::PetscMat, idxm::Array{PetscInt}, idxn::Array{PetscInt},
                      vals::Array{PetscScalar}, flag::Integer=INSERT_VALUES)

  for i=1:length(idxm)
    idxm[i] -= 1
  end

  for i=1:length(idxn)
    idxn[i] -= 1
  end

  err = MatSetValues(mat, idxm, idxn, vals, flag)

  for i=1:length(idxm)
    idxm[i] += 1
  end

  for i=1:length(idxn)
    idxn[i] += 1
  end

  return err
end

"""
  Method for AbstractMatrix
"""
function set_values1!{T}(mat::AbstractMatrix, idxm::Array{PetscInt},
                         idxn::Array{PetscInt}, vals::Array{T},
                         flag::Integer=INSERT_VALUES)

  if flag == INSERT_VALUES
    for i=1:length(idxn)
      for j=1:length(idxm)
        mat[idxm[j], idxn[i]] = vals[j, i]
      end
    end
  elseif flag == ADD_VALUES
    for i=1:length(idxn)
      for j=1:length(idxm)
        mat[idxm[j], idxn[i]] += vals[j, i]
      end
    end
  end

  return PetscErrorCode(0)
end

"""
  Method for `SparseMatrixCSC`
"""
function set_values1!{T}(mat::SparseMatrixCSC, idxm::Array{PetscInt},
                         idxn::Array{PetscInt}, vals::Array{T},
                         flag::Integer=INSERT_VALUES)

  if flag == INSERT_VALUES
    for i=1:length(idxn)
      for j=1:length(idxm)
        mat[idxm[j], idxn[i]] = vals[j, i]
      end
    end
  elseif flag == ADD_VALUES  # optimized += implementation
    # hoist
    colptr = mat.colptr
    rowval = mat.rowval
    nzval = mat.nzval
    for i=1:length(idxn)
      row_start = colptr[idxn[i]]
      row_end = colptr[idxn[i]+1] - 1
      rowvals_extract = unsafe_aview(rowval, row_start:row_end)
      for j=1:length(idxm)
        idx = searchsortedfirst(rowvals_extract, idxm[j])
        idx = row_start + idx - 1
        nzval[idx] += vals[j, i]
      end
    end
  end

  return PetscErrorCode(0)
end

"""
  Like [`set_values1!`](@ref), but retrieves values.  See that function for
  the meanings of the arguments. Note that Petsc does
  not support getting values for the non-local block of the matrix

  **Inputs**

   * mat: a matrix, can be a Petsc matrix or a julia matrix

  **Inputs/Outputs**

   * idxm
   * idxn
   * vals

  Aliasing restrictions: idxm and idxn cannot alias
"""
function get_values1!(mat::PetscMat, idxm::Array{PetscInt},
                      idxn::Array{PetscInt},
                      vals::Array{PetscScalar})

  for i=1:length(idxm)
    idxm[i] -= 1
  end

  for i=1:length(idxn)
    idxn[i] -= 1
  end

  err = MatGetValues(mat, idxm, idxn, vals)

  for i=1:length(idxm)
    idxm[i] += 1
  end

  for i=1:length(idxn)
    idxn[i] += 1
  end

  return err
end

"""
  Method for `AbstractMatrix`
"""
function get_values1!(mat::AbstractMatrix, idxm::Array{PetscInt},
                      idxn::Array{PetscInt}, vals::Array)

  for i=1:length(idxn)
    for j=1:length(idxm)
      vals[j, i] = mat[idxm[j], idxn[i]]
    end
  end

  return PetscErrorCode(0)
end

"""
  Begin matrix assembly for `PetscMat`.  No-op for Julia matrices

  **Inputs**

   * mat: AbstractMatrix
   * flag: type of matrix assembly (see [`MatAssemblyBegin`](@ref)
"""
function assembly_begin(mat::PetscMat, flg::Integer)
  MatAssemblyBegin(mat, flg)
end

function assembly_begin(mat::AbstractMatrix, flg::Integer)
end

"""
  Counterpart of [`assembly_end`](@ref)

  **Inputs**

   * mat: AbstractMatrix
   * flg: type of matrix assembly
"""
function assembly_end(mat::PetscMat, flg::Integer)
  MatAssemblyEnd(mat, flg)
end

function assembly_end(mat::AbstractMatrix, flg::Integer)
end

#TODO: transpose, inplace and/out out of place
import Base: transpose, transpose!

# PetscView?
"""
  Print a non-Petsc matrix to a given IO (a Julia IO, not a Petsc IO)

  **Inputs**

   * A: the matrix
   * f: an IO, default STDOUT
"""
function PetscView(A::AbstractMatrix, f::IO=STDOUT)
  println(f, "A = \n", A)
end

"""
  Size of local part of matrix

  **Inputs**

   * A: AbstractMatrix
"""
function size_local(A::PetscMat)
  return MatGetLocalSize(A)
end

function size_local(A::AbstractMatrix)
  return size(A)
end

"""
  Global size of matrix, same as size_local() for serial matrices

  **Inputs**

   * A: AbstractMatrix
"""
function size_global(A::PetscMat)
  return MatGetSize(A)
end

function size_global(A::AbstractMatrix)
  return size(A)
end

"""
  Returns the rows owned by this process (1-based)

  **Inputs**

   * A: a matrix

  **Outputs**

   * rng: a UnitRange
"""
function local_indices(A::PetscMat)
  low, high = MatGetOwnershipRange(A)
  return (low+1):high
end

function local_indices(A::AbstractMatrix)
  return 1:size(A, 1)
end

"""
  Fill the matrix with zeros. The sparsity pattern of the matrix (if applicable)
  should be defined before this function is called

  **Inputs**

   * A: AbstractMatrix
"""
function fill_zero!(A::PetscMat)
  MatZeroEntries(A)
end

function fill_zero!(A::AbstractMatrix)
  fill!(A, 0.0)
end

function fill_zero!(A::SparseMatrixCSC)
  fill!(A.nzval, 0.0)
end

import Base.scale!
"""
  scale! for Petsc matrix
"""
function scale!(A::PetscMat, a::Number)
  MatScale(A, PetscScalar(a))
end

"""
  Adds the specified value to the diagonal of the matrix

  **Inputs**

   * A: the matrix
   * a: the value
"""
function diagonal_shift!(A::PetscMat, a::Number)
  MatShift(A, PetscScalar(a))
end

function diagonal_shift!(A::AbstractMatrix, a::Number)

  nmax = min(size(A, 1), size(A, 2))
  @simd for i=1:nmax
    A[i, i] += a
  end
end
#TODO; optimized implementation for SparseMatrixCSC

import Base: norm, vecnorm

"""
  Norm for Petsc matrices, 1, 2, and infinity norms supported
"""
function norm(A::PetscMat, p::Real)
  if p == 1
    _p = NORM_1
  elseif p == 2
    _p = NORM_2
  elseif p == Inf
    _p = NORM_INFINITY
  end
  MatNorm(A, _p)
end

"""
  Frobenius norm, consistent with Julias interface
"""
function vecnorm(A::PetscMat)
  MatNorm(A, NORM_FROBENIUS)
end

import Base: A_mul_B!, At_mul_B!

"""
  Computes b = A*x
"""
function A_mul_B!(b::PetscVec, A::PetscMat, x::PetscVec)
  MatMult(A, x, b)
  return b
end

"""
  Computes b = A.'*x
"""
function At_mul_B!(b::PetscVec, A::PetscMat, x::PetscVec)
  MatMultTranspose(A, x, b)
  return b
end


text:
using Minesweeper

#import GUI functions
import_gui()

# Creates game in GUI.
# Default game is: Game(dims=(10,10), n_mines=15)
# Pass game object for custom game
# Start new game: File->New Game
# Game setup: File->Setup
# Flagging: switch on and off using toggle slider in top left corner
generate_gui()

# configure game in File-> Setup or pass custom game
game = Game(dims=(10,20), n_mines=35)
generate_gui(game)

ext: jl


text: println("Hello AD from Binder!")

max_stars_repo_name: amaurydeldicque/my-first-binder

max_issues_repo_name: amaurydeldicque/my-first-binder

max_forks_repo_name: amaurydeldicque/my-first-binder


text:
using Match
import Base.Meta.isexpr

# argument such as int[2] were translated to NTuple{UInt8, 16}
# but it seems they are never referenced so its fine
const as_ntuple = [
    :avcodec_align_dimensions2
    :av_display_rotation_get
    :av_display_rotation_set
    :av_display_matrix_flip
    :av_image_fill_max_pixsteps
    :av_image_fill_linesizes
    :av_image_fill_pointers
    :av_image_alloc
    :av_image_copy
    :av_image_copy_uc_from
    :av_image_fill_black
    :av_read_image_line2
    :av_read_image_line
    :av_write_image_line2
    :av_write_image_line
    :av_tree_find
    :av_xtea_init
    :av_xtea_le_init
    :avcodec_align_dimensions2
]

function rewrite(e::Expr)
    @match e begin
        Expr(:function, [fncall, body])  =>  rewrite_fn(e, fncall, body)
        Expr(:const, [arg])              =>  rewrite_const(e)
        _                                =>  e
    end
end

function rewrite_const(e)
    body = e.args[]
    lhs, rhs = body.args
    name = string(lhs)
    if match(r"^AV_PIX_FMT_[\w\d]+$", name) |> !isnothing
        func, be, le = rhs.args
        func == :AV_PIX_FMT_NE || return e
        rhs = Expr(:macrocall, Symbol("@AV_PIX_FMT_NE"), nothing, be, le)
        body.args[2] = rhs
    end
    e
end

# Rewrite function signatures to be more friendly
function rewrite_fn(e, fncall, body, use_strpack=false)
    fncall.args[1] isa Symbol || return e

    parms = Any[]
    content = Any[]

    push!(parms, fncall.args[1])  # function name

    # Add explicit conversions for Ptr, UInt32/Int32, and va_list
    for call_arg in fncall.args[2:end]
        @match call_arg begin
            # Don't type Ptr{x} types
            Expr(:(::), [sym, Expr(:curly, [:Ptr, _])]) => push!(parms, sym)

            # Type all integers as Integer
            Expr(:(::), [sym, (:UInt32 || :Cuint || :Int32 || :Cint)]) => (sym; push!(parms, :($sym::Integer)))

            # va_list
            Expr(:(::), [sym, :va_list]) => push!(parms, sym)

            # Everything else is unchanged
            _ => push!(parms, call_arg)
        end
    end

    if !isexpr(body, :block)
        # body just consists of the ccall
        push!(content, body)
    else
        append!(content, body.args)
    end

    fncall = Expr(:call, parms...)
    body = Expr(:block, content...)

    new = Expr(e.head, fncall, body)

    return new
end


text:
module HierarchicalMatrices

using Base

include("BBox.jl")
include("BinTree.jl")
include("ABtDecomposition.jl")
include("HMatrix.jl")

end # module

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices


text:
using GNSSTools
using Statistics
using PyPlot
pygui(true)

t = Array(0:0.1:100)
carrier = sin.((2π/10).*t)

# I channel
ranging_code_I = l5i_codes[1][1:10]
ranging_code_I[ranging_code_I .== 0] .= -1
neuman_code_I = nh10[5:9]
neuman_code_I[neuman_code_I .== 0] .= -1
navbits_I = [1, 1, -1]

# Q channel
ranging_code_Q = l5q_codes[1][1:10]
ranging_code_Q[ranging_code_Q .== 0] .= -1
neuman_code_Q = nh20[4:7]
neuman_code_Q[neuman_code_Q .== 0] .= -1
navbits_Q = [-1, -1, 1]

t_length = 1e-3
signal = definesignal(define_l1ca_code_type(), 5e6, 2*t_length; f_if=1.25e6,
                      include_thermal_noise=false, include_phase_noise=false,
                      include_adc=false, include_carrier=true)
generatesignal!(signal)
final_signal = real.(deepcopy(signal.data[1:200]))
final_signal = final_signal ./ maximum(final_signal)

signal = definesignal(define_l1ca_code_type(), 5e6, 2*t_length; f_if=1.25e6,
                      include_thermal_noise=true, include_phase_noise=true,
                      include_adc=false, include_carrier=true)
generatesignal!(signal)
final_signal_noise = real.(deepcopy(signal.data[1:200]))
final_signal_noise = final_signal_noise ./ maximum(final_signal_noise)

ranging_code_multiple = [ranging_code_I; ranging_code_I]

ylims = [-1.50, 1.50]
fig = figure(figsize=(6,3))
ax1 = fig.add_subplot(4,2,1)
ax1.step(ranging_code_I, "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax2 = fig.add_subplot(4,2,2)
ax2.step(ranging_code_multiple, "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax3 = fig.add_subplot(4,2,3)
ax3.plot(final_signal[1:100], "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax4 = fig.add_subplot(4,2,4)
ax4.plot(final_signal[101:200], "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax5 = fig.add_subplot(4,2,5)
ax5.plot(final_signal, "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax6 = fig.add_subplot(4,2,6)
ax6.plot(final_signal_noise[1:100], "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax7 = fig.add_subplot(4,2,7)
ax7.plot(final_signal_noise[101:200], "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax8 = fig.add_subplot(4,2,8)
ax8.plot(final_signal_noise, "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
# ax4 = fig.add_subplot(4,2,2)
# ax4.step(ranging_code_Q, "k")
# ylim(ylims)
# tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
# ax5 = fig.add_subplot(4,2,4)
# ax5.step(neuman_code_Q, "k")
# ylim(ylims)
# tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
# subplots_adjust(hspace=0.3, wspace=0.4, top=0.91, left=0.08, right=0.93, bottom=0.1)
savefig("figures/ch3_figure_essentials_signals.svg", dpi=300)

n0 = 3125
n0_offset = -25
n0_range = 100
signal = definesignal(define_l1ca_code_type(), 5e6, 1; f_if=1.25e6,
                      include_thermal_noise=true, include_phase_noise=true,
                      include_adc=false, include_carrier=true,
                      code_start_idx=n0, CN0=40)
generatesignal!(signal)
replica = definesignal(define_l1ca_code_type(), 5e6, 1e-3)
fd_est, n0_est, SNR_est, corr_result1 = courseacquisition(signal, replica, 1;
                                                         return_corrresult=true,
                                                         M=1,
                                                         start_idx=1)
fd_est, n0_est, SNR_est, corr_result2 = courseacquisition(signal, replica, 1;
                                                         return_corrresult=true,
                                                         M=1,
                                                         start_idx=1+5000)
fd_est, n0_est, SNR_est, corr_result3 = courseacquisition(signal, replica, 1;
                                                         return_corrresult=true,
                                                         M=10,
                                                         start_idx=1)

corr_crop1 = corr_result1[:,n0-n0_range+n0_offset:n0+n0_range+n0_offset]
corr_crop2 = corr_result2[:,n0-n0_range+n0_offset:n0+n0_range+n0_offset]
corr_crop3 = corr_result3[:,n0-n0_range+n0_offset:n0+n0_range+n0_offset]

mplot3d = PyPlot.PyObject(PyPlot.axes3D)
fig = figure(figsize=(4,5))
ax1 = fig.add_subplot(3, 1, 1, projection="3d")
x, y = meshgrid(Array(1:size(corr_crop1)[2]), Array(1:size(corr_crop1)[1]))
surf(x, y, corr_crop1, antialiased=true, shade=true, alpha=0.9, cmap="coolwarm")
# ax1.grid(false)
# ax1.set_xticks([])
# ax1.set_yticks([])
# ax1.set_zticks([])
axis("off")
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)

ax3 = fig.add_subplot(3, 1, 2, projection="3d")
x, y = meshgrid(Array(1:size(corr_crop2)[2]), Array(1:size(corr_crop2)[1]))
surf(x, y, corr_crop2, antialiased=true, shade=true, alpha=0.9, cmap="coolwarm")
axis("off")
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)

ax3 = fig.add_subplot(3, 1, 3, projection="3d")
x, y = meshgrid(Array(1:size(corr_crop3)[2]), Array(1:size(corr_crop3)[1]))
surf(x, y, corr_crop3, antialiased=true, shade=true, alpha=0.9, cmap="coolwarm")
axis("off")
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)

savefig("figures/ch3_figure_essentials_acquisition_noncoherent.svg", dpi=300)

fd_est, n0_est, SNR_est, corr_result4 = courseacquisition(signal, replica, 1;
                                                         return_corrresult=true,
                                                         M=2,
                                                         start_idx=1)
replica2 = definesignal(define_l1ca_code_type(), 5e6, 2e-3)
fd_est, n0_est, SNR_est, corr_result5 = courseacquisition(signal, replica2, 1;
                                                         return_corrresult=true,
                                                         M=1,
                                                         start_idx=1)

fd_est, n0_est, SNR_est, corr_result6 = courseacquisition(signal, replica, 1;
                                                         return_corrresult=true,
                                                         M=10,
                                                         start_idx=1)
replica2 = definesignal(define_l1ca_code_type(), 5e6, 10e-3)
fd_est, n0_est, SNR_est, corr_result7 = courseacquisition(signal, replica2, 1;
                                                         return_corrresult=true,
                                                         M=1,
                                                         start_idx=1)

ylims = [-0.50, 1.50]
fig = figure(figsize=(6,3))
ax1 = fig.add_subplot(4,2,1)
result = corr_result1[argmax(corr_result1)[1],:]
ax1.plot(result ./ maximum(result), "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax2 = fig.add_subplot(4,2,2)
result = corr_result2[argmax(corr_result2)[1],:]
ax2.plot(result ./ maximum(result), "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax3 = fig.add_subplot(4,2,3)
result = corr_result5[argmax(corr_result5)[1],1:5000]
ax3.plot(result ./ maximum(result), "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax4 = fig.add_subplot(4,2,4)
result = corr_result4[argmax(corr_result4)[1],:]
ax4.plot(result ./ maximum(result), "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax5 = fig.add_subplot(4,2,5)
result = corr_result7[argmax(corr_result7)[1],1:5000]
ax5.plot(result ./ maximum(result), "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
ax6 = fig.add_subplot(4,2,6)
result = corr_result6[argmax(corr_result6)[1],:]
ax6.plot(result ./ maximum(result), "k")
ylim(ylims)
tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
savefig("figures/ch3_figure_essentials_signals_correlation_results.svg", dpi=300)


text:
struct PCA
    Scores::Array{Float64,2}
    Loadings::Array{Float64,2}
    Values::Array
    algorithm::String
end

"""
    PCA_NIPALS(X; Factors = minimum(size(X)) - 1, tolerance = 1e-7, maxiters = 200)

Compute's a PCA from `x` using the NIPALS algorithm with a user specified number of latent variables(`Factors`).
The tolerance is the minimum change in the F norm before ceasing execution. Returns a PCA object.

"""
function PCA_NIPALS(X; Factors = minimum(size(X)) - 1, tolerance = 1e-7, maxiters = 200)
    tolsq = tolerance * tolerance
    #Instantiate some variables up front for performance...
    Xsize = size(X)
    Tm = zeros( ( Xsize[1], Factors ) )
    Pm = zeros( ( Factors, Xsize[2] ) )
    t = zeros( ( 1, Xsize[1] ) )
    p = zeros( ( 1, Xsize[2] ) )
    #Set tolerance to floating point precision
    Residuals = copy(X)
    for factor in 1:Factors
        lastErr = sum(abs.(Residuals)); curErr = tolerance + 1;
        t = Residuals[:, 1]
        iterations = 0
        while (abs(curErr - lastErr) > tolsq) && (iterations < maxiters)
            p = Residuals' * t
            p = p ./ sqrt.( p' * p )
            t = Residuals * p
            #Track change in Frobenius norm
            lastErr = curErr
            curErr = sqrt(sum( ( Residuals - ( t * p' ) ) .^ 2))
            iterations += 1
        end
        Residuals -= t * p'
        Tm[:,factor] = t
        Pm[factor,:] = p
    end
    #Find singular values/eigenvalues
    EigVal = sqrt.( LinearAlgebra.diag( Tm' * Tm ) )
    #Scale loadings by singular values
    Tm = Tm * LinearAlgebra.Diagonal( 1.0 / EigVal )
    return PCA(Tm, Pm, EigVal, "NIPALS")
end

"""
    PCA(X; Factors = minimum(size(X)) - 1)

Compute's a PCA from `x` using LinearAlgebra's SVD algorithm with a user specified number of latent variables(`Factors`).
Returns a PCA object.

"""
function PCA(Z; Factors = minimum(size(Z)) - 1)
    svdres = LinearAlgebra.svd(Z)
    return PCA(svdres.U[:, 1:Factors], svdres.Vt[1:Factors, :], svdres.S[1:Factors], "SVD")
end

"""
    (T::PCA)(Z::Array; Factors = length(T.Values), inverse = false)

Calling a PCA object on new data brings the new data `Z` into or out of (`inverse` = true) the PCA basis.

"""
function (T::PCA)(Z::Array; Factors = length(T.Values), inverse = false)
    if (inverse)
        return Z * (Diagonal(T.Values[1:Factors]) * T.Loadings[1:Factors,:])
    else
        return Z * (Diagonal( 1.0 ./ T.Values[1:Factors]) * T.Loadings[1:Factors,:])'
    end
end

struct LDA
    Scores::Array{Float64,2}
    Loadings::Array{Float64,2}
    Values::Array
end

"""
    LDA(X, Y; Factors = 1)

Compute's a LinearDiscriminantAnalysis transform from `x` with a user specified number of latent variables(`Factors`).
Returns an LDA object.
"""
function LDA(X, Y; Factors = 1)
    (Obs, ClassNumber) = size( Y )
    @assert Factors < ClassNumber
    Variables = size( X )[ 2 ]
    #Instantiate some variables...
    ClassMeans = zeros( ClassNumber, Variables )
    ClassSize = zeros( ClassNumber )
    WithinCovariance = zeros( Variables, Variables )
    BetweenCovariance = zeros( Variables, Variables  )
    ClassCovariance = zeros( Variables, Variables  )
    for class in 1 : ClassNumber
        Members = Y[ :, class ] .== 1
        ClassSize[class] = sum(Members)
        ClassMeans[class,:] = StatsBase.mean(X[Members,:], dims = 1)
    end
    GlobalMean = StatsBase.mean(X, dims = 1)

    for class in 1 : ClassNumber
        Members = Y[ :, class ] .== 1
        #calculate the between class covariance matrix
        Diff = ClassMeans[class,:] .- GlobalMean'
        BetweenCovariance .+= ClassSize[class] .* ( Diff * Diff' )
        #calculate the within class covariance matrix
        for member in findall(Members .== true)
            MeanCentered = X[member,:] .- ClassMeans[class, : ]
            WithinCovariance .+= ( MeanCentered * MeanCentered' )
        end
    end
    #Calculate the discriminant axis'
    #eig = LinearAlgebra.eigen(Base.inv(WithinCovariance) * BetweenCovariance)
    eig = LinearAlgebra.eigen(WithinCovariance \ BetweenCovariance)
    if any( imag.( eig.values ) .> 1e-1)
        println("Warning: Some eigenvalues found to have complex contributions > 0.1")
    end
    #Maybe reccomend to the user to do pca first or centerscale or both?
    ReVals = real.(eig.values)
    Sorted = sortperm( ReVals, rev = true)
    Contributions = ReVals[Sorted] .>= 1e-9
    Loadings = real.(eig.vectors[:,Sorted[ Contributions ][1:Factors] ] )
    #Project the X data into the LDA basis
    Scores = X * Loadings
    return LDA( Scores, Loadings, ReVals[ Sorted[ Contributions][1:Factors] ])
end

"""
    ( model::LDA )( Z; Factors = length(model.Values) )

Calling a LDA object on new data brings the new data `Z` into the LDA basis.

"""
function ( model::LDA )( Z; Factors = length(model.Values) )
     Projected = Z * model.Loadings[:,1:Factors]
end

"""
    HLDA(X, YHOT; K = 1, Factors = 1)

Compute's a Hierarchical LinearDiscriminantAnalysis transform from `x` with a user
specified number of latent variables(`Factors`). The adjacency matrices are created
from `K` nearest neighbors.

Returns an LDA object.
*Note: this can be used with any other LDA functions such as Gaussian discriminants
or explained variance*.

Lu D, Ding C, Xu J, Wang S. Hierarchical Discriminant Analysis. Sensors (Basel).
2018 Jan 18;18(1). pii: E279. doi: 10.3390/s18010279.
"""
function HLDA(X, YHOT; K = 1, Factors = 2)
    Dist = SquareEuclideanDistance( X )
    ICAM = OutOfClassAdjacencyMatrix(Dist, YHOT, K)
    ICD = LinearAlgebra.Diagonal(sum( ICAM, dims = 1) )
    ICS = ICD .- ICAM

    Decomp = eigen(X' * ICS * X)
    ReVals = real.(Decomp.values)
    Sorted = sortperm( ReVals, rev = false)
    Loadings = real.(Decomp.vectors[ : , Sorted[ 1:Factors ] ] )

    Dist = SquareEuclideanDistance( X * Loadings[:,1:Factors] )
    OCAM = InClassAdjacencyMatrix(Dist, YHOT, K)
    OCD = LinearAlgebra.Diagonal(sum( OCAM, dims = 1) )
    OCS = OCD .- OCAM

    Decomp = eigen(X' * OCS * X)
    ReVals = real.(Decomp.values)
    Sorted = sortperm( ReVals, rev = true)
    Contributions = ReVals[Sorted] .>= 1e-16
    if Factors > sum(Contributions)
        Factors = sum(Contributions)
    end
    Loadings = real.(Decomp.vectors[ : , Sorted[ Contributions ][ 1 : Factors ] ] )

    return LDA( X * Loadings,  Loadings, ReVals[ Sorted[ Contributions][1:Factors] ] )
end

#Untested...
struct CanonicalCorrelationAnalysis
    U
    V
    r
end

"""
    CanonicalCorrelationAnalysis(A, B)

Returns a CanonicalCorrelationAnalysis object which contains (U, V, r) from Arrays A and B.
Currently Untested for correctness but should compute....
"""
function CanonicalCorrelationAnalysis(A, B)
    (Obs,Vars) = size(A);;
    CAA = (1/Obs) .* A * A'
    CBB = (1/Obs) .* B * B'
    CAB = (1/Obs) .* A * B'
    maxrank = min( LinearAlgebra.rank( A ), LinearAlgebra.rank( B ) )
    CAAInvSqrt = CAA ^ -0.5
    CBBInvSqrt = CBB ^ -0.5
    singvaldecomp = LinearAlgebra.svd( CAAInvSqrt * CAB * CBBInvSqrt )
    Aprime = CAAInvSqrt * singvaldecomp.U[ :,1 : maxrank ]
    Bprime = CAAInvSqrt * singvaldecomp.Vt[ :,1 : maxrank ]
    return CanonicalCorrelationAnalysis(Aprime' * A, Bprime' * B, singvaldecomp.S[1 : maxrank] )
end

"""
    findpeaks( vY; m = 3)

Finds the indices of peaks in a vector vY with a window span of `2m`.
Original R function by Stas_G:(https://stats.stackexchange.com/questions/22974/how-to-find-local-peaks-valleys-in-a-series-of-data)
This version is based on a C++ variant by me.
"""
function findpeaks( vY; m = 3)
    @assert length(size(vY)) == 1
    sze = size(vY)[1];
    (i,q) = (0,0);#generic iterator, second generic iterator
    (lb,rb) = (0,0);#left bound, right bound
    ret = [];
    for i in 1 : ( sze - 2 )
        #Find all regions with negative laplacian between neighbors
        if (sign( vY[ i + 2 ]  - vY[ i + 1 ] ) - sign( vY[ i + 1 ]  - vY[ i ] ) ) < 0
            #Now assess all regions with negative laplacian between neighbors...
            lb = i - m - 1;# define left bound of vector
            lb = (lb < 1) ? 1 : lb
            rb = i + m + 1# define right bound of vector
            rb = (rb >= (sze-2)) ? rb = (sze-3) : rb
            #We have found a peak by our criterion
            if !any(vY[lb:rb] .> vY[i+1])
                push!( ret, i + 1 );
            end #End if found peak
        end #End if laplace condition
    end #End loop
    return ret
end

function RecursiveAlignment(aligned, reference; maxlags::Int = 800, lookahead::Int = 0, minlength::Int = 20, mincorr::Float64 = 0.05)
    if length(aligned) < minlength
        return aligned;
    end
    #Calculate cross correlation between the spectrum and the reference
    lagspan = min( length( aligned ) - 1, maxlags )
    r = StatsBase.crosscor( aligned, reference, -lagspan : lagspan  )
    maxi = argmax( r ) #Find largest value
    if r[ maxi ] <= mincorr # is the largest correlation junk?
        return aligned
    end
    lag = (maxi > lagspan/2) ? maxi - lagspan - 1 : maxi - 1;

    if lag == 0
        if lookahead <= 0
            return aligned
        end
        lookahead -= 1
    end

    #Pad spectra based on lags
    if !( ( lag == 0 ) || ( lag >= length( aligned ) ) )
        if lag > 0
    	   aligned = vcat( ones( lag ) * aligned[ 1 ], aligned[ 1 : (length(aligned) - lag) ] );
        elseif lag < 0
    	   lag = abs( lag );
    	   aligned = vcat( aligned[ ( lag + 1 ) : length( aligned ) ], ones( lag ) * aligned[ end ] );
        end
    end
    #Find new branch point
    middle = round( length( aligned ) / 2) |> Int
    quarter = floor( middle / 4 ) |> Int
    mini = argmin( aligned[ (middle - quarter ) : ( middle + quarter ) ] ) |> Int
    midpnt = mini + middle - quarter
    #Calculate branched components
    LHS = RecursiveAlignment(aligned[ 1 : midpnt ], reference[ 1 : midpnt ];
                        maxlags = maxlags, lookahead = lookahead, minlength = minlength, mincorr = mincorr);
    RHS = RecursiveAlignment(aligned[ (midpnt + 1) : end ], reference[ (midpnt + 1) : end ];
                        maxlags = maxlags, lookahead = lookahead, minlength = minlength, mincorr = mincorr);
    return vcat( LHS, RHS )
end

"""
    RAFFT(raw, reference; maxlags::Int = 500, lookahead::Int = 1, minlength::Int = 20, mincorr::Float64 = 0.05)

RAFFT corrects shifts in the `raw` spectral bands to be similar to those in a given `reference` spectra through
the use of "recursive alignment by FFT". It returns an array of corrected spectra/chromatograms. The number of maximum lags can be
specified, the `lookahead` parameter ensures that additional recursive executions are performed so the first solution
found is not preemptively accepted, the minimum segment length(`minlength`) can also be specified if FWHM are estimable,
and the minimum cross correlation(`mincorr`) for a match can dictate whether peaks were found to align or not.

*Note* This method works best with flat baselines because it repeats last known values when padding aligned spectra.
It is highly efficient, and in my tests does a good job, but other methods definitely exist. Let me know if other peak Alignment
methods are important for your work-flow, I'll see if I can implement them.

Application of Fast Fourier Transform Cross-Correlation for the Alignment of Large Chromatographic and Spectral Datasets
Jason W. H. Wong, Caterina Durante, and, Hugh M. Cartwright. Analytical Chemistry 2005 77 (17), 5655-5661
"""
function RAFFT(raw, reference; maxlags::Int = 500, lookahead::Int = 1, minlength::Int = 20, mincorr::Float64 = 0.05)
    corrected = zeros( size( raw ) )
    for sample in 1 : size( raw )[ 1 ]
        corrected[ sample, : ] = RecursiveAlignment( raw[ sample, : ], reference;
                                                    maxlags = maxlags,
                                                    lookahead = lookahead,
                                                    minlength = minlength,
                                                    mincorr = mincorr)
    end
    return corrected
end

struct COW
   Nodes::Array
   Lengths::Array
end

"""
    COW( A, B; segments = 20, slack = 1,
               maxslack = Int( floor( length(A) / segments ) ) - 2 )

COW makes a CorrelationOptimizedWarping object which warp corrects spectra `A`
to reference `B`. The user can select the number of segments, slack size,
and optionally the maximum slack parameter.

Note: Not fully tested.

"Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping"
Nielsen, N. P. V.; Carstensen, J. M.; Smedsgaard, J.Journal of Chromatography A. 1998,805, 17–35.
"""
function COW( A, B;
                  segments = 20, slack = 1,
                  maxslack = Int( floor( length(A) / segments ) ) - 2 )
   @warn "This code (COW) is mostly untested. Please report any bugs!"
   vars = length( A )
   seglen = Int( floor( vars / segments ) )
   #@assert( vars == length( B ), "spectra should have the same dimensionality" )
   @assert(seglen > slack, "Slack cannot exceed the length of segments ($seglen)")
   @assert(maxslack < (seglen - 2), "Maximum slack cannot exceed the length of segments ($seglen - 2)")
   if seglen < 10
      @warn "the segment length is less then 10 bins. This could be a bad idea!"
   end
   seginit = collect(1 : seglen : vars)
   seginit[end] = vars
   extent = ( (length(seginit) - segments) == 1 ) ? (1:segments) : (1:(segments - 1))
   goalsegments = [ B[ seginit[i]:seginit[(i + 1)] ] for i in extent ]
   goallengths = length.(goalsegments)
   slacks = (-slack):slack
   nodes = copy(seginit)
   corr, oldstate = zeros( segments -1 ), zeros(Int, segments - 2 )
   highscore = -1.1 * segments
   bestslacks = zeros( segments-2 )
   iterate = Iterators.product( repeat( [ slacks ], segments - 2 )... )
   for i in iterate
      #Check change
      curstate = reverse( i ) #The reverse is easier to think about for me...
      difference = findfirst( oldstate .!= curstate )
      #Compute state
      update = seginit[(difference+1):(end-1)] .+ curstate[difference:end]
      nodes[ (difference+1):(end-1) ] = collect(update)
      #Check validity of segments
      if all(nodes[2:(end-1)] .<= vars) && all( abs.(nodes[2:(end-1)] .- seginit[2:(end-1)]) .<= maxslack )
         #Get score
         for d in difference : ( segments - 1 )
            span = nodes[ d ]: (nodes[ d + 1 ] - 1)
            resampled = LinearResample( A[ span ], goallengths[d] )
            corr[ d ] = Statistics.cor( resampled, goalsegments[d] )
         end
         #Store the best score
         if sum(corr) > highscore
            highscore = sum(corr)
            bestslacks = curstate
         end
      end
      #Update state
      oldstate = curstate
   end
   #Warp to best state
   warped = copy( A )
   update = seginit[2:(segments-1)] .+  bestslacks
   nodes = copy(seginit)
   nodes[ 2:(segments-1) ] = update
   return COW(nodes, length.(goalsegments) )
end

"""
    (moo::COW)(A)

Applies a learned warping from a COW object to spectrum A.
"""
function (moo::COW)(A)
   warped = []
   for s in 1 : (length(moo.Nodes) - 1)
      if s == 1
         span = 1 : (moo.Nodes[ s + 1 ] - 1 )
      else
         span = (moo.Nodes[ s ]+1) : (moo.Nodes[ s + 1 ] - 1 )
      end
      fudge = (s == (length(moo.Nodes) - 1)) ? -1 : 0
      warped = vcat(warped, LinearResample( A[ span ], moo.Lengths[ s ] + fudge ) )
   end
   return warped
end

"""
    DynamicTimeWarping( A, B )

Applies Dynamic Time Warping of spectrum `A` so it maps to reference `B`.

Returns the DTW distance, shortest path, and the DTW cost matrix. This does not
automatically rescale `A` to `B`, the way inwhich you want to do that is up to the user
for now.

Sakoe, Hiroaki; Chiba, Seibi (1978). "Dynamic programming algorithm optimization for spoken word recognition".
IEEE Transactions on Acoustics, Speech, and Signal Processing. 26 (1): 43–49. doi:10.1109/tassp.1978.1163055.
"""
function DynamicTimeWarping( A, B )
    #TODO: Add different distances
    DTW = SquareEuclideanDistance( A, B )
    Na, Nb = size( DTW )
    DTW[1,:] = cumsum(DTW[1,:])
    DTW[:,1] = cumsum(DTW[:,1])
    for i = 2:Na, j = 2:Nb
        options = DTW[i-1, j-1]
        if i < Na
            options = vcat( options, DTW[i-1, j] )
        end
        if j < Nb
            options = vcat( options, DTW[i, j-1] )
        end
        DTW[i, j] += minimum(options)
    end
    shortestpath = [(Na, Nb)]
    cost = 0.0
    for step in 2 : Na
        lastposition = shortestpath[end]
        if lastposition[1] == 1
            push!(shortestpath, (1, lastposition[2] - 1) )
        elseif lastposition[1] == 1
            push!(shortestpath, (lastposition[1] - 1, 1) )
        else
            legalmoves = [ lastposition .- 1,
                            (lastposition[1], lastposition[2] - 1),
                            (lastposition[1] - 1, lastposition[2]) ]
            bestmove = argmin( DTW[CartesianIndex.( legalmoves )] )
            push!(shortestpath, legalmoves[bestmove] )
        end
    end
    if shortestpath[end] != (1,1)
        push!(shortestpath, (1,1))
    end
    cost = sum( DTW[ CartesianIndex.( shortestpath ) ] )
    return cost, shortestpath, DTW
end

"""
    AssessHealth( X )

Returns a somewhat detailed Dict containing information about the 'health' of a dataset. What is included is the following:
    - PercentMissing: percent of missing entries (includes nothing, inf / nan) in the dataset
    - EmptyColumns: the columns which have only 1 value
    - RankEstimate: An estimate of the rank of X
    - (optional)Duplicates: returns the rows of duplicate observations
"""
function AssessHealth( X; checkduplicates = true )
    (obs, vars) = size( X )
    NumberMissing = sum( map( x -> isa(x, Missing) || isa(x, Nothing) || isnan(x) || isinf(x), X ) )
    PercentMissing = NumberMissing / (obs * vars)
    empties = []
    for c in 1 : vars
        if length(unique(X[:,c])) == 1
            push!(empties, c)
        end
    end
    RankEst = (NumberMissing == 0) ? LinearAlgebra.rank(X) : missing
    InfoDict = Dict("PercentMissing" => PercentMissing,
                    "EmptyColumns" => empties,
                    "RankEstimate" => RankEst)
    if checkduplicates
        if NumberMissing == 0
            dupes = []
            for o1 in 1:obs
                for o2 in (o1 + 1):obs
                    if X[o1,:] == X[o2,:]
                        push!(dupes, o1)
                    end
                end
            end
            InfoDict["Duplicates"] = dupes
        else
            InfoDict["Duplicates"] = missing
        end
    end
    return InfoDict
end


text:
abstract Enum

Base.typemin{T<:Enum}(x::Type{T}) = T(0)
function Base.show(io::IO, x::Enum)
    T = typeof(x)
    if T == Enum
        print(io, "Enum")
    else
        print(io, string("enum ", T.name, ' ', '(', join(names(T), ", "), ')'))
    end
end
function Base.convert{T<:Integer}(::Type{T},x::Enum)
    if x.n < typemin(T) || x.n > typemax(T)
        error(InexactError())
    end
    convert(T, x.n)
end
function Base.convert{T<:Enum}(::Type{T},x::Integer)
    if x < typemin(T).n || x > typemax(T).n
        error(InexactError())
    end
    T(x)
end
function Base.show{T<:Enum}(io::IO, x::T)
    print(io, T.name, '(')
    first = true
    for s = T
        if s[2] == x.n
            if !first
                print(io, " or ")
            end
            print(io, s[1])
            first = false
        end
    end
    if first
        print(io, "invalid enum constant")
    end
    print(io, ')')
end
Base.names{T<:Enum}(x::Type{T}) = [s[1] for s in T]

macro enum(T,syms...)
    assert(!isempty(syms))
    vals = Array((Symbol,Int64),0)
    i::Int64 = -1
    lo::Int64 = typemax(Int64)
    hi::Int64 = typemin(Int64)
    hasexpr = false
    signed = false
    for s in syms
        if i == typemax(Int64)
            error("@expr enum has a max size of Int64")
        end
        i += 1
        if isa(s,Symbol)
            push!(vals, (s,i))
        elseif isa(s,Expr) && s.head == :(=) && length(s.args) == 2 && isa(s.args[1],Symbol) && isa(s.args[2],Integer)
            if abs(s.args[2]) > typemax(Int64)
                error("@expr argument too large")
            end
            i = s.args[2]
            push!(vals, (s.args[1], i))
            hasexpr = true
        else
            error(string("invalid syntax in @expr macro call: ",s))
        end
        if i < 0 && !signed
            signed = true
        end
        if i < lo
            lo = i
        end
        if i > hi
            hi = i
        end
    end
    if signed
        if max(abs(lo),abs(hi)) < typemax(Int32)
            enumT = Int32
        else
            enumT = Int64
        end
    else
        if max(abs(lo),abs(hi)) < typemax(Uint32)
            enumT = Uint32
        else
            enumT = Uint64
        end
    end
    blk = quote
        immutable $(esc(T)) <: $(esc(Enum))
            n::$(enumT)
            $(esc(T))(n::Integer) = new(n)
        end
        if $lo != 0
            $(esc(:(Base.typemin)))(x::Type{$(esc(T))}) = $(esc(T))($lo)
        end
        $(esc(:(Base.typemax)))(x::Type{$(esc(T))}) = $(esc(T))($hi)
        $(esc(:(Base.start)))(x::Type{$(esc(T))}) = $(esc(:(Base.start)))($vals)
        $(esc(:(Base.next)))(x::Type{$(esc(T))},s) = $(esc(:(Base.next)))($vals,s)
        $(esc(:(Base.done)))(x::Type{$(esc(T))},s) = $(esc(:(Base.done)))($vals,s)
        $(esc(:(Base.length)))(x::Type{$(esc(T))}) = $(length(vals))
        $(esc(:(Base.isless)))(x::$(esc(T)),y::$(esc(T))) = isless(x.n, y.n)
    end
    for (sym,i) in vals
        push!(blk.args, :(const $(esc(sym)) = $(esc(T))($i)))
    end
    push!(blk.args, :nothing)
    blk.head = :toplevel
    return blk
end


text:
using BenchmarkTools

include("../../src/part/isInside.jl")  # loads functions to be tested

function onSegment_TEST()
  p = Point2D(0.0,0.0)
  r = Point2D(1.0,1.0)
  q = Point2D(0.5,0.5)

  #@code_warntype onSegment(p,q,r)
  @btime onSegment($p,$q,$r)

  return  onSegment(p,q,r)
end

function orientation_TEST()
  p   = Point2D(0.0,0.0)
  q   = Point2D(1.0,0.0)
  r   = Point2D(1.0,1.0)
  mid = Point2D(0.5,0.0)

  #@code_warntype orientation(p,q,r)
  @btime orientation($p,$q,$r)

  # CCW
  println("should be 2: ",orientation(p,q,r))

  # CW
  println("should be 1: ",orientation(r,q,p))

  # colinear
  println("should be 0: ",orientation(p,q,mid))

  nothing
end

function  doIntersect_TEST()
  p1 = Point2D(0.0,0.0)
  q1 = Point2D(1.0,0.0)
  p2 = Point2D(0.5,0.5)
  q2 = Point2D(0.5,-0.5)

  #@code_warntype doIntersect(p1,q1,p2,q2)
  @btime doIntersect($p1,$q1,$p2,$q2)

  nothing
end

function isInside_Point2D_TEST()
  n = 4
  p1 = Point2D(0.0,0.0)
  p2 = Point2D(1.0,0.0)
  p3 = Point2D(1.0,1.0)
  p4 = Point2D(0.0,1.0)
  p  = Point2D(rand(),rand())
  myextreme = [100000.0, p[2]]

  polygon = [p1,p2,p3,p4]

  #@code_warntype isInside(polygon, n, p)
  @btime isInside($polygon, $n, $p; extreme=$myextreme)

  return isInside(polygon,n,p)
end
function isInside_TEST()
  n = 4
  p1 = [0.0,0.0]
  p2 = [1.0,0.0]
  p3 = [1.0,1.0]
  p4 = [0.0,1.0]
  p = [50+rand(),rand()]
  myextreme = [100000.0, p[2]]

  polygon = [p1,p2,p3,p4]

  #@code_warntype isInside(polygon, n, p)
  @btime isInside($polygon, $n, $p; extreme=$myextreme)

  return isInside(polygon,n,p)
end

max_stars_repo_name: pseastham/eFEMpart


text:
@testset "Copy" begin
    function test_copy(A, seq)
        @test convert(String, copy(LongSequence{A}(seq))) == seq
    end

    for len in [1, 10, 16, 32, 1000, 10000]
        test_copy(DNAAlphabet{4}, random_dna(len))
        test_copy(RNAAlphabet{4}, random_rna(len))
        test_copy(AminoAcidAlphabet, random_aa(len))

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_copy(DNAAlphabet{2}, random_dna(len, probs))
        test_copy(RNAAlphabet{2}, random_rna(len, probs))
    end

    seq = dna"ACGTACGTACGT"
    subseq = seq[3:6]
    @test copy(subseq) == dna"GTAC"
end

@testset "Concatenation" begin
    function test_concatenation(A, chunks)
        parts = UnitRange{Int}[]
        for i in 1:lastindex(chunks)
            start = rand(1:length(chunks[i]))
            stop = rand(start:length(chunks[i]))
            push!(parts, start:stop)
        end
        str = string([chunk[parts[i]] for (i, chunk) in enumerate(chunks)]...)
        seq = *([LongSequence{A}(chunk)[parts[i]] for (i, chunk) in enumerate(chunks)]...)
        @test convert(String, seq) == uppercase(str)
    end

    for _ in 1:100
        n = rand(1:10)
        chunks = [random_dna(rand(1:100)) for _ in 1:n]
        test_concatenation(DNAAlphabet{4}, chunks)

        chunks = [random_rna(rand(1:100)) for _ in 1:n]
        test_concatenation(RNAAlphabet{4}, chunks)

        chunks = [random_aa(rand(1:100)) for _ in 1:n]
        test_concatenation(AminoAcidAlphabet, chunks)

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        chunks = [random_dna(rand(1:100), probs) for _ in 1:n]
        test_concatenation(DNAAlphabet{2}, chunks)

        chunks = [random_rna(rand(1:100), probs) for _ in 1:n]
        test_concatenation(RNAAlphabet{2}, chunks)
    end
end

@testset "Repetition" begin
    function test_repetition(A, chunk)
        start = rand(1:length(chunk))
        stop = rand(start:length(chunk))
        n = rand(1:10)
        str = chunk[start:stop] ^ n
        seq = LongSequence{A}(chunk)[start:stop] ^ n
        @test convert(String, seq) == uppercase(str)
    end

    for _ in 1:10
        chunk = random_dna(rand(1:100))
        test_repetition(DNAAlphabet{4}, chunk)

        chunk = random_rna(rand(1:100))
        test_repetition(RNAAlphabet{4}, chunk)

        chunk = random_aa(rand(1:100))
        test_repetition(AminoAcidAlphabet, chunk)

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        chunk = random_dna(rand(1:100), probs)
        test_repetition(DNAAlphabet{2}, chunk)

        chunk = random_rna(rand(1:100), probs)
        test_repetition(RNAAlphabet{2}, chunk)
    end
end

@testset "Length" begin
    for len in [0, 1, 2, 3, 10, 16, 32, 1000, 10000]
        seq = LongDNASeq(random_dna(len))
        @test length(seq) === lastindex(seq) === len

        seq = LongRNASeq(random_rna(len))
        @test length(seq) === lastindex(seq) === len

        seq = LongAminoAcidSeq(random_aa(len))
        @test length(seq) === lastindex(seq) === len
    end

    @test length(char"いろはabc") === 6
end

@testset "Access" begin
    dna_seq = dna"ACTG"

    @test dna_seq[1] === DNA_A
    @test dna_seq[2] === DNA_C
    @test dna_seq[3] === DNA_T
    @test dna_seq[4] === DNA_G

    # Access indexes out of bounds
    @test_throws BoundsError dna_seq[-1]
    @test_throws BoundsError dna_seq[0]
    @test_throws BoundsError dna_seq[5]

    @test dna"ACTGNACTGN"[1:5] == dna"ACTGN"
    @test dna"ACTGNACTGN"[5:1] == dna""

    rna_seq = rna"ACUG"
    @test rna_seq[1] === RNA_A
    @test rna_seq[2] === RNA_C
    @test rna_seq[3] === RNA_U
    @test rna_seq[4] === RNA_G

    # Access indexes out of bounds
    @test_throws BoundsError rna_seq[-1]
    @test_throws BoundsError rna_seq[0]
    @test_throws BoundsError rna_seq[5]

    @test rna"ACUGNACUGN"[1:5] == rna"ACUGN"
    @test rna"ACUGNACUGN"[5:1] == rna""

    @test aa"KSAAV"[3] == AA_A
    @test char"いろはにほ"[3] == 'は'
end

@testset "Equality" begin
    @testset "DNA" begin
        a = dna"ACTGN"
        b = dna"ACTGN"
        @test a == b
        @test dna"ACTGN" == dna"ACTGN"
        @test dna"ACTGN" != dna"ACTGA"
        @test dna"ACTGN" != dna"ACTG"
        @test dna"ACTG"  != dna"ACTGN"

        a = dna"ACGTNACGTN"
        b = dna"""
        ACGTN
        ACGTN
        """
        @test a == b
    end

    @testset "RNA" begin
        a = rna"ACUGN"
        b = rna"ACUGN"
        @test a == b
        @test rna"ACUGN" == rna"ACUGN"
        @test rna"ACUGN" != rna"ACUGA"
        @test rna"ACUGN" != rna"ACUG"
        @test rna"ACUG"  != rna"ACUGN"

        a = rna"ACUGNACUGN"
        b = rna"""
        ACUGN
        ACUGN
        """
        @test a == b
    end

    @testset "AminoAcid" begin
        a = aa"ARNDCQEGHILKMFPSTWYVX"
        b = aa"ARNDCQEGHILKMFPSTWYVX"
        @test a == b
        @test a != aa"ARNDCQEGHILKMFPSTWYXV"
        @test a != aa"ARNDCQEGHLKMFPSTWYVX"

        b = aa"""
        ARNDCQEGHI
        LKMFPSTWYV
        X
        """
        @test a == b
    end
end


text:
@testset "Literals" begin
    dv = @data []
    @test isequal(dv, DataArray([], Bool[]))
    @test typeof(dv) == DataVector{Any}

    dv = @data Float64[]
    @test isequal(dv, DataArray(Float64[], Bool[]))
    @test typeof(dv) == DataVector{Float64}

    dv = @data [1, missing, 3]
    @test isequal(dv,
                  DataArray([1, 0, 3],
                            [false, true, false]))

    dv = @data [1 missing 3]
    @test isequal(dv,
                  DataArray([1 0 3],
                            [false true false]))

    dv = @data Float64[1, missing, 3]
    @test isequal(dv,
                  DataArray(Float64[1, 0, 3],
                            [false, true, false]))
    @test typeof(dv) == DataVector{Float64}

    dv = @data Float64[1 missing 3]
    @test isequal(dv,
                  DataArray(Float64[1 0 3],
                            [false true false]))
    @test typeof(dv) == DataMatrix{Float64}

    dv = @data [missing, missing]
    @test isequal(dv, DataArray(Any, 2))
    @test typeof(dv) == DataVector{Any}

    dv = @data [missing missing]
    @test isequal(dv, DataArray(Any, 1, 2))
    @test typeof(dv) == DataMatrix{Any}

    dm = @data [1 missing; 3 4]
    @test isequal(dm,
                  DataArray([1 0; 3 4],
                            [false true; false false]))

    dm = @data Float64[1 missing; 3 4]
    @test isequal(dm,
                  DataArray(Float64[1 0; 3 4],
                            [false true; false false]))
    @test typeof(dm) == DataMatrix{Float64}

    dm = @data [missing missing; missing missing]
    @test isequal(dm, DataArray(Any, 2, 2))
    @test typeof(dm) == DataMatrix{Any}

    pdv = @pdata [1, missing, 3]
    @test isequal(pdv,
                  PooledDataArray([1, 0, 3],
                                  [false, true, false]))

    pdv = @pdata Float64[1, missing, 3]
    @test isequal(pdv,
                  PooledDataArray(Float64[1, 0, 3],
                                  [false, true, false]))
    @test typeof(pdv) == PooledDataArray{Float64,UInt32,1}

    pdv = @pdata [1 missing 3]
    @test isequal(pdv,
                  PooledDataArray([1 0 3],
                                  [false true false]))

    pdv = @pdata Float64[1 missing 3]
    @test isequal(pdv,
                  PooledDataArray(Float64[1 0 3],
                                  [false true false]))
    @test typeof(pdv) == PooledDataArray{Float64,UInt32,2}

    pdm = @pdata [1 missing; 3 4]
    @test isequal(pdm,
                  PooledDataArray([1 0; 3 4],
                                  [false true; false false]))

    pdm = @pdata Float64[1 missing; 3 4]
    @test isequal(pdm,
                  PooledDataArray(Float64[1 0; 3 4],
                                  [false true; false false]))
    @test typeof(pdm) == PooledDataArray{Float64,UInt32,2}

    pdm = @pdata [1 missing;
                3 4]
    @test isequal(pdm,
                  PooledDataArray([1 0; 3 4],
                                  [false true; false false]))

    dv1 = @data zeros(4)
    dv2 = @data ones(4)
    dv3 = @data rand(4)

    dm1 = @data zeros(4, 4)
    dm2 = @data ones(4, 4)
    dm3 = @data rand(4, 4)

    pdv1 = @pdata zeros(4)
    pdv2 = @pdata ones(4)
    pdv3 = @pdata rand(4)

    pdm1 = @pdata zeros(4, 4)
    pdm2 = @pdata ones(4, 4)
    pdm3 = @pdata rand(4, 4)

    mixed1 = @data ["x", 1, 1.23, missing]
    mixed2 = @data [missing, "x", 1, 1.23, missing]

    @test isequal(mixed1, DataArray(Any["x", 1, 1.23, 0],
                                    [false, false, false, true]))
    @test isequal(mixed2, DataArray(Any[missing, "x", 1, 1.23, 0],
                                    [true, false, false, false, true]))

    x = 5.1
    ex = :([1, 2, 3])
    DataArrays.parsedata(ex)
    @test isequal(@data([1, 2, 3]),
                  DataArray([1, 2, 3], [false, false, false]))

    ex = :([1, 2, 3.0])
    DataArrays.parsedata(ex)
    @test isequal(@data([1, 2, 3.0]),
                  DataArray([1, 2, 3.0], [false, false, false]))

    ex = :([1, 2, x])
    DataArrays.parsedata(ex)
    @test isequal(@data([1, 2, x]),
                  DataArray([1, 2, x], [false, false, false]))

    ex = :([1, 2, missing])
    DataArrays.parsedata(ex)
    @test isequal(@data([1, 2, missing]),
                  DataArray([1, 2, 1], [false, false, true]))

    ex = :([1, 2, x, missing])
    DataArrays.parsedata(ex)
    @test isequal(@data([1, 2, x, missing]),
                  DataArray([1, 2, x, 1], [false, false, false, true]))

    # Matrices
    ex = :([1 2; 3 4])
    DataArrays.parsedata(ex)
    @data([1 2; 3 4])
    @test isequal(@data([1 2; 3 4]),
                  DataArray([1 2; 3 4], [false false; false false]))

    ex = :([1 2; 3.0 4])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 2; 3.0 4]),
                  DataArray([1 2; 3.0 4], [false false; false false]))

    ex = :([1 2; x x])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 2; x x]),
                  DataArray([1 2; x x], [false false; false false]))

    ex = :([1 2; missing missing])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 2; missing missing]),
                  DataArray([1 2; 1 1], [false false; true true]))

    ex = :([1 2; x missing])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 2; x missing]),
                  DataArray([1 2; x 1], [false false; false true]))

    # Complex vector expressions
    ex = :([1 + 1, 2 + 2, x * x, missing])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 + 1, 2 + 2, x * x, missing]),
                  DataArray([1 + 1, 2 + 2, x * x, 1],
                            [false, false, false, true]))

    ex = :([sin(1), cos(2) + cos(2), exp(x * x), sum([1, 1, 1])])
    DataArrays.parsedata(ex)
    @test isequal(@data([sin(1),
                         cos(2) + cos(2),
                         exp(x * x),
                         sum([1, 1, 1])]),
                  DataArray([sin(1),
                             cos(2) + cos(2),
                             exp(x * x),
                             sum([1, 1, 1])],
                            [false, false, false, false]))

    ex = :([1 + 1im, 2 + 2im])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 + 1im, 2 + 2im]),
                  DataArray([1 + 1im, 2 + 2im],
                            [false, false]))

    # Complex matrix expressions
    ex = :([1 + 1 2 + 2; x * x missing])
    DataArrays.parsedata(ex)
    @test isequal(@data([1 + 1 2 + 2; x * x missing]),
                  DataArray([1 + 1 2 + 2; x * x 1],
                            [false false; false true]))

    ex = :([sin(1) cos(2) + cos(2);
            exp(x * x) sum([1, 1, 1])])
    DataArrays.parsedata(ex)
    @test isequal(@data([sin(1) cos(2) + cos(2);
                        exp(x * x) sum([1, 1, 1])]),
                  DataArray([sin(1) cos(2) + cos(2);
                             exp(x * x) sum([1, 1, 1])],
                            [false false;
                             false false]))

    @test isequal(DataArrays.fixargs(:([1, 2, missing, x]).args, -1),
                  (Any[1, 2, -1, :x], Any[false, false, true, false]))

    @test isequal(DataArrays.findstub_vector(:([1, 2, missing, x])), 1)
    @test isequal(DataArrays.findstub_vector(:([missing, missing, missing, x])), :x)

    # Lots of variables
    a, b, c, d = 1, 2, 3, 4
    @data [a, b, c, d]
end


text:
############################################################
# joMatrix - miscaleneous constructor-only operators #######
############################################################

# identity operators: joDirac
include("joMatrixConstructors/joDirac.jl")

# identity operators: joEye
include("joMatrixConstructors/joEye.jl")

# diagonal operators: joDiag
include("joMatrixConstructors/joDiag.jl")

# matrix of ones
include("joMatrixConstructors/joOnes.jl")

# matrix of zeros
include("joMatrixConstructors/joZeros.jl")

# matrix of constants
include("joMatrixConstructors/joConstants.jl")

# vector conversion operators: joReal joImag joConj
include("joMatrixConstructors/joVecConvert.jl")

# sinc interpolation
include("joMatrixConstructors/joSincInterp.jl")

# 1D linear interpolation
include("joMatrixConstructors/joLinInterp1D.jl")


text:
# run Julia REPL with: JULIA_NUM_THREADS=4 julia
module SingletonExample

export AppKey

# AppKey contains an app id and encryption key
struct AppKey
    appid::String
    value::UInt128
end

# placeholder for AppKey object.  Start with nothing until initialized.
const appkey = Ref{AppKey}()

# initializer
function construct()
    global appkey
    if !isassigned(appkey)
        ak = AppKey("myapp", rand(UInt128))
        println("constructing $ak")
        appkey[] = ak
    end
    return nothing
end

function test()
    construct()
end

function test_multithreading()
    println("Number of threads: ", Threads.nthreads())
    Threads.@threads for i in 1:8
        construct()
    end
end

end #module

# first test
using .SingletonExample
SingletonExample.test()

# Quit REPL, evaluate the module above, and run this test function.
SingletonExample.test_multithreading()

# test
#=
julia> SingletonExample.test()
constructing AppKey("myapp", 0x2b20f8bf74d7175ec25ad70669faa434)

julia> SingletonExample.test_multithreading()
Number of threads: 4
constructing AppKey("myapp", 0xaa0c9aa1746439568ec7f2ff99edb107)
constructing AppKey("myapp", 0xb9078c91ce3e764322d4d5c980b011d0)
constructing AppKey("myapp", 0x23c3507affbd65cab8c9fb25d5c3db29)
constructing AppKey("myapp", 0x3e3f51b5fecc028c1007b8e182fa2d74)
=#

# How do we ensure a single construction???

# run Julia REPL with: JULIA_NUM_THREADS=4 julia
module SingletonExample2

export AppKey

# object that needs to be initialized as appkey
struct AppKey
    appid::String
    value::UInt128
end

# placeholder for AppKey object.  Start with nothing until initialized.
const appkey = Ref{AppKey}()

# create a ReentrantLock
const appkey_lock = Ref(ReentrantLock())

# change construct() function to acquire lock before
# construction, and release it after it's done.
function construct()
    global appkey
    global appkey_lock
    lock(appkey_lock[])
    try
        if !isassigned(appkey)
            ak = AppKey("myapp", rand(UInt128))
            println("constructing $ak")
            appkey[] = ak
        else
            println("skipped construction")
        end
    finally
        unlock(appkey_lock[])
        return appkey[]
    end
end

function test_multithreading()
    println("Number of threads: ", Threads.nthreads())
    Threads.@threads for i in 1:8
        construct()
    end
end

end #module

using .SingletonExample2
SingletonExample2.test_multithreading()

#=
julia> F.test_multithreading()
Number of threads: 4
constructing appkey object Main.F.AppKey(0xb8)
skipped construction.
skipped construction.
skipped construction.
skipped construction.
skipped construction.
skipped construction.
skipped construction.
=#


text:
module PlutoGrid

using HypertextLiteral: @htl, JavaScript
using Tables
using DataFrames

export readonly_table, editable_table, create_dataframe

const MODIFIED_COL_INDICATOR = "modified_column_"
const MODIFIED_COL_REGEX = r"^(?!modified_column_).+"

"""
	readonly_table(df; sortable=true, filterable=true, pagination=false)

Shows a non-editable table in Pluto.

`df`: DataFrame or any Tables.jl compatible data source (the latter is internally converted to a DataFrame)

`sortable`: enables sorting in the table GUI for all columns (enabled by default)

`filterable`: enables filtering in the table GUI for all columns (enabled by default)

`resizable`: allows resizing of column widths in GUI (enabled by default)

`pagination`: enables pagination of the table (disabled by default)

`height`: vertical size of the table in Pluto in pixel (default: 600)

"""
function readonly_table(df:: DataFrame; filterable:: Bool=true, kwargs...)
	column_defs = _make_col_defs(df; filterable)
	data = _prepare_data(df)
    return _create_table(column_defs, data; filterable, kwargs...)
end

readonly_table(df; kwargs...) = readonly_table(DataFrame(df); kwargs...)

"""
	@bind edits editable_table(df, [editable_cols]; sortable=true, filterable=true, pagination=false, height=600, return_only_modified=false)

Shows an editable table in Pluto. In case of user edits in the table,

`df`: DataFrame or any Tables.jl compatible data source (the latter is internally converted to a DataFrame)

`editable_cols`: (optional) Names of columns which are editable. If not given, all columns are editable.

`sortable`: enables sorting in the table GUI for all columns (enabled by default)

`filterable`: enables filtering in the table GUI for all columns (enabled by default)

`resizable`: allows resizing of column widths in GUI (enabled by default)

`pagination`: enables pagination of the table (disabled by default)

`height`: vertical size of the table in Pluto in pixel (default: 600)

"""
function editable_table(df:: DataFrame, editable_cols:: AbstractVector{<: AbstractString}=collect(names(df)); filterable:: Bool=true, return_only_modified:: Bool=false, kwargs...)
	column_defs = _make_col_defs(df; filterable, editable_cols)
	data = _prepare_data(df)
	return _create_table(column_defs, data; editable=true, filterable, kwargs...)
end

editable_table(df, editable_cols; kwargs...) = editable_table(DataFrame(df), editable_cols; kwargs...)
editable_table(df; kwargs...) = editable_table(DataFrame(df); kwargs...)

function _create_table(column_defs:: AbstractVector{<: AbstractDict}, data:: AbstractVector;
	sortable=true, filterable=true, resizable=true, pagination=false, height:: Integer=600,
	editable=false, insert=true, delete=true, auto_confirm=false)

	edit_button = @htl("""
		<button
			id="update_grid"
			type="button">
				Confirm
		</button>
	""")

	edit_button_callbacks = JavaScript("""
	div.querySelector("button#update_grid").addEventListener("click", (e) => {
		div.value = rowData; // return complete table
		div.dispatchEvent(new CustomEvent("input"));
		e.preventDefault();
		div.querySelector("button#update_grid").style.background='green';
		})
	""")

	insert_button = @htl("""
		<button
		id="insert_row"
		type="button">
			Insert Row
		</button>
	""")

	insert_new_row_callback = JavaScript("""
	div.querySelector("button#insert_row").addEventListener("click", (e) => {
		const row = Object.assign({}, rowData[rowData.length - 1]);
		gridOptions.rowData.push(row);
		gridOptions.api.setRowData(gridOptions.rowData);
		const data = gridOptions.rowData[gridOptions.rowData.length -1]
		gridOptions.columnApi.getAllColumns().forEach(col => {
			data[$(MODIFIED_COL_INDICATOR) + col.colId] = true;
			});

		gridOptions.api.refreshCells({force: true});

		div.querySelector("button#update_grid").style.background='red';
		})
	""")

	delete_button = @htl("""
		<button
		id="delete_row"
		type="button">
			Delete Row
		</button>
	""")

	delete_row_callback = JavaScript("""
	div.querySelector("button#delete_row").addEventListener("click", (e) => {

		const selectedRow = gridOptions.api.getFocusedCell()
		const id = gridOptions.rowData[selectedRow.rowIndex].i

		gridOptions.rowData.splice(selectedRow.rowIndex, 1)
		gridOptions.api.setRowData(gridOptions.rowData)

		div.querySelector("button#update_grid").style.background='red';
		})
	""")

	checkbox_renderer = JavaScript("""
	// source: https://stackoverflow.com/a/62173238/14693778
	function CheckboxRenderer() {}

		CheckboxRenderer.prototype.init = function(params) {
		  this.params = params;

		  this.eGui = document.createElement('input');
		  this.eGui.type = 'checkbox';
		  this.eGui.checked = params.value;

		  this.checkedHandler = this.checkedHandler.bind(this);
		  this.eGui.addEventListener('click', this.checkedHandler);
		}

		CheckboxRenderer.prototype.checkedHandler = function(e) {
		  let checked = e.target.checked;
		  let colId = this.params.column.colId;
		  this.params.node.setDataValue(colId, checked);
		}

		CheckboxRenderer.prototype.getGui = function(params) {
		  return this.eGui;
		}

		CheckboxRenderer.prototype.destroy = function(params) {
		  this.eGui.removeEventListener('click', this.checkedHandler);
		}
	""")

	return @htl("""
<div id="myGrid" style="height: $(height)px;" class="ag-theme-alpine">
$(editable ? edit_button : "")
$((editable && insert) ? insert_button : "")
$((editable && delete) ? delete_button : "")

<script src="https://unpkg.com/ag-grid-community/dist/ag-grid-community.min.js"></script>
<script>

function numberParser(params) {
	return Number(params.newValue);
};

var div = currentScript.parentElement;
// set default output value
div.value = null;

const columnDefs = $(column_defs);
const rowData = $(_transfer_data(data));

$checkbox_renderer

$(editable ? edit_button_callbacks : JavaScript(""))
$((editable && insert) ? insert_new_row_callback : JavaScript(""))
$((editable && delete) ? delete_row_callback : JavaScript(""))

// let the grid know which columns and what data to use
const gridOptions = {
  columnDefs: columnDefs,
  rowData: rowData,
  components: { checkboxRenderer: CheckboxRenderer },
  defaultColDef: {
    filter: $(sortable),
    sortable: $(filterable),
	resizable: $(resizable),
	cellStyle: params => {
		// source: https://stackoverflow.com/questions/65273946/ag-grid-highlight-cell-logic-not-working-properly
      if (
        params.data[$(MODIFIED_COL_INDICATOR)+ params.column.colDef.field]
      ) {
        return { 'color': 'red', 'background-color': 'yellow' };
      } else {
        return null;
      }
    },
  },
  pagination: $(pagination),
  undoRedoCellEditing: true,
  onCellValueChanged: (params) => {
		if (params.oldValue === params.newValue) {
		return;
		}
		const column = params.column.colDef.field;
		params.data[$(MODIFIED_COL_INDICATOR) + column] = true;
		params.api.refreshCells({
		force: true,
		columns: [column],
		rowNodes: [params.node]
		});
		div.querySelector("button#update_grid").style.background='red';
	},
	onCellKeyPress: (event) => {
		if ((event.event.keyCode == 13) && (event.event.shiftKey == true)) {
			div.querySelector("button#update_grid").click();
		}
	},
};
new agGrid.Grid(div, gridOptions);

$(auto_confirm ? JavaScript("""div.querySelector("button#update_grid").click();""") : JavaScript(""))

</script>
</div>
""")
end

function _make_col_defs(df; filterable=true, editable_cols=String[])
	setdiff(editable_cols, names(df)) == [] || error("not all columns defined as editable are in input data")

	column_defs = Dict[]
	for c in names(df)
		col_dict = Dict{String, Any}("field" => c)

		col_is_editable = c ∈ editable_cols
		col_dict["editable"] = col_is_editable

		# special types and filters for specific element types
		if eltype(df[!, c]) <: Bool
			col_dict["cellRenderer"] = JavaScript("""'checkboxRenderer'""")

		elseif eltype(df[!, c]) <: Number
			col_dict["type"] = "numericColumn"
			filterable && (col_dict["filter"] = "agNumberColumnFilter")
			col_is_editable && (col_dict["valueParser"] = JavaScript("numberParser"))
		end
		# ToDo: add type / filter for dates

		push!(column_defs, col_dict)
	end
    return column_defs
end

_prepare_data(df) = [NamedTuple(row) for row in Tables.rows(df)]

_transfer_data(data) = if isdefined(Main, :PlutoRunner) && isdefined(Main.PlutoRunner, :publish_to_js)
	# faster data transfer using MsgPack
	JavaScript(Main.PlutoRunner.publish_to_js(data))
else
	data
end

"""
	create_dataframe(x)

Utility function to create a DataFrame from a `@bind` variable of a `editable_table`.
"""
create_dataframe(::Nothing; kwargs...) = DataFrame()
create_dataframe(::Missing; kwargs...) = DataFrame()

function create_dataframe(x:: Dict; drop_modified_indicator=true)
	df = DataFrame(x)
	drop_modified_indicator ? df[:, MODIFIED_COL_REGEX] : df
end

function create_dataframe(x:: AbstractVector; drop_modified_indicator=true)
	df = DataFrame()
	for row in x
		push!(df, row; cols=:union)
	end
	drop_modified_indicator ? df[:, MODIFIED_COL_REGEX] : df
end

# precompilation
let
	df = DataFrame(x=1:10, y=10:-1:1)
	readonly_table(df)
	editable_table(df)
end

end


text:
# ######################
# A DataFrame way to get file

import Base: getindex, view

export getindex

function Base.getindex(file::JDFFile, rows, col::String)
    JDF.load(file; cols = [col])[rows, 1]
end

Base.getindex(file::JDFFile, rows, cols::AbstractVector{String}) = begin
    JDF.load(file; cols = cols)[rows, :]
end

Base.view(file::JDFFile, rows, cols) = getindex(file, rows, cols)

getindex(file::JDFFile, rows, cols) = JDF.load(file)[rows, cols]

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using jFEMTools
import Tensors: Vec,  ⋅
using SparseArrays
import WriteVTK
const jF = jFEMTools

# We solve a simple 3D Poisson PDE with homogeneous Dirichlet boundary conditions
# using ContinuousLagrange finite elements on a Hexahedral mesh

# Problem:

# Δu = f  in Ω
# u  = 0  in ∂Ω

# We start  generating a simple grid with 20x20 triangular elements
# using `unitSquareMesh2`. The generator defaults to the unit square,
# so we don't need to specify the corners of the domain.
mesh = hyper_rectagle_mesh2(HexahedronCell,(10,10,10))

# ### Initiate function Spaces
dim = jF.getdim(mesh)
P1 = ContinuousLagrange(:Hexahedron,1)
Wh = FEMFunctionSpace(mesh, P1, 1)

# Declare variables
u_h = TrialFunction(Wh)
#v_h = TestFunction(Wh)

# ### Degrees of freedom
# Next we need to define a `DofHandler`, which will take care of numbering
# and distribution of degrees of freedom for our approximated fields.
# We create the `DofHandler` and then add a single field called `u_h`.
dh = DofHandler(mesh,[u_h])

# Now that we have distributed all our dofs we can create our tangent matrix,
# using `create_sparsity_pattern`. This function returns a sparse matrix
# with the correct elements stored.
K = jF.create_sparsity_pattern(dh);

# ### Boundary conditions
dbc = Dirichlet(dh, u_h, "boundary", 0.0)

# ### RHS function
f(x::Vec{dim}) where {dim} = 2*π^2*sin(π*x[1])*sin(π*x[2])*sin(π*x[3])

# ### Assembling the linear system
# Now we have all the pieces needed to assemble the linear system, $K u = f$.
function doassemble(Wh, K::SparseMatrixCSC, dh::jF.DofHandler)
  # We allocate the element stiffness matrix and element force vector
  # just once before looping over all the cells instead of allocating
  # them every time in the loop.
  n_basefuncs = jF.getnbasefunctions(Wh)
  Ke = zeros(n_basefuncs, n_basefuncs)
  fe = zeros(n_basefuncs)
  b = zeros(jF.ndofs(dh))
  cell_dofs = Vector{Int}(undef, jF.ndofs_per_cell(dh))

  # Next we define the global force `f` and
  # create an assembler. The assembler
  # is just a thin wrapper around `f` and `K` and some extra storage
  # to make the assembling faster.
  assembler = jF.start_assemble(K, b)

  # It is now time to loop over all the cells in our mesh
  @inbounds for (cellcount, cell) in enumerate(CellIterator(mesh))
      # Always remember to reset the function space local data
      fill!(Ke, 0)
      fill!(fe, 0)
      jF.reinit!(Wh, cell)

      # It is now time to loop over all the quadrature points in the cell and
      # assemble the contribution to `Ke` and `fe`. The integration weight
      # can be queried from `cellvalues` by `getdetJdV`.
      for q_point in 1:jF.getnquadpoints(Wh)
          dΩ = jF.getdetJdV(Wh, q_point)
          fh = jF.function_value(f, Wh, cell, q_point)
          # For each quadrature point we loop over all the (local) shape functions.
          # We need the value and gradient of the testfunction `v` and also the gradient
          # of the trial function `u`.
          for i in 1:n_basefuncs
              v  = jF.shape_value(Wh, q_point, i)
              ∇v = jF.shape_gradient(Wh, q_point, i)
              fe[i] += fh*v * dΩ
              for j in 1:n_basefuncs
                  ∇u = jF.shape_gradient(Wh, q_point, j)
                  Ke[i, j] += (∇v ⋅ ∇u) * dΩ
              end
          end
      end
      # The last step in the element loop is to assemble `Ke` and `fe`
      # into the global `K` and `f` with `assemble!`.
      jF.celldofs!(cell_dofs, dh, cell)
      jF.assemble!(assembler, cell_dofs, fe, Ke)
  end
  return K, b
end

# ### Solution of the system
# The last step is to solve the system. First we call `doassemble`
# to obtain the global stiffness matrix `K` and force vector `f`.
K, b = doassemble(Wh, K, dh);

# To account for the boundary conditions we use the `apply!` function.
# This modifies elements in `K` and `f` respectively, such that
# we can get the correct solution vector `u` by using `\`.
apply!(K,b,dbc);
u = K \ b;

# Plot approximation
vi = jFEMTools.vertexdofs(dh, u_h);
vtk_file = vtk_grid("poisson3D", mesh)
vtk_file["u", WriteVTK.VTKPointData()] = u[vi]
outfiles = WriteVTK.vtk_save(vtk_file)

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents, Unitful
import Unitful: Time, ms, s, minute, hr

resetClock!(𝐶)

clk = Clock(t0=60, unit=s)
@test tau(clk) == 60s
@test clk.time == 60
@test clk.Δt == 0.01
@test clk.unit == s

resetClock!(clk)
@test tau(clk) == 0
@test clk.Δt == 0.01
@test clk.unit == NoUnits

clk = Clock(t0=60, unit=s)
@test_warn "deleted time unit" sync!(clk)
@test tau(clk) == 0
@test clk.Δt == 0.01
@test clk.unit == NoUnits

a = [0.0]
incr!(x) = x[1] += 1
event!(clk, fun(incr!, a), every, 1)
DiscreteEvents.incr!(clk)
@test a[1] == 1
@test clk.time == 1
event!(clk, fun(DiscreteEvents.stop!, clk), 5.1)
clk.state = DiscreteEvents.Undefined()
run!(clk, 9)
@test clk.state == DiscreteEvents.Halted()
@test a[1] == 5
@test clk.time == 5.1
@test clk.end_time == 10
resume!(clk)
@test clk.time == 10
@test a[1] == 10

@test_warn  "undefined transition" DiscreteEvents.init!(clk)
clk.state = DiscreteEvents.Undefined()
DiscreteEvents.init!(clk)
@test clk.state == DiscreteEvents.Idle()


text:
using Test, DataStructures

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/Classic.jl")
    using .Classic
else
    # invoked during test
    using Classic
end

"""
    init!(g::Classic.Graph{V,E}) where {V,E<:Classic.Edge}

Create an example graph with unweighted edges
"""
function init!(g::Classic.Graph{V,E}) where {V,E<:Classic.Edge}
    #cc = Classic

    cities = Vector{String}(["Seattle", "San Francisco",
    "Los Angeles", "Riverside", "Phoenix", "Chicago", "Boston",
    "New York", "Atlanta", "Miami", "Dallas", "Houston", "Detroit",
    "Philadelphia", "Washington",
    ])

    append!(g,cities)

    add!(g, "Seattle", "Chicago")
    add!(g, "Seattle","San Francisco")
    add!(g, "San Francisco", "Riverside")
    add!(g, "San Francisco", "Los Angeles")
    add!(g, "Los Angeles", "Riverside")
    add!(g, "Los Angeles", "Phoenix")
    add!(g, "Riverside", "Phoenix")
    add!(g, "Riverside", "Chicago")
    add!(g, "Phoenix", "Houston")
    add!(g, "Dallas", "Chicago")
    add!(g, "Dallas", "Atlanta")
    add!(g, "Dallas", "Houston")
    add!(g, "Houston", "Atlanta")
    add!(g, "Houston", "Miami")
    add!(g, "Atlanta", "Chicago")
    add!(g, "Atlanta", "Washington")
    add!(g, "Atlanta", "Miami")
    add!(g, "Miami", "Washington")
    add!(g, "Chicago", "Detroit")
    add!(g, "Detroit", "Boston")
    add!(g, "Detroit", "Washington")
    add!(g, "Detroit", "Washington")
    add!(g, "Detroit", "New York")
    add!(g, "Boston", "New York")
    add!(g, "New York", "Philadelphia")
    add!(g, "Philadelphia", "Washington")

    return cities
end

"""
    init!(g::Classic.WeightedGraph{V,E}) where {V,E<:Classic.WeightedEdge}

Create a example graph with weighted edges
"""
function init!(g::Classic.WeightedGraph{V,E}) where {V,E<:Classic.WeightedEdge}
    #cc = Classic

    cities = Vector{String}(["Seattle", "San Francisco",
    "Los Angeles", "Riverside", "Phoenix", "Chicago", "Boston",
    "New York", "Atlanta", "Miami", "Dallas", "Houston", "Detroit",
    "Philadelphia", "Washington",
    ])

    append!(g,cities)

    add!(g, "Seattle", "Chicago", 1737);
    add!(g, "Seattle", "San Francisco", 678);
    add!(g, "San Francisco", "Riverside", 386);
    add!(g, "San Francisco", "Los Angeles", 348);
    add!(g, "Los Angeles", "Riverside", 50);
    add!(g, "Los Angeles", "Phoenix", 357);
    add!(g, "Riverside", "Phoenix", 307);
    add!(g, "Riverside", "Chicago", 1704);
    add!(g, "Phoenix", "Dallas", 887);
    add!(g, "Phoenix", "Houston", 1015);
    add!(g, "Dallas", "Chicago", 805);
    add!(g, "Dallas", "Atlanta", 721);
    add!(g, "Dallas", "Houston", 225);
    add!(g, "Houston", "Atlanta", 702);
    add!(g, "Houston", "Miami", 968);
    add!(g, "Atlanta", "Chicago", 588);
    add!(g, "Atlanta", "Washington", 543);
    add!(g, "Atlanta", "Miami", 604);
    add!(g, "Miami", "Washington", 923);
    add!(g, "Chicago", "Detroit", 238);
    add!(g, "Detroit", "Boston", 613);
    add!(g, "Detroit", "Washington", 396);
    add!(g, "Detroit", "New York", 482);
    add!(g, "Boston", "New York", 190);
    add!(g, "New York", "Philadelphia", 81);
    add!(g, "Philadelphia", "Washington", 123);

    return cities
end

function test_SimpleEdge(io::IO)
    #cc = Classic

    e = SimpleEdge(2, 7)
    @test 2 == e[1]
    @test 7 == e[2]
    @test "2 -> 7" == string(e)
    @test 2 == e[1]
    @test 7 == e[2]

    #r = Classic.reverse(e)
    r = Base.reverse(e)
    @test e[1] == r[2]
    @test e[2] == r[1]
end

function test_UnweightedGraph(io::IO)
    #cc = Classic

    @test_throws TypeError UnweightedGraph{Int64,Int64}()
    @test UnweightedGraph{Int64,SimpleEdge}() isa Any

    # V::Int64, E::Edge
    g = UnweightedGraph{Int64,SimpleEdge}()
    @test isempty(g.vertices)
    @test isempty(g.edges_lists)

    g = UnweightedGraph{Int64,SimpleEdge}([1,2,3])
    @test 3 == length(g.vertices)
    @test [1,2,3] == g.vertices
    @test 3 == length(g.edges_lists)

    # V::String, E::Edge
    g = UnweightedGraph{String,SimpleEdge}()
    @test isempty(g.vertices)
    @test isempty(g.edges_lists)

    g = UnweightedGraph{String,SimpleEdge}(["apple", "orange", "guava"])
    @test 3 == length(g.vertices)
    @test ["apple", "orange", "guava"] == g.vertices
    @test 3 == length(g.edges_lists)

    # Learn types
    @test UnweightedGraph{Int64,Edge} <: Graph{Int64,Edge}
    @test isa(UnweightedGraph{Int64,Edge}(), Graph{Int64,Edge})
end

"""

Sovle shortest route of `UnweightedGraph`
"""
function test_bfs(io::IO)
    #cc = Classic

    g = UnweightedGraph{String,SimpleEdge}()
    cities = init!(g)
    @test length(cities) == length(g.vertices)
    @test length(cities) == length(g.edges_lists)

    for edges in g.edges_lists
        @test !isempty(edges)
    end

    #show(g)

    # `Classic.` is necessary because `is_goal` is ambiguous.
    is_goal(pt) = Classic.is_goal("Miami", pt)
    neighbor_of(v) = Classic.neighbor_of(g, v)

    node = bfs("Boston", is_goal, neighbor_of)
    @test ["Boston", "Detroit", "Washington", "Miami"] == node_to_path(node)
    #println("Shortest route: $(node_to_path(node))")

    nothing
end

function test_WeightedEdge(io::IO)
    #cc = Classic

    e = WeightedEdge((3,7),13.0)
    @test 3 == e[1]
    @test 7 == e[2]
    @test 13 == weight(e)

    e = WeightedEdge(3,7,13.0)
    @test 3 == e[1]
    @test 7 == e[2]
    @test 13 == weight(e)

    #e = Classic.reverse(e)
    e = reverse(e)
    @test 3 == e[2]
    @test 7 == e[1]
    @test 13 == weight(e)

    # sum() and +operator
    edges = Vector{WeightedEdge}([
        WeightedEdge(3,7,13.0),
        WeightedEdge(4,8,17.0),
        WeightedEdge(5,9,19.0)
    ])
    @test 49 == sum(edges)

end

function test_WeightedGraph(io::IO)
    #cc = Classic

    @test_throws TypeError Classic.WeightedGraph{Int64,Int64}()
    @test WeightedGraph{Int64,WeightedEdge}() isa Any # @test_nothrows

    g = WeightedGraph{Int64,WeightedEdge}([1,2,3])
    @test 3 == length(g.vertices)
    @test [1,2,3] == g.vertices
    @test 3 == length(g.edges_lists)

    g = WeightedGraph{String,WeightedEdge}(["apple", "orange", "guava"])
    @test 3 == length(g.vertices)
    @test ["apple", "orange", "guava"] == g.vertices
    @test 3 == length(g.edges_lists)

    add!(g, "apple", "orange", 13)
    @test 1 == length(g.edges_lists[1])
    @test 1 == length(g.edges_lists[2])
    @test 0 == length(g.edges_lists[3])
end

function test_mst(io::IO)
    #cc = Classic

    city_graph = WeightedGraph{String,WeightedEdge}()
    cities = init!(city_graph)
    @test length(cities) == length(city_graph.vertices)
    @test length(cities) == length(city_graph.edges_lists)

    #show(city_graph)

    mst_path = mst(city_graph, 1)

    #print_weighted_path(city_graph, mst_path)
    print(io, city_graph, mst_path)
    @test 14 == length(mst_path)
    #println("Total Weight: $(sum(mst_path))")
end

function test_DijkstraNode(io::IO)
    #cc = Classic

    node = DijkstraNode(13, 0.1)
    @test 13 == node.index
    @test 0.1 == node.distance

    que = PriorityQueue{DijkstraNode,Float64}()
    enqueue!(que, DijkstraNode(17, 1.2), 1.2)
    enqueue!(que, DijkstraNode(19, 2.3), 2.3)
    @test 2 == length(que)

    # the new node is considered duplicated because
    # its index is 17 even its distance is different
    # this is implemented by isequal
    @test_throws ArgumentError enqueue!(que, DijkstraNode(17, 3.3), 1.1)

    @test 2 == length(que)
    node = dequeue!(que)
    @test 1.2 == node.distance
    node = dequeue!(que)
    @test 2.3 == node.distance
end

function test_dijkstra(io::IO)
    # Create a default graph
    city_graph = WeightedGraph{String,WeightedEdge}()
    cities = init!(city_graph)
    @test length(cities) == length(city_graph.vertices)
    @test length(cities) == length(city_graph.edges_lists)

    # Just a reminder
    # struct DijkstraResult
    #     distances::Vector{Float64}
    #     path_db::Dict{Int64,WeightedEdge}
    # end
    dijkstra_result = dijkstra(city_graph, "Los Angeles")
    @test DijkstraResult == typeof(dijkstra_result)
    println(io, dijkstra_result)
    #@show dijkstra_result

    #ditance_db = Dict{V,Float64}()
    distance_db = distances_to_distance_db(city_graph, dijkstra_result.distances)
    println(io, "Distance from Los Angeles to:")
    println(io, distance_db)
    #@show distance_db

    # path::Vector{E}
    #path = path_db_to_path(
    path = shortest_path(
            dijkstra_result.path_db,
            index_of(city_graph, "Los Angeles"),
            index_of(city_graph, "Boston")
        )
    #@show path

    #print_weighted_path(city_graph, path)
    print(io, stdout, city_graph, path)
    println(io, "Total Weight: $(sum(path))")
    # TODO show better

    @test nothing != dijkstra_result
end

io = devnull
#io = stdout

@testset "Unweighted Graph" begin
    test_SimpleEdge(io)
    test_UnweightedGraph(io)
    test_bfs(io)
end

@testset "Weighted Graph" begin
    test_WeightedEdge(io)
    test_WeightedGraph(io)
    test_mst(io)
end

@testset "Dijkstra" begin
    test_DijkstraNode(io)
    test_dijkstra(io)
end

nothing
# print("""
# # --------------- #
# # Test completed! #
# # --------------- #
# """)


text:
mutable struct WrappedAPI
    source
    structs
    funcs
    consts
    enums
    misc
    bags
    extended_vk_constructors
end

# not efficient
vars(w_api) = OrderedDict([k => v for field ∈ getproperty.(Ref(w_api), [:structs, :funcs, :consts, :enums]) for (k, v) ∈ field])

Base.show(io::IO, w_api::WrappedAPI) = print(io, "Wrapped API with $(length(w_api.structs)) structs, $(length(w_api.funcs)) functions, $(length(w_api.consts)) consts and $(length(w_api.enums)) enums wrapped from $(w_api.source)")

abstract type Wrapper end

default_spacing(::SDefinition) = "\n"^2
default_spacing(decl::FDefinition) = decl.short && length(generate(decl)) < 150 ? "\n" : "\n"^2
default_spacing(::CDefinition) = "\n"
default_spacing(::EDefinition) = "\n"

unsolved_dependencies(sdef, decls) = collect(values(decls)[findall(keys(decl) .== setdiff(type_dependencies(sdef), keys(decls)))])
are_dependencies_solved(def, decls) = isempty(unsolved_dependencies(def, decls))

function write_api!(io::IO, def::Declaration; spacing)
    write(io, generate(def) * spacing(def))
end

exports(symbols) = "export $(join_args(symbols))"

pre_wrap_code = """
abstract type VulkanStruct end
abstract type Handle <: VulkanStruct end
abstract type ReturnedOnly <: VulkanStruct end
abstract type Bag end
struct BagEmpty <: Bag end
const EmptyBag = BagEmpty()

Base.cconvert(T::Type, x::VulkanStruct) = x
Base.cconvert(T::Type{<: Ptr}, x::AbstractArray{<: VulkanStruct}) = getproperty.(x, :vks)
Base.cconvert(T::Type{<: Ptr}, x::AbstractArray{<: Handle}) = getproperty.(x, :handle)
Base.cconvert(T::Type{<: Ptr}, x::VulkanStruct) = Ref(x.vks)
Base.cconvert(T::Type{<: Ptr}, x::Handle) = x
Base.unsafe_convert(T::Type, x::VulkanStruct) = x.vks
Base.unsafe_convert(T::Type{Ptr{Nothing}}, x::Handle) = x.handle

Base.broadcastable(x::VulkanStruct) = Ref(x) # indicate that VulkanStructs behave as scalars for broadcasting

"""

"""Write a wrapped API to files in dest_dir.

Spacing options can be controlled by providing the corresponding argument with a function that takes a Declaration type as argument.
"""
function Base.write(w_api::WrappedAPI, destfile; spacing=default_spacing)
    decls = OrderedDict((vcat(w_api.consts, w_api.enums, w_api.structs, w_api.bags)...)...)
    decls_order = resolve_dependencies(decls)
    check_dependencies(decls, decls_order)
    open(destfile, "w+") do io
        write(io, pre_wrap_code)
    end
    for fdef ∈ values(w_api.extended_vk_constructors)
        open(destfile, "a+") do io
            write_api!(io, fdef; spacing)
        end
    end
    for decl ∈ getindex.(Ref(decls), decls_order)
        open(destfile, "a+") do io
            write_api!(io, decl; spacing)
        end
    end
    open(destfile, "a+") do io
        write(io, "\n\n" * join(w_api.misc, "\n") * "\n\n")
    end
    open(destfile, "a+") do io
        write(io, "\n\n")
        write_api!.(Ref(io), collect(values(w_api.funcs)); spacing)
        write(io, "\n\n" * exports(unique(vcat((map(x -> x.name, filter(x -> x.name ∉ ["Base.convert", "Base", extension_types...], vcat(map(collect ∘ values, [decls, w_api.funcs])...)))), collect(Iterators.flatten(map.(Ref(x -> strip.(first.(split.(x, "=")))), getproperty.(values(w_api.enums), :fields))))))))
    end

    format(destfile)
    nothing
end

stype_splice(fdef) = getproperty(vk, Symbol(stypes[fdef.name]))
name_transform(decl::Declaration) = name_transform(decl.name, typeof(decl))

# arguments whose value is always predetermined by the function signature
# they are dropped as argument and replaced wherever necessary

include("wrapping/struct_logic.jl")
include("wrapping/annotations.jl")
include("wrapping/constructor_logic.jl")
include("wrapping/function_logic.jl")

function wrap!(w_api::WrappedAPI)
    sdefs = collect(values(api.structs))
    errors = OrderedDict()
    wrap!(w_api, sdefs, errors)
    wrap!(w_api, values(api.funcs), errors)
    wrap!(w_api, values(api.consts), errors)
    wrap!(w_api, values(api.enums), errors)
    length(errors) == 0 ? @info("API successfully wrapped.") : @warn("API wrapped with $(length(errors)) errors:")
    for (field, msg) ∈ errors
        println("\t\e[31;1;1m$field: $msg\e[m")
    end
    w_api
end

function wrap!(w_api, objects, errors)
    for obj ∈ objects
        try
            wrap!(w_api, obj)
        catch e
            msg = hasproperty(e, :msg) ? e.msg : "$(typeof(e))"
            errors[obj.name] = msg
            println("\e[31;1;1m$(obj.name): $msg\e[m")
            rethrow(e)
            typeof(e) ∈ [ErrorException, AssertionError] ? continue : rethrow(e)
        end
    end
end

function wrap!(w_api, edef::EDefinition)

    new_edef = EDefinition(remove_vk_prefix(edef.name), remove_vk_prefix.(edef.fields), edef.with_begin_block, isnothing(edef.type) ? nothing : remove_vk_prefix(edef.type), edef.enum_macro)
    w_api.enums[new_edef.name] = new_edef
    w_api.funcs["convert_$(edef.name)"] = FDefinition("Base.convert(T::Type{$(new_edef.name)}, e::$(edef.name)) = T(UInt(e))")
    w_api.funcs["convert_$(new_edef.name)"] = FDefinition("Base.convert(T::Type{$(edef.name)}, e::$(new_edef.name)) = T(UInt(e))")
end

function wrap!(w_api, cdef::CDefinition)
    new_cdef = CDefinition(remove_vk_prefix(cdef.name), is_literal(cdef.value) ? cdef.value : remove_vk_prefix(cdef.value))
    w_api.consts[new_cdef.name] = new_cdef
end

function wrap(library_api::API)
    global api = library_api
    w_api = WrappedAPI(api.source, OrderedDict{String, SDefinition}(), OrderedDict{String, FDefinition}(), OrderedDict{String, CDefinition}(), OrderedDict{String, EDefinition}(), String[], OrderedDict{String, SDefinition}(), OrderedDict{String, FDefinition}())
    wrap!(w_api)
end


text:
# see also the macro versions in datasets.jl

# -------------------------------------------------------
# To add a new dataset assuming it has a header and is, at path
# `data/newdataset.csv`
#
# 1. start by loading it with CSV:
#
#   fpath = joinpath("datadir", "newdataset.csv")
#   data = CSV.read(fpath, copycols=true,
#                   categorical=true)
#
# 2. load it with DelimitedFiles and Tables
#
#   data_raw, data_header = readdlm(fpath, ',', header=true)
#   data_table = Tables.table(data_raw; header=Symbol.(vec(data_header)))
#
# 3. retrieve the conversions:
#
#   for (n, st) in zip(names(data), scitype_union.(eachcol(data)))
#     println(":$n=>$st,")
#   end
#
# 4. copy and paste the result in a coerce
#
#  data_table = coerce(data_table, ...)
#
# -------------------------------------------------------

const DATA_DIR = joinpath(MODULE_DIR, "..", "data")

const COERCE_BOSTON = (:Chas => Count,)

const COERCE_REDUCED_AMES = (
    :target       => Continuous,
    :OverallQual  => OrderedFactor,
    :GrLivArea    => Continuous,
    :Neighborhood => Multiclass,
    :x1stFlrSF    => Continuous,
    :TotalBsmtSF  => Continuous,
    :BsmtFinSF1   => Continuous,
    :LotArea      => Continuous,
    :GarageCars   => Count,
    :MSSubClass   => Multiclass,
    :GarageArea   => Continuous,
    :YearRemodAdd => Count,
    :YearBuilt    => Count)

const COERCE_AMES = (
    :Id             => Count,
    :MSSubClass     => Multiclass,
    :MSZoning       => Multiclass,
    :LotFrontage    => Continuous,
    :LotArea        => Continuous,
    :Street         => Multiclass,
    :LotShape       => Multiclass,
    :LandContour    => Multiclass,
    :LotConfig      => Multiclass,
    :LandSlope      => OrderedFactor,
    :Neighborhood   => Multiclass,
    :Condition1     => Multiclass,
    :Condition2     => Multiclass,
    :BldgType       => Multiclass,
    :HouseStyle     => Multiclass,
    :OverallQual    => OrderedFactor,
    :OverallCond    => OrderedFactor,
    :YearBuilt      => Count,
    :YearRemodAdd   => Count,
    :RoofStyle      => Multiclass,
    :RoofMatl       => Multiclass,
    :Exterior1st    => Multiclass,
    :Exterior2nd    => Multiclass,
    :MasVnrType     => Multiclass,
    :MasVnrArea     => Continuous,
    :ExterQual      => OrderedFactor,
    :ExterCond      => OrderedFactor,
    :Foundation     => Multiclass,
    :BsmtQual       => OrderedFactor,
    :BsmtCond       => OrderedFactor,
    :BsmtExposure   => OrderedFactor,
    :BsmtFinType1   => Multiclass,
    :BsmtFinSF1     => Continuous,
    :BsmtFinType2   => Multiclass,
    :BsmtFinSF2     => Continuous,
    :BsmtUnfSF      => Continuous,
    :TotalBsmtSF    => Continuous,
    :Heating        => Multiclass,
    :HeatingQC      => OrderedFactor,
    :CentralAir     => Multiclass,
    :Electrical     => Multiclass,
    :x1stFlrSF      => Continuous,
    :x2ndFlrSF      => Continuous,
    :LowQualFinSF   => Continuous,
    :GrLivArea      => Continuous,
    :BsmtFullBath   => Count,
    :BsmtHalfBath   => Count,
    :FullBath       => Count,
    :HalfBath       => Count,
    :BedroomAbvGr   => Count,
    :KitchenAbvGr   => Count,
    :KitchenQual    => OrderedFactor,
    :TotRmsAbvGrd   => Count,
    :Functional     => OrderedFactor,
    :Fireplaces     => Count,
    :FireplaceQu    => OrderedFactor,
    :GarageType     => Multiclass,
    :GarageYrBlt    => Count,
    :GarageFinish   => Multiclass,
    :GarageCars     => Count,
    :GarageArea     => Continuous,
    :GarageQual     => OrderedFactor,
    :GarageCond     => OrderedFactor,
    :PavedDrive     => Multiclass,
    :WoodDeckSF     => Continuous,
    :OpenPorchSF    => Continuous,
    :EnclosedPorch  => Continuous,
    :x3SsnPorch     => Continuous,
    :ScreenPorch    => Continuous,
    :PoolArea       => Continuous,
    :MiscVal        => Continuous,
    :MoSold         => Multiclass,
    :YrSold         => Count,
    :SaleType       => Multiclass,
    :SaleCondition  => Multiclass,
    :target         => Continuous)

const COERCE_IRIS = (
    :sepal_length => Continuous,
    :sepal_width  => Continuous,
    :petal_length => Continuous,
    :petal_width  => Continuous,
    :target       => Multiclass)

const COERCE_CRABS = (
    :sp    => Multiclass,
    :sex   => Multiclass,
    :index => Count,
    :FL    => Continuous,
    :RW    => Continuous,
    :CL    => Continuous,
    :CW    => Continuous,
    :BD    => Continuous)

typeof(COERCE_CRABS)

const COERCE_SMARKET = (
    :Year=>Continuous,
    :Lag1=>Continuous,
    :Lag2=>Continuous,
    :Lag3=>Continuous,
    :Lag4=>Continuous,
    :Lag5=>Continuous,
    :Volume=>Continuous,
    :Today=>Continuous,
    :Direction=>Multiclass{2})

"""
load_dataset(fpath, coercions)

Load one of standard dataset like Boston etc assuming the file is a
comma separated file with a header.

"""
function load_dataset(fname::String, coercions::Tuple)
    fpath = joinpath(DATA_DIR, fname)
    data_raw, data_header = readdlm(fpath, ',', header=true)
    data_table = Tables.table(data_raw; header=Symbol.(vec(data_header)))
    return coerce(data_table, coercions...; tight=true)
end

load_boston()       = load_dataset("Boston.csv", COERCE_BOSTON)
load_reduced_ames() = load_dataset("reduced_ames.csv", COERCE_REDUCED_AMES)
function load_ames()
    data = load_dataset("ames.csv", COERCE_AMES)
    levels!(data.LandSlope, ["Gtl", "Mod", "Sev"])
    levels!(data.ExterQual, ["Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.ExterCond, ["Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.HeatingQC, ["Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.KitchenQual, ["Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.BsmtQual, ["_NA", "Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.BsmtCond, ["_NA", "Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.BsmtExposure, ["_NA", "No", "Mn", "Av", "Gd"])
    levels!(data.FireplaceQu, ["None", "Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.GarageQual, ["_NA", "Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.GarageCond, ["_NA", "Po", "Fa", "TA", "Gd", "Ex"])
    levels!(data.Functional, ["Typ", "Min1", "Min2", "Mod", "Maj1", "Maj2",
                              "Sev", "Sal"])
    return data
end
load_iris()         = load_dataset("iris.csv", COERCE_IRIS)
load_crabs()        = load_dataset("crabs.csv", COERCE_CRABS)
function load_smarket()
    data1 = load_dataset("smarket.csv", COERCE_SMARKET)
    return merge(data1, (Year=Dates.Date.(data1.Year),))
end

"""Load a well-known public regression dataset with `Continuous` features."""
macro load_boston()
    quote
        y, X = unpack(load_boston(), ==(:MedV), x->x != :Chas)
        (X, y)
    end
end

"""Load a reduced version of the well-known Ames Housing task"""
macro load_reduced_ames()
    quote
        y, X = unpack(load_reduced_ames(), ==(:target), x-> true)
        (X, y)
    end
end

"""Load the full version of the well-known Ames Housing task."""
macro load_ames()
    quote
        y, X = unpack(load_ames(), ==(:target), x->x != :Id)
        (X, y)
    end
end

"""Load a well-known public classification task with nominal features."""
macro load_iris()
    quote
        y, X = unpack(load_iris(), ==(:target), x-> true)
        (X, y)
    end
end

"""Load a well-known crab classification dataset with nominal features."""
macro load_crabs()
    quote
        y, X = unpack(load_crabs(), ==(:sp), x-> !(x in [:sex, :index]))
        (X, y)
    end
end

""" Load S&P Stock Market dataset, as used in (An Introduction to
Statistical Learning with applications in
R)[https://rdrr.io/cran/ISLR/man/Smarket.html](https://rdrr.io/cran/ISLR/man/Smarket.html),
by Witten et al (2013), Springer-Verlag, New York."""
macro load_smarket()
    quote
        y, X = unpack(load_smarket(), ==(:Direction), x-> true)
        (X, y)
    end
end