text:
using FastGaussQuadrature: gausslegendre

"""
	cache = POOrthogonalCollocationCache(Ntst::Int, m::Int, Ty = Float64)

Structure to hold the cache for the collocation method.

# Arguments
- `Ntst` number of time steps
- `m` degree of the polynomials
- `Ty` type of the time variable
"""
struct POOrthogonalCollocationCache{T}
	Ntst::Int
	degree::Int
    lagrange_vals::Matrix{T}
    lagrange_driv::Matrix{T}
	gauss_nodes::Vector{T}
    gauss_weight::Vector{T}
	# TODO how do we ensure that eltype(mesh) is what we expect?
	mesh::Vector{T} 		# τs, we need a vector here for mesh adaptation
	mesh_coll::LinRange{T} 	# σs
	full_mesh::Vector{T}
end

function POOrthogonalCollocationCache(Ntst::Int, m::Int, Ty = Float64)
	τs = LinRange{Ty}( 0, 1, Ntst + 1) |> collect
	σs = LinRange{Ty}(-1, 1, m + 1)
	L, ∂L = getL(σs)
	zg, wg = gausslegendre(m)
	prob = POOrthogonalCollocationCache(Ntst, m, L, ∂L, zg, wg, τs, σs, zeros(Ty, 1 + m * Ntst))
	# put the mesh where we removed redundant timing
	prob.full_mesh .= getTimes(prob)
	return prob
end

@inline Base.eltype(pb::POOrthogonalCollocationCache) = eltype(pb.lagrange_vals)
@inline Base.size(pb::POOrthogonalCollocationCache) = (pb.degree, pb.Ntst)
@inline getLs(pb::POOrthogonalCollocationCache) = (pb.lagrange_vals, pb.lagrange_driv)
@inline getMesh(pb::POOrthogonalCollocationCache) = pb.mesh
@inline getMeshColl(pb::POOrthogonalCollocationCache) = pb.mesh_coll
getMaxTimeStep(pb::POOrthogonalCollocationCache) = maximum(diff(getMesh(pb)))
τj(σ, τs, j) = τs[j] + (1 + σ)/2 * (τs[j+1] - τs[j])
# get the sigma corresponding to τ in the interval (𝜏s[j], 𝜏s[j+1])
σj(τ, τs, j) = -(2*τ - τs[j] - τs[j + 1])/(-τs[j + 1] + τs[j])

# code from Jacobi.lagrange
function lagrange(i, x, z)
    nz = length(z)
    l = one(z[1])
	for k in 1:(i-1)
        l = l * (x - z[k]) / (z[i] - z[k])
    end
    for k in (i+1):nz
        l = l * (x - z[k]) / (z[i] - z[k])
    end
    return l
end

dlagrange(i, x, z) = ForwardDiff.derivative(x -> lagrange(i, x, z), x)

# accept a range, ie σs = LinRange(-1, 1, m + 1)
function getL(σs::AbstractVector)
	m = length(σs) - 1
	zs, = gausslegendre(m)
	L = zeros(m, m+1); ∂L = zeros(m, m + 1)
	for j in 1:m+1
		for i in 1:m
			L[i, j]  =  lagrange(j, zs[i], σs)
			∂L[i, j] = dlagrange(j, zs[i], σs)
		end
	end
	return (;L, ∂L)
end

"""
$(SIGNATURES)

Return the times at which is evaluated the collocation problem.
"""
function getTimes(pb::POOrthogonalCollocationCache)
	m, Ntst = size(pb)
	Ty = eltype(pb)
	ts = Matrix{Ty}(undef, Ntst, m + 1)
	tsvec = Ty[0]
	τs = pb.mesh
	σs = pb.mesh_coll
	for j in 1:Ntst
		for l in 1:m+1
			ts[j, l] = τj(σs[l], τs, j)
			l>1 && push!(tsvec, τj(σs[l], τs, j))
		end
	end
	return vec(tsvec)
end

function updateMesh!(pb::POOrthogonalCollocationCache, mesh)
	pb.mesh .= mesh
	pb.full_mesh .= getTimes(pb)
end
####################################################################################################

"""
	pb = PeriodicOrbitOCollProblem(kwargs...)

This composite type implements an orthogonal collocation (at Gauss points) method of piecewise polynomials to locate periodic orbits. More details (maths, notations, linear systems) can be found [here](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/periodicOrbitCollocation/).

## Arguments
- `F` vector field specified as a function of two arguments `F(x,p)`
- `J` is the jacobian of `F` at `(x, p)`. It can assume three forms.
    1. Either `J` is a function and `J(x, p)` returns a `::AbstractMatrix`. In this case, the default arguments of `contParams::ContinuationPar` will make `continuation` work.
    2. Or `J` is a function and `J(x, p)` returns a function taking one argument `dx` and returning `dr` of the same type as `dx`. In our notation, `dr = J * dx`. In this case, the default parameters of `contParams::ContinuationPar` will not work and you have to use a Matrix Free linear solver, for example `GMRESIterativeSolvers`,
    3. Or `J` is a function and `J(x, p)` returns a variable `j` which can assume any type. Then, you must implement a linear solver `ls` as a composite type, subtype of `AbstractLinearSolver` which is called like `ls(j, rhs)` and which returns the solution of the jacobian linear system. See for example `examples/SH2d-fronts-cuda.jl`. This linear solver is passed to `NewtonPar(linsolver = ls)` which itself passed to `ContinuationPar`. Similarly, you have to implement an eigensolver `eig` as a composite type, subtype of `AbstractEigenSolver`.
- `ϕ::AbstractVector` used to set a section for the phase constraint equation
- `xπ::AbstractVector` used in the section for the phase constraint equation
- `N::Int` dimension of the state space
- `coll_cache::POOrthogonalCollocationCache` cache for collocation. See docs of `POOrthogonalCollocationCache` .

## Methods

Here are some useful methods you can apply to `pb`

- `length(pb)` gives the total number of unknowns
- `size(pb)` returns the triplet `(N, m, Ntst)`
- `getMesh(pb)` returns the mesh `0 = τ0 < ... < τNtst+1 = 1`. This is useful because this mesh is born to vary by automatic mesh adaptation
- `getMeshColl(pb)` returns the (static) mesh `0 = σ0 < ... < σm+1 = 1`
- `getTimes(pb)` returns the vector of times (length `1 + m * Ntst`) at the which the collocation is applied.
- `generateSolution(pb, orbit, period)` generate a guess from a function `t -> orbit(t)` which approximates the periodic orbit.
- `POOcollSolution(pb, x)` return a function interpolating the solution `x` using a piecewise polynomials function

# Orbit guess
You will see below that you can evaluate the residual of the functional (and other things) by calling `pb(orbitguess, p)` on an orbit guess `orbitguess`. Note that `orbitguess` must be of size 1 + N * (1 + m * Ntst) where N is the number of unknowns in the state space and `orbitguess[end]` is an estimate of the period ``T`` of the limit cycle.

Note that you can generate this guess from a function using `generateSolution`.

# Functional
 A functional, hereby called `G`, encodes this problem. The following methods are available

- `pb(orbitguess, p)` evaluates the functional G on `orbitguess`

"""
@with_kw_noshow struct PeriodicOrbitOCollProblem{TF, TJ, vectype, Tmass, Tcache <: POOrthogonalCollocationCache} <: AbstractPeriodicOrbitProblem
	# Function F(x, par)
	F::TF = nothing

	# Jacobian of F w.r.t. x
	J::TJ = nothing

	# variables to define a Section for the phase constraint equation
	ϕ::vectype = nothing
	xπ::vectype = nothing

	# dimension of the problem in case of an AbstractVector
	N::Int = 0

	# whether the time discretisation is adaptive
	adaptmesh::Bool = false

	# whether the problem is nonautonomous
	isautonomous::Bool = true

	# mass matrix
	massmatrix::Tmass = nothing

	# collocation cache
	coll_cache::Tcache = nothing
end

# trivial constructor
function PeriodicOrbitOCollProblem(Ntst, m, N = 0)
	cache = POOrthogonalCollocationCache(Ntst, m)
	PeriodicOrbitOCollProblem(; N = N, coll_cache = cache)
end

# TODO rename this in num_mesh? or meshSize
@inline getMeshSize(pb::PeriodicOrbitOCollProblem) = pb.coll_cache.Ntst

@inline length(pb::PeriodicOrbitOCollProblem) = ( (n, m, Ntst) = size(pb); return n * (1 + m * Ntst) )
# the size is (n, m, Ntst)
@inline Base.size(pb::PeriodicOrbitOCollProblem) = (pb.N, size(pb.coll_cache)...)
@inline Base.eltype(pb::PeriodicOrbitOCollProblem) = eltype(pb.coll_cache)
getLs(pb::PeriodicOrbitOCollProblem) = getLs(pb.coll_cache)

# these functions extract the time slices components
getTimeSlices(x::AbstractVector, N, degree, Ntst) = @views reshape(x[1:end-1], N, (degree) * Ntst + 1)
# array of size Ntst ⋅ (m+1) ⋅ n
getTimeSlices(pb::PeriodicOrbitOCollProblem, x) = getTimeSlices(x, size(pb)...)
getTimes(pb::PeriodicOrbitOCollProblem) = getTimes(pb.coll_cache)
getMesh(pb::PeriodicOrbitOCollProblem) = getMesh(pb.coll_cache)
getMeshColl(pb::PeriodicOrbitOCollProblem) = getMeshColl(pb.coll_cache)
getMaxTimeStep(pb::PeriodicOrbitOCollProblem) = getMaxTimeStep(pb.coll_cache)
updateMesh!(pb::PeriodicOrbitOCollProblem, mesh) = updateMesh!(pb.coll_cache, mesh)

function Base.show(io::IO, pb::PeriodicOrbitOCollProblem)
	N, m, Ntst = size(pb)
	println(io, "┌─ Collocation problem for periodic orbits")
	println(io, "├─ type            : Vector{", eltype(pb), "}")
	println(io, "├─ time slices     : ", Ntst)
	println(io, "├─ degree          : ", m)
	println(io, "├─ space dimension : ", pb.N)
	println(io, "└─ # unknowns      : ", pb.N * (1 + m * Ntst))
end

"""
$(SIGNATURES)

This function generates an initial guess for the solution of the problem `pb` based on the orbit `t -> orbit(t)` for t ∈ [0,1] and the period `period`.
"""
function generateSolution(pb::PeriodicOrbitOCollProblem, orbit, period)
	n, _m, Ntst = size(pb)
	ts = getTimes(pb)
	Nt = length(ts)
	ci = zeros(eltype(pb), n, Nt)
	for (l, t) in pairs(ts)
		ci[:, l] .= orbit(t * period)
	end
	return vcat(vec(ci), period)
end

@views function (prob::PeriodicOrbitOCollProblem)(u::AbstractVector, pars)
	uc = getTimeSlices(prob, u)
	T = getPeriod(prob, u, pars)
	result = zero(u)
	resultc = getTimeSlices(prob, result)
	functionalColl!(prob, resultc, uc, T, getLs(prob.coll_cache), pars)
	# add  the phase condition
	result[end] = dot(u[1:end-1], prob.ϕ) - dot(prob.xπ, prob.ϕ)
	return result
end

# function or collocation problem
@views function functionalColl!(pb, out, u, period, (L, ∂L), pars)
	Ty = eltype(u)
	n, ntimes = size(u)
	m = pb.coll_cache.degree
	Ntst = pb.coll_cache.Ntst
	# on veut faire des slices a t fixes, donc des gj[:, j], c'est le plus rapide
	gj = zeros(Ty, n, m)
	∂gj = zeros(Ty, n, m)
	uj = zeros(Ty, n, m+1)
	mesh = getMesh(pb)
	rg = 1:m+1
	for j in 1:Ntst
		uj .= u[:, rg]
		gj 	.=  uj * L'
		∂gj .=  uj * ∂L'
		# mul!(gj, uj, L')
		# mul!(∂gj, uj, ∂L')
		# compute the collocation residual
		for l in 1:m
			out[:, rg[l]] .= ∂gj[:, l] .- (period * (mesh[j+1]-mesh[j]) / (2)) .* pb.F(gj[:, l], pars)
		end
		rg = rg .+ m
	end
	# add the periodicity condition
	out[:, end] .= u[:, end] .- u[:, 1]
end

"""
$(SIGNATURES)

Compute the full periodic orbit associated to `x`. Mainly for plotting purposes.
"""
@views function getPeriodicOrbit(prob::PeriodicOrbitOCollProblem, u::AbstractVector, p)
	T = getPeriod(prob, u, p)
	ts = getTimes(prob)
	uc = getTimeSlices(prob, u)
	n, m, Ntst = size(prob)
	return SolPeriodicOrbit(t = ts .* T, u = uc)
end

# function needed for automatic Branch switching from Hopf bifurcation point
function reMake(prob::PeriodicOrbitOCollProblem, F, dF, par, hopfpt, ζr::AbstractVector, orbitguess_a, period; orbit = t->t)
	M = length(orbitguess_a)
	N = length(ζr)

	n, m, Ntst = size(prob)
	nunknows = N * (1 + m*Ntst)

	# update the problem
	probPO = setproperties(prob, N = N, F = F, J = dF, ϕ = zeros(nunknows), xπ = zeros(nunknows))

	probPO.ϕ[1:N] .= ζr
	probPO.xπ[1:N] .= hopfpt.x0

	# append period at the end of the initial guess
	orbitguess = generateSolution(probPO, t -> orbit(2pi*t/period), period)

	return probPO, orbitguess
end

####################################################################################################
const DocStrjacobianPOColl = """
- `jacobianPO` Specify the choice of the linear algorithm, which must belong to `(:autodiffDense, :FiniteDifferences)`. This is used to select a way of inverting the jacobian dG
    - For `:autodiffDense`. Same as for `:autodiffMF` but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using `options`.
    - For `:FiniteDifferencesDense`, same as for `:autodiffDense` but we use Finite Differences to compute the jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.
"""

"""
$(SIGNATURES)

This is the Newton Solver for computing a periodic orbit using (Standard / Poincaré) Shooting method.
Note that the linear solver has to be apropriately set up in `options`.

# Arguments

Similar to [`newton`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref). These two problems have specific options to be tuned, we refer to their link for more information and to the tutorials.

- `prob` a problem of type `<: AbstractShootingProblem` encoding the shooting functional G.
- `orbitguess` a guess for the periodic orbit. See [`ShootingProblem`](@ref) and See [`PoincareShootingProblem`](@ref) for information regarding the shape of `orbitguess`.
- `par` parameters to be passed to the functional
- `options` same as for the regular [`newton`](@ref) method.

# Optional argument
$DocStrjacobianPOColl

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(prob::PeriodicOrbitOCollProblem, orbitguess, par, options::NewtonPar;
		jacobianPO = :autodiffDense, kwargs...)
	@assert jacobianPO in
			(:autodiffDense, ) "This jacobian is not defined. Please chose another one."

	if jacobianPO == :autodiffDense
		jac = (x, p) -> ForwardDiff.jacobian(z -> prob(z, p), x)
	end

	return newton(prob, jac, orbitguess, par, options; kwargs...)
end

"""
$(SIGNATURES)

This is the continuation method for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

# Arguments

Similar to [`continuation`](@ref) except that `prob` is either a [`ShootingProblem`](@ref) or a [`PoincareShootingProblem`](@ref). By default, it prints the period of the periodic orbit.

# Optional argument
- `δ = 1e-8` used for finite differences
- `jacobianPO` Specify the choice of the linear algorithm, which must belong to `[:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferences]`. This is used to select a way of inverting the jacobian dG
- `updateSectionEveryStep = 0` updates the section every `updateSectionEveryStep` step during continuation

# Choices for `jacobianPO`
$DocStrjacobianPOColl
"""
function continuation(
	prob::PeriodicOrbitOCollProblem,
	orbitguess, par, lens::Lens, _contParams::ContinuationPar,
	_linearAlgo::AbstractBorderedLinearSolver;
	jacobianPO = :autodiffDense,
	updateSectionEveryStep = 0, kwargs...)
	@assert jacobianPO in
			(:autodiffDense,) "This jacobian is not defined. Please chose another one."

	jac = (x, p) -> FloquetWrapper(prob, ForwardDiff.jacobian(z -> prob(z, p), x), x, p)

	# we have to change the Bordered linearsolver to cope with our type FloquetWrapper
	linearAlgo = @set _linearAlgo.solver = FloquetWrapperLS(_linearAlgo.solver)

	options = _contParams.newtonOptions
	contParams = @set _contParams.newtonOptions.linsolver = FloquetWrapperLS(options.linsolver)

	if computeEigenElements(contParams)
		contParams = @set contParams.newtonOptions.eigsolver = FloquetLUColl(contParams.newtonOptions.eigsolver, length(prob), prob.N)
	end

	# change the user provided finalise function by passing prob in its parameters
	_finsol = modifyPOFinalise(prob, kwargs, updateSectionEveryStep)
	_recordsol = modifyPORecord(prob, kwargs, par, lens)
	_plotsol = modifyPOPlot(prob, kwargs)

	branch, u, τ = continuation(
					prob, jac,
					orbitguess, par, lens,
					contParams, linearAlgo;
					kwargs...,
					recordFromSolution = _recordsol,
					finaliseSolution = _finsol,
					plotSolution = _plotsol)
	return setproperties(branch; type = :PeriodicOrbit, functional = prob), u, τ
end

"""
$(SIGNATURES)

Compute the maximum of the periodic orbit associated to `x`.
"""
function getMaximum(prob::PeriodicOrbitOCollProblem, x::AbstractVector, p)
	sol = getPeriodicOrbit(prob, x, p).u
	return maximum(sol)
end

# this function updates the section during the continuation run
@views function updateSection!(prob::PeriodicOrbitOCollProblem, x, par; stride = 0)
	n, m, Ntst = size(prob)
	xc = getTimeSlices(prob, x)
	T = getPeriod(prob, x, par)

	# update the reference point
	prob.xπ .= x[1:end-1]

	# update the normals
	ϕc = @views reshape(prob.ϕ, n, (m) * Ntst + 1)
	for ii in 1: (m) * Ntst + 1
		ϕc[:, ii] .= prob.F(xc[:, ii], par) ./ (m*Ntst)
	end
	return true
end
####################################################################################################
# mesh adaptation routine

# iterated derivatives
∂(f) = x -> ForwardDiff.derivative(f, x)
∂(f, n) = n == 0 ? f : ∂(∂(f), n-1)

struct POOcollSolution{Tpb, Tx}
	pb::Tpb
	x::Tx
end

@views function (sol::POOcollSolution)(t0)
	n,m,Ntst = size(sol.pb)
	xc = getTimeSlices(sol.pb, sol.x)

	T = getPeriod(sol.pb, sol.x, 0)
	t = t0 / T

	mesh = getMesh(sol.pb)
	indτ = searchsortedfirst(mesh, t) - 1
	if indτ <= 0
		return sol.x[1:n]
	elseif indτ > Ntst
		return xc[:, end]
	end
	# println("--> ", t, " belongs to ", (mesh[indτ], mesh[indτ+1])) # waste lots of ressources
	@assert mesh[indτ] <= t <= mesh[indτ+1] "Please open an issue on the website of BifurcationKit.jl"
	σ = σj(t, mesh, indτ)
	# @assert -1 <= σ <= 1 "Strange value of $σ"
	σs = getMeshColl(sol.pb)
	out = zeros(typeof(t), sol.pb.N)
	rg = (1:m+1) .+ (indτ-1)*m
	for l in 1:m+1
		out .+= xc[:, rg[l]] .* lagrange(l, σ, σs)
	end
	out
end


text:
# This file was generated, do not modify it. # hide
rms(y[test], MLJ.predict(tree, rows=test))

ext: jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using DocOpt

function compute(src, dst, N)
  local z, y, x
  (x,y,z) = size(src)
  for n = 1:N
    for k = 1:z
      for j = 1:y
        for i = 1:x
          if (i == 1 || i == x || j == 1 || j == y || k == 1 || k == z)
            dst[i,j,k] = src[i,j,k]
          else
            dst[i,j,k] = (src[i,j,k-1] + src[i,j,k+1] + src[i,j-1,k] + src[i,j+1,k] + src[i-1,j,k] + src[i+1,j,k]) / 6
          end
        end
      end
    end
    tmp=src
    src=dst
    dst=tmp
  end
  return src
end

function initialize(x, y, z)
  src = Array(Float32, x, y, z)
  rand!(src)
  dst = copy(src)
  return src, dst
end

function checksum(src)
  local z, y, x
  (x, y, z) = size(src)
  sum(src) / (x * y * z)
end

function main()
    doc = """laplace-3d.jl

Laplace 6-point 3D stencil.

Usage:
  laplace-3d.jl -h | --help
  laplace-3d.jl [--size=<size>] [--iterations=<iterations>]

Options:
  -h --help                  Show this screen.
  --size=<size>              Specify a 3d array size (<size> x <size> x <size>) [default: 300].
  --iterations=<iterations>  Specify a number of iterations [default: 100].
"""
    arguments = docopt(doc)

    size = parse(Int, arguments["--size"])
    iterations = parse(Int, arguments["--iterations"])

    srand(0)

    function laplace_3d(iterations)
        local src, dst
        println("Run laplace-3d with size ", size, "x", size, "x", size, " for ", iterations, " iterations.")
        (src, dst) = initialize(size, size, size)
        tic()
        compute(src, dst, 0)
        println("SELFPRIMED ", toq())

        tic()
        dst = compute(src, dst, iterations)
        println("SELFTIMED ", toq())
        println("checksum: ", checksum(dst))
    end

    laplace_3d(iterations)

end

main()


text:
module PaddedMatrices

using VectorizationBase, ArrayInterface,
    SLEEFPirates, VectorizedRNG,
    LoopVectorization, LinearAlgebra,
    Random, Base.Threads#, StackPointers#,
    # SpecialFunctions # Perhaps there is a better way to support erf?

using VectorizationBase: align, gep, AbstractStridedPointer, AbstractSIMDVector, vnoaliasstore!, staticm1,
    static_sizeof, lazymul, vmul, vadd, vsub, StridedPointer, gesp, zero_offsets, pause,
    CACHE_COUNT, NUM_CORES, CACHE_INCLUSIVITY, zstridedpointer
using LoopVectorization: maybestaticsize, mᵣ, nᵣ, preserve_buffer, CloseOpen
using ArrayInterface: StaticInt, Zero, One, OptionallyStaticUnitRange, size, strides, offsets, indices,
    static_length, static_first, static_last, axes,
    dense_dims, DenseDims, stride_rank, StrideRank
# using Threads: @spawn
# import ReverseDiffExpressionsBase:
    # RESERVED_INCREMENT_SEED_RESERVED!, ∂getindex,
    # alloc_adjoint, uninitialized, initialized, isinitialized
# import LoopVectorization: isdense

# using Parameters: @unpack

export @StrideArray, @gc_preserve, # @Constant,
    AbstractStrideArray, AbstractStrideVector, AbstractStrideMatrix,
    StrideArray, StrideVector, StrideMatrix,
    PtrArray,# PtrVector, PtrMatrix,
    # ConstantArray, ConstantVector, ConstantMatrix, allocarray,
    matmul!, matmul_serial!, mul!, *ˡ, StaticInt,
    matmul, matmul_serial
# LazyMap,

include("type_declarations.jl")
include("l3_cache_buffer.jl")
include("size_and_strides.jl")
include("adjoints.jl")
include("stridedpointers.jl")
include("indexing.jl")
include("initialization.jl")
include("views.jl")
include("rand.jl")
include("kernels.jl")
include("blas.jl")
include("broadcast.jl")
include("miscellaneous.jl")

# Commented, because I don't want this to be the only doc string.
# """
# To find a mode, define methods for `logdensity` and logdensity_and_gradient!` dispatching on obj, and evaluating at the position `q`.

# logdensity(obj, q, [::StackPointer])
# ∂logdensity!(∇, obj, q, [::StackPointer])

# These must return a value (eg, a logdensity). logdensity_and_gradient! should store the gradient in ∇.
# """
function logdensity end
function ∂logdensity! end

@generated function calc_factors(::Val{nc} = Val{NUM_CORES}()) where {nc}
    t = Expr(:tuple)
    for i ∈ nc:-1:1
        d, r = divrem(nc, i)
        iszero(r) && push!(t.args, (i, d))
    end
    t
end
const CORE_FACTORS = calc_factors()

const BCACHE = Float64[]
"""
Length is one less than `Base.nthreads()`
"""
const MULTASKS = Task[]
_nthreads() = min(NUM_CORES, length(MULTASKS));

function runfunc(t::Task, tid)
    t.sticky = true
    ccall(:jl_set_task_tid, Cvoid, (Any, Cint), t, tid)
    push!(@inbounds(Base.Workqueues[tid+1]), t)
    ccall(:jl_wakeup_thread, Cvoid, (Int16,), tid % Int16)
    t
end
runfunc(func, tid) = runfunc(Task(func), tid)
function runfunc!(ft, tid)
    @inbounds MULTASKS[tid] = runfunc(ft, tid)
    nothing
end

function __init__()
    resize!(BCACHE, BSIZE * BCACHE_COUNT)
    _nt = nthreads() - 1
    _nt > 1 && resize!(MULTASKS, _nt)
    if _nt < NUM_CORES
        msg = string(
            "Your system has $NUM_CORES physical cores, but `PaddedMatrices.jl` only has ",
            "$(_nt > 1 ? "$(_nt) threads" : "1 thread") available ",
            "(it doesn't spawn tasks on `threadid() == 1`). ",
            "You should start Julia with at least $(NUM_CORES + 1) threads.",
            "",
        )
    end
end

# include("precompile.jl")
# _precompile_()

end # module


text:
module CellListMapBenchmarks

import Pkg
using Dates
using CellListMap

include("./namd/simulate.jl")
using .Simulation

export version
export namd10k, namd100k

# function to get the CellListMap version used
get_version()=filter(x-> x.second.name == "CellListMap", Pkg.dependencies()) |> x -> first(x)[2].version

# Package dir
const package_dir = "$(@__DIR__)/.."

# Namd dir
const namd_dir = "/home/leandro/programs/namd/NAMD_2.14_Linux-x86_64-multicore"

# Create benchmark dir
function create_dir(version,name)
    dir="$package_dir/results/$version/$name"
    mkpath(dir)
    println("Results directory for $version/$name created.")
    return dir
end

# Run namd 10k benchmark
function namd10k()
    hostname=gethostname()
    name="namd10k"
    version = get_version()
    dir = create_dir(version,name)
    working_dir = Simulation.working_dir
    params = Simulation.Params(
        x0=Simulation.getcoor(working_dir*"/ne10k_initial.pdb"),
        trajfile=working_dir*"/ne10k_traj.xyz",
        nsteps=5000
    )
    np = Threads.nthreads()
    t_cl = @elapsed Simulation.simulate(params)
    t_namd = @elapsed run(`$namd_dir/namd2 +p$np $working_dir/ne10k.namd`)
    log = open(dir*"/$hostname-$(np)-threads.dat","w")
    println(log,"Number of threads: $np")
    println(log,"CellListMap: $t_cl")
    println(log,"Namd 2.14:   $t_namd")
    close(log)
end

# Run namd 100k benchmark
function namd100k()
    hostname=gethostname()
    name="namd100k"
    version = get_version()
    dir = create_dir(version,name)
    working_dir = Simulation.working_dir
    params = Simulation.Params(
        x0=Simulation.getcoor(working_dir*"/ne100k_initial.pdb"),
        trajfile=working_dir*"/ne100k_traj.xyz",
        nsteps=1000,
        box = Box([ 99.88, 99.88, 99.88 ], 12., lcell=2)
    )
    np = Threads.nthreads()
    t_cl = @elapsed Simulation.simulate(params)
    t_namd = @elapsed run(`$namd_dir/namd2 +p$np $working_dir/ne100k.namd`)
    log = open(dir*"/$hostname-$(np)-threads.dat","w")
    println(log,"Number of threads: $np")
    println(log,"CellListMap: $t_cl")
    println(log,"Namd 2.14:   $t_namd")
    close(log)
end

end


text:
using PartialRejectionSampling
using Random

using RCall;
@rimport spatstat as rspatstat;

# Create the point process
β₀ = 0.1
r = 0.02  # interaction range = 2*radius
b = β₀ / (π * (r/2)^2)

c, w = zeros(2), 1.0
win = PRS.SquareWindow(c, w)

hc = PRS.HardCorePointProcess(b, r, win)

# Time comparison PRS / spatstat
rng = Random.MersenneTwister(123)

@time PRS.generate_sample_prs(rng, hc);
@time rspatstat.rHardcore(beta=hc.β, R=hc.r);

# K-plot with spatstat
R"""
library(spatstat)
sampl_spatstat = rHardcore(beta=$b, R=$r)
K <- Kest(sampl_spatstat)
plot(K)
"""

# K-plot with PRS
sampl_PRS = hcat(PRS.generate_sample_prs(rng, hc)...)

R"""
library(spatstat)

tmp = $(hc.window.c) + $(hc.window.w)
xrange <- c($(hc.window.c[1]), tmp[1])
yrange <- c($(hc.window.c[2]), tmp[2])
window <- owin(xrange, yrange)

sampl <- ppp($(sampl_PRS[1, :]), $(sampl_PRS[2, :]), window)
K <- Kest(sampl)
plot(K)
"""


text:
# function memory_saved(dt,n)

#     # total size and levels
#     tsize = sum([ismissing(x) ? 0 : length(x.value) for x in dt.data[n]])
#     s = Set(dt.data[n][:].values)

#     # for refs size
#     refsize = length(s) < typemax(UInt8) ? 1 : length(s) < typemax(UInt16) ? 2 : 4

#     # approximate size in categorical arrays
#     rsize = sum(length.(collect(s))) + refsize*dt.rows

#     return tsize*.8 > rsize ? true : false
# end

function get_numbytes(typelist,nvar)
    nb = Vector{UInt16}(undef,nvar)
    for i in 1:nvar
        for i in 1:nvar
            if 0 < typelist[i] < 2045
                nb[i] = typelist[i]
            elseif typelist[i] == 32768
                nb[i] = 8 # 8 bytes
            elseif typelist[i] == 65526
                nb[i] = 8 # double
            elseif typelist[i] == 65527
                nb[i] = 4 # float
            elseif typelist[i] == 65528
                nb[i] = 4 # long
            elseif typelist[i] == 65529
                nb[i] = 2 # int
            elseif typelist[i] == 65530
                nb[i] = 1 # byte
            end
        end
        return nb
    end
end

function strtonull(str)

	n = findfirst(x->isequal('\u0',x) || isvalid(Char,x) == false,str)
	if n == nothing
		return str
	end
	if n == 1
		return ""
	end
	return str[1:n-1]
end

function alloc_array(vtype,vfmt,nobs::Int64)

    # create an Array for the relevant type
    if 0 <= vtype < 2045 || vtype == 32768 # string variable
        return Vector{Union{Missing,String}}(undef,nobs)
    elseif vtype == 65526
        if vfmt == "%d" || vfmt[1:3] == "%td"
            return Vector{Union{Missing,Date}}(undef,nobs)
        elseif vfmt[1:3] == "%tc" || vfmt[1:3] == "%tC"
            return Vector{Union{Missing,DateTime}}(undef,nobs)
        else
            return Vector{Union{Missing,Float64}}(undef,nobs)
        end
    elseif vtype == 65527
        if vfmt == "%d" || vfmt[1:3] == "%td"
            return Vector{Union{Missing,Date}}(undef,nobs)
        elseif vfmt[1:3] == "%tc" || vfmt[1:3] == "%tC"
            return Vector{Union{Missing,DateTime}}(undef,nobs)
        else
            return Vector{Union{Missing,Float32}}(undef,nobs)
        end
    elseif vtype == 65528
        if vfmt == "%d" || vfmt[1:3] == "%td"
            return Vector{Union{Missing,Date}}(undef,nobs)
        elseif vfmt[1:3] == "%tc" || vfmt[1:3] == "%tC"
            return Vector{Union{Missing,DateTime}}(undef,nobs)
        else
            return Vector{Union{Missing,Int32}}(undef,nobs)
        end
    elseif vtype == 65529
        if vfmt == "%d" || vfmt[1:3] == "%td"
            return Vector{Union{Missing,Date}}(undef,nobs)
        else
            return Vector{Union{Missing,Int16}}(undef,nobs)
        end
    elseif vtype == 65530
        return Vector{Union{Missing,Int8}}(undef,nobs)
    end

    error(vtype, " is not a valid variable type in Stata.")
end

"""
	read_stata(fn::String; chunks::Int=10, read_labels=false)

converts a stata datafile `fn` to Julia DataFrame. The original data file bigger than 100MB will be read in `chunks` (default = 10)
to save memory. If `read_labels` is set to `true`, it will not convert the data but will extract labels (both variable labels
and value labels) from the stata data file and convert them to Julia [Lables](https://github.com/mwsohn/Labels.jl).
"""
function read_stata(fn::String; chunks::Int=10, read_labels=false)

    fh = open(fn,"r")

    # dta file (<stata_dta><header><release>)
    header = String(read(fh,67))
    if header[2:10] != "stata_dta"
        error("Not a version 13 or 14 data file")
    end

    # data format version
    release = parse(Int16,header[29:31])
	if release == 117 # version 13
		len_varname = 33
		len_format = 49
		len_labelname = 33
		len_varlabel = 81
	elseif release == 118 # version 14
		len_varname = 129
		len_format = 57
		len_labelname = 129
		len_varlabel = 321
	else
		error("Can't convert data format version ",release,".")
	end

    # byte order: LSF or MSF
    byteorder = header[53:55]
    if byteorder == "MSF"
        error("Big-endian data are not supported yet.")
    end

    # number of variables
    skip(fh,3) # <K>
    nvar = read(fh,Int16)

    # number of observations
    skip(fh,7) #</K><N>
    nobs = Int(read(fh,Int32))

    # dataset label length
	if release == 117
		skip(fh,11)
		dslabel_len = read(fh,Int8)
    elseif release == 118
		skip(fh,15)
		dslabel_len = read(fh,Int16)
	end

    # read the label
    if dslabel_len > 0
        dslabel = String(read(fh,dslabel_len))
    end

    # time stamp
    skip(fh,19)
    timestamplen = read(fh,Int8)
    if timestamplen > 0
        timestamp = String(read(fh,timestamplen))
    end

    # map
    skip(fh,26) # </timestamp></header><map>
    statamap =Vector{Int64}(undef,14)
    read!(fh,statamap)

    # variable types
    skip(fh,22) # </map><variable_types>
    typelist = Vector{UInt16}(undef,nvar)
    read!(fh,typelist)
    #print(Int.(typelist))

    # variable names
    skip(fh,27)
    varlist = Vector{Symbol}(undef,nvar)
    for i in 1:nvar
        varlist[i] = Symbol(strtonull(String(read(fh,len_varname))))
    end

    # sort list
    skip(fh,21) # </varnames><sortlist>
    for i in 1:nvar
        srt = read(fh,Int16)
    end

    # formats
    skip(fh,22) # </sortlist><formats> + 2 (2 bytes left over from the previous sequence)
    fmtlist = Vector{String}(undef,nvar)
    for i in 1:nvar
        fmtlist[i] = strtonull(String(read(fh,len_format)))
    end

    # value label names
    skip(fh,29) # </formats><value_label_names>
    valuelabels = Vector{String}(undef,nvar)
    numvlabels = 0
    for i in 1:nvar
        valuelabels[i] = strtonull(String(read(fh,len_labelname)))

        # count the number of value labels
        if length(valuelabels[i]) > 0
            numvlabels += 1
        end
    end

    # variable labels
    skip(fh,37) # </value_label_names><variable_labels>
    varlabels = Vector{String}(undef,nvar)
    for i in 1:nvar
        varlabels[i] = strtonull(String(read(fh,len_varlabel)))
    end

    # characteristics - we will not import them
    # read until we hit '</characteristics>'
    skip(fh,35) # </variable_labels><characteristics>
    while (true)
        readuntil(fh,'<')
        if String(copy(read(fh,5))) == "data>"
            break
        end
    end

    # nubmer of bytes for each variable
    numbytes = get_numbytes(typelist,nvar)

    # total length of each observation
    rlen = sum(numbytes)

    # number of bytes to skip in IOBuffer
    numskip = zeros(Int,nvar)
    numskip[1] = 0
    for i in 2:length(numbytes)
        numskip[i] = numbytes[i-1] + numskip[i-1]
    end

	# save the start position of the data section
	data_pos = position(fh)

	# skip the data section for now
	skip(fh,rlen*nobs)

    # if there is strLs, read them now
    skip(fh,7) # </data>
    tst = String(read(fh,7))
    if tst == "<strls>"
        # read strLs
        #
        # strL. Stata 13 introduced long strings up to 2 billon characters. strLs are
        # separated by "GSO".
        # (v,o): Position in the data.frame.
        # t:     129/130 defines whether or not the strL is stored with a binary 0.
        # len:   length of the strL.
        # strl:  long string.
        strls = OrderedDict()
        t = OrderedDict()
        while (String(read(fh,3)) == "GSO")
            v = read(fh,Int32)
            if release == 117
                o = read(fh,Int32)
            elseif release == 118
                o = read(fh,Int64)
            end
            t[(v,o)] = read(fh,UInt8)
            len = read(fh,UInt32)
            strls[(v,o)] = String(read(fh,len))
        end
    end

	tst = String(read(fh,5))
    if tst == "trls>"
        skip(fh,14)
    else
        error("Wrong position")
    end

	if read_labels == true

		# read value labels
    	# loop through all value labels
    	# define a Dict first

    	value_labels = Dict()

    	# numvlabels is defined above at the header section
    	for i in 1:numvlabels

			skipstr = String(read(fh,5))
			if skipstr != "<lbl>"
				break
			end
        	len = read(fh,Int32)
        	labname = Symbol(strtonull(String(read(fh,len_labelname))))

        	skip(fh,3) # padding
	        numvalues = read(fh,Int32) # number of entries
	        txtlen = read(fh,Int32) # length of value label text
	        value_labels[labname] = Dict()

	        #
	        offset = Vector{Int32}(undef,numvalues)
			read!(fh,offset) # offset
			values = Vector{Int32}(undef,numvalues)
	        read!(fh,values) # values
	        valtext = String(read(fh,txtlen)) # text table

	        for k in 1:numvalues
	            if k == numvalues
	                offset_end = txtlen
	            else
	                offset_end = offset[k+1]
	            end
	            if values[k] != ""
	                value_labels[labname][values[k]] = strtonull(valtext[offset[k]+1:offset_end])
	            end
	        end
	        skip(fh,6) # </lbl>
	    end

		variable_dict = Dict()
		lblname_dict = Dict()
	    for i in 1:nvar
	        variable_dict[varlist[i]] = varlabels[i]
	        lblname_dict[varlist[i]] = Symbol(valuelabels[i])
	    end

		return Label(variable_dict,value_labels,lblname_dict)
	end

	# read data now
	seek(fh,data_pos)

	# if the data size < 100MB, then
	# slurp the entire data section into memory
	# otherwise, we will read the data by smaller batches
	io = IOBuffer()
	if rlen*nobs < 100_000_000
		write(io,read(fh,rlen*nobs))
		seek(io,0)
		rdf = _read_dta(io,release,rlen,nobs,nvar,varlist,typelist,fmtlist,numskip,strls)
	else
		len = max(100000,ceil(Int,nobs/chunks))
		totlen = nobs
		rdf = DataFrame()
		for nread in 1:ceil(Int,nobs/len)
			if len < totlen
				totlen -= len
			else
				len = totlen
			end
			seek(io,0)
			write(io,read(fh,rlen*len))
			seek(io,0)

			rdf = vcat(rdf,_read_dta(io,release,rlen,len,nvar,varlist,typelist,fmtlist,numskip,strls))
		end
	end

	# close the file
	close(fh)

	# free memory
	# GC.gc()

	return rdf
end

function read_labels(fn::String)
	return read_stata(fn,read_labels=true)
end

function _read_dta(io, release, rlen, len, nvar,varlist,typelist,fmtlist,numskip,strls)

	df = DataFrame()

	dataitemf32::Float32 = 0.
	dataitemf64::Float64 = 0.
	dataitemi8::Int8 = 0
	dataitemi16::Int16 = 0
	dataitemi32::Int32 = 0
	v::Int32 = 0
	o::Int64 = 0
	z = zeros(UInt8,8)

	# interate over the number of variables
	for j in 1:nvar

		df[!,varlist[j]] = alloc_array(typelist[j],fmtlist[j],len)

		# len == number of observation in the batch
		for i in 1:len

			seek(io,numskip[j] + (i-1)*rlen)

			if 0 <= typelist[j] < 2045
				df[i,j] = strtonull(String(read(io,typelist[j])))
				# if empty string, return missing
				if df[i,j] == ""
					df[i,j] = missing
				end
			elseif typelist[j] == 32768 # long string
				if release == 117
					v = read(io,Int32)
					o = read(io,Int32)
				elseif release == 118
					z = read(io,8)
					v = reinterpret(Int16,z[1:2])[1]
					o = (reinterpret(Int64,z)[1] >> 16)
				end
				if (v,o) == (0,0)
					df[i,j] = missing
				else
					df[i,j] = strtonull(strls[(v,o)])
				end
			elseif typelist[j] == 65526
				dataitemf64 = read(io,Float64)
				if dataitemf64 > 8.9884656743e307
					df[i,j] = missing
				elseif fmtlist[j] == "%d" || fmtlist[j][1:3] == "%td"
					# convert it to Julia date
					df[i,j] = Date(1960,1,1) + Dates.Day(round(Int,dataitemf64))
				elseif fmtlist[j][1:3] == "%tc" || fmtlist[j][1:3] == "%tC"
					df[i,j] = DateTime(1960,1,1,0,0,0) + Dates.Millisecond(round(Int,dataitemf64))
				else
					df[i,j] = dataitemf64
				end
			elseif typelist[j] == 65527
				dataitemf32 = read(io,Float32)
				if dataitemf32 > 1.70141173319e38
					df[i,j] = missing
				elseif fmtlist[j] == "%d" || fmtlist[j][1:3] == "%td"
					# convert it to Julia date
					df[i,j] = Date(1960,1,1) + Dates.Day(round(Int,dataitemf32))
				elseif fmtlist[j][1:3] == "%tc" || fmtlist[j][1:3] == "%tC"
					df[i,j] = DateTime(1960,1,1,0,0,0) + Dates.Millisecond(round(Int,dataitemf32))
				else
					df[i,j] = dataitemf32
				end
			elseif typelist[j] == 65528
				dataitemi32 = read(io,Int32)
				if dataitemi32 > 2147483620
					df[i,j] = missing
				elseif fmtlist[j] == "%d" || fmtlist[j][1:3] == "%td"
					# convert it to Julia date
					df[i,j] = Date(1960,1,1) + Dates.Day(dataitemi32)
				elseif fmtlist[j][1:3] == "%tc" || fmtlist[j][1:3] == "%tC"
					df[i,j] = DateTime(1960,1,1,0,0,0) + Dates.Millisecond(dataitemi32)
				else
					df[i,j] = dataitemi32
				end
			elseif typelist[j] == 65529
				dataitemi16 = read(io,Int16)
				if dataitemi16 > 32740
					df[i,j] = missing
				elseif fmtlist[j] == "%d" || fmtlist[j][1:3] == "%td"
					# convert it to Julia date
					df[i,j] = Date(1960,1,1) + Dates.Day(dataitemi16)
				else
					df[i,j] = dataitemi16
				end
			elseif typelist[j] == 65530
				dataitemi8 = read(io,Int8)
				if dataitemi8 > 100
					df[i,j] = missing
				else
					df[i,j] = dataitemi8
				end
			end
		end
		# strls will be converted to categorical regardless of `categorize` option
		if typelist[j] == 32768
			categorical!(df,varlist[j])
		end
	end

	return df
end

# """
#     df, label = read_stata(fn::String)

# Converts a Stata datafile into a Julia DataFrame. It produces two memory objects (`df` and `label`).
# The first is the dataframe. The second is a [Lables](https://github.com/mwsohn/Labels.jl) object
# that provide functionality to attach variable and value labels.
# """
# function read_stata(fn::String)

#     dt = ReadStat.read_dta(fn)
#     df = readstat2dataframe(dt)

#     return df,get_labels(dt)
# end

#############################################################################
#
# DataFrame tools
#
#############################################################################
"""
    dfcompress!(df::DataFrame)

Reduce `df`'s memory use by changing the eltype of each column in the `df` to
the type that can accommodate the larget and the smallest values within the same
integer or float class.
"""
function dfcompress!(df::DataFrame)
    for v in propertynames(df)

        # if Array is empty after all missings are dropped
        # drop it from the df
        if length(df[!,v]) == sum(ismissing.(df[!,v]))
            delete!(df,v)
            println(v, " was empty, now deleted")
            continue
        end

        # get the original eltype
        eltype_old = nonmissingtype(eltype(df[!,v]))

        # if string, continue
        if eltype_old == String
            continue
        end

        # compress
        df[!,v] = Stella.acompress(df[:,v])

        if eltype_old != nonmissingtype(eltype(df[!,v]))
            println(v, " was ", eltype_old, ", now ", nonmissingtype(eltype(df[!,v])))
        end
    end
end

"""
    acompress(da::AbstractVector)

Compresses a vector to a smallest numeric type that can hold without loss of information.
"""
function acompress(da::AbstractVector)

    # get the original eltype
    eltype_old = nonmissingtype(eltype(da))

    # string variable - do not compress
    if eltype_old == String
        return da
    end

    # copy the vector
    da2 = copy(da)

    if  eltype_old <: Integer
        # get minimum and maximum values
        varmin = minimum(collect(skipmissing(da2)))
        varmax = maximum(collect(skipmissing(da2)))
        # test if Int8
        if eltype_old != Int8 && varmin >= typemin(Int8) && varmax <= typemax(Int8)
            return convert(Vector{Union{Missing,Int8}},da2)
        elseif eltype_old != Int16 && varmin >= typemin(Int16) && varmax <= typemax(Int16)
            return convert(Vector{Union{Missing,Int16}},da2)
        else
            return da2
        end
    elseif eltype_old <: AbstractFloat
        # first test if the floats are integer numbers
        if sum(isinteger.(collect(skipmissing(da2))) .== false) == 0
            return acompress(convert(Vector{Union{Missing,Int64}},da2))
        end

        # get minimum and maximum values
        if eltype_old == Float64 && minimum(collect(skipmissing(da2))) >= floatmin(Float32) && maximum(collect(skipmissing(da2))) <= floatmax(Float32)
            return convert(Vector{Union{Missing,Float32}},da2)
        else
            return da2
        end
    end
end

function atype(df::DataFrame,v::Symbol)
    # Array type = DA for DataArray, CA for Categorical Array, and UV for Union Vector
    if isdefined(Main,:CategoricalArrays) && typeof(df[!,v]) <: CategoricalArray
        return string("CA (", replace(string(eltype(df[!,v].refs)),"UInt" => ""), ")")
    elseif isdefined(Main,:DataArrays) && typeof(df[!,v]) <: DataArray
         return "DA"
    elseif isdefined(Main,:PooledArrays) && typeof(df[!,v]) <: PooledArray
         return "PA"
    elseif isa(eltype(df[!,v]),Union)
        return "UV" # Union Vector
    else
        return "VC"
    end
end

function etype(df::DataFrame,v::Symbol)
    # Eltype
    if typeof(df[!,v]) <: CategoricalArray
        eltyp = string(eltype(df[!,v].pool.levels))
        if in(eltyp,["String","AbstractString"])
            eltyp = string("Str",getmaxwidth(df[!,v].pool.levels))
        end
    else
        eltyp = string(nonmissingtype(eltype(df[!,v])))
        if in(eltyp,["String","AbstractString"])
            eltyp = string("Str",getmaxwidth(df[!,v]))
        elseif eltyp == "Dates.Date"
            eltyp = "Date"
        end
    end

    return eltyp
end

function eltype2(df::DataFrame,v::Symbol)
    if typeof(df[!,v]) <: CategoricalArray
        return eltype(df[!,v].pool.index)
    end
    return nonmissingtype(eltype(df[!,v]))
end

"""
    desc(df::DataFrame,varnames::Symbol...; labels::Union{Nothing,Label}=nothing)

Displays variables in a dataframe much like `showcols`. It can display additional
attributes such as variable labels, value labels and display formats (not used in Julia)
if an optional `labels` is specified. It mimics Stata's `describe` command.
`labels` is automatically converted from a stata file by `read_stata` function. Or one can
be easily created as described in [Labels](https://github.com/mwsohn/Labels.jl).
"""
function desc(df::DataFrame,varnames::Symbol...; labels::Union{Nothing,Label}=nothing)

    if length(varnames) == 0
        varnames = propertynames(df)
    end

    # number of variables
    numvar = length(varnames)

    # length of variable names and value labels
    varlen = length.(string.(varnames)) # length of variable names
    lablen = zeros(Int,length(varnames)) # length of value labels
    for (i,v) in enumerate(varnames)
        lablen[i] = labels != nothing && lblname(labels,v) != nothing ? length(string(lblname(labels,v))) : 0
    end

    # width for variable names
    maxval = maximum(varlen)
    maxval = maxval < 8 ? 8 : maxval

    # width for variable types
    maxatype = 8
    maxeltype = 7
    maxmiss = 7

    # width for value label names
    maxlab = maximum(lablen)
    maxlab = maxlab < 11 ? 11 : maxlab

    # width for display formats
    # maxformat = maximum(forlen)
    # maxformat = maxformat < 6 ? 6 : maxformat

    # width for the variable index - minimum 3 spaces
    maxobs = length(string(numvar))
    maxobs = maxobs < 3 ? 3 : maxobs
    print(lpad("Num",maxobs),"  ")

    # variable name
    print(rpad("Variable",maxval),"  ")

    # type of array
    print(rpad("Atype",maxatype),"  ")

    # type (eltype)
    print(rpad("Eltype",maxeltype),"  ")

    # percent missing
    print(rpad("Missing",maxmiss),"  ")

    if labels != nothing
        # label
        print("Description")
    end

    print("\n")

    # dashes for a line by itself -- assume 30 characters for "label"
    numdashes = maxobs+maxval+maxatype+maxmiss+maxeltype+4
    if labels == nothing
        println(repeat("-",numdashes+4))
    else
        println(repeat("-",numdashes+30))
    end

    nrows = size(df,1)

    # output dataframe
    dfv = DataFrame(Variable = varnames)
    dfv[!,:ArrayType] = Vector{String}(undef,size(dfv,1))
    dfv[!,:Eltype] = Vector{String}(undef,size(dfv,1))
    dfv[!,:Missing] = Vector{String}(undef,size(dfv,1))
    dfv[!,:Lblname] = Vector{String}(undef,size(dfv,1))
    dfv[!,:Description] = Vector{String}(undef,size(dfv,1))

    for (i,v) in enumerate(collect(varnames))

        # variable name
        varstr = string(v)

	# Array Type
        dfv[i,:ArrayType] = atype(df,v)

        # Eltype
        dfv[i,:Eltype] = etype(df,v)

        # percent missing
	nmiss = nmissing(df[!,v])
        dfv[i,:Missing] = string(round(100 * nmiss/nrows,digits=1),"%")

        print(lpad(string(i),maxobs),"  ",rpad(varstr,maxval),"  ",
            rpad(dfv[i,:ArrayType],maxatype),"  ",
            rpad(dfv[i,:Eltype],maxeltype),"  ",
            lpad(dfv[i,:Missing],maxmiss),"  ")

        if labels != nothing
            dfv[i,:Lblname] = lblname(labels,v) == nothing ? "" : string(lblname(labels,v))
            dfv[i,:Description] = varlab(labels,v)
            print(dfv[i,:Description])
        end

        print("\n")
    end

    return dfv
end

function nmissing(s::AbstractArray)
    return sum(ismissing.(s))
end

function getmaxwidth(s::AbstractArray)
    if isa(s, CategoricalArray) && nonmissingtype(eltype(s)) <: CategoricalString
	return maximum(length.(s.pool.levels))
    end

    if nmissing(s) == size(s,1)
	return 0
    end

    return  maximum(length.(collect(skipmissing(s))))
end

"""
    ds(df::DataFrame, typ::Type, args...)
    ds(df::DataFrame, re::Regex)

Returns an array of Symbols of the columns whose eltype matches the `typ` or the regular expression `re`.
For `String` variables, up to two additional arguments can be specified to indicate the minimum and maximum widths
to be returned. In Julia, string variables does not have an associated widths as in Stata.
So, in this function, the width is measured by the length of the longest string. The function follows `Number` type
structure. For example, to list all numeric variables, use `Number`; for all integer variables, use `Integer`; for
all floating point numbers, use `AbstractFloat`. This function is useful to subset a DataFrame using
eltypes or column names. This function emulates the Stata `ds` command.

# Example
```
julia> showcols(aapl)
6081×7 DataFrames.DataFrame
│ Col # │ Name      │ Eltype  │ Missing │
├───────┼───────────┼─────────┼─────────┤
│ 1     │ Date      │ String  │ 0       │
│ 2     │ Open      │ Float64 │ 0       │
│ 3     │ High      │ Float64 │ 0       │
│ 4     │ Low       │ Float64 │ 0       │
│ 5     │ Close     │ Float64 │ 0       │
│ 6     │ Volume    │ Int64   │ 0       │
│ 7     │ Adj_Close │ Float64 │ 0       │

julia> ds(aapl,Int64)
1-element Array{Symbol,1}:
 :Volume

julia> ds(aapl,Float64)
5-element Array{Symbol,1}:
 :Open
 :High
 :Low
 :Close
 :Adj_Close

julia> ds(aapl,Number)
6-element Array{Symbol,1}:
 :Open
 :High
 :Low
 :Close
 :Volume
 :Adj_Close

julia> ds(aapl,String)
1-element Array{Symbol,1}:
 :Date

julia> ds(aapl,r".*lose")
2-element Array{Symbol,1}:
 :Close
 :Adj_Close

```
"""
function ds(df::DataFrame, typ::Type, args...)

    dslist = Vector{Symbol}()

    for v in propertynames(df)
        if typ == String && nonmissingtype(eltype(df[!,v])) == String
            if length(args) == 0
                push!(dslist,v)
                continue
            end

            maxlen = getmaxwidth(df[!,v])
            if length(args) == 1
                if args[1] <= maxlen
                    push!(dslist,v)
                end
            elseif length(args) == 2
                if args[1] <= maxlen <= args[2]
                    push!(dslist,v)
                end
            end
        elseif nonmissingtype(eltype(df[!,v])) == typ
            push!(dslist,v)
        elseif (typ == Integer && nonmissingtype(eltype(df[!,v])) <: Integer) || (typ == AbstractFloat && nonmissingtype(eltype(df[!,v])) <: AbstractFloat) || (typ == Number && nonmissingtype(eltype(df[!,v])) <: Number)
            push!(dslist,v)
        end
    end

    return dslist

end
function ds(df::DataFrame,re::Regex)
    dslist = Array{Symbol,1}()

    for v in propertynames(df)
        if occursin(re,string(v))
            push!(dslist,v)
        end
    end

    return dslist
end

"""
    pickone(df::DataFrame,groupvars::Vector{Symbol})

Create a vector that identifies one record in a `groupvars` in the `df` DataFrame.
This function creates a variable similar to the Stata command
`egen byte pickone = tag(groupvar)`. `groupvars` can be an array of Symbols or a single
Symbol.
"""
function pickone(df::DataFrame,groupvars::Vector{Symbol})
    df[:_____obs_____] = collect(1:size(df,1))
    done = zeros(Int8,size(df,1))
    for subdf in groupby(df, groupvars)
        done[subdf[1,:_____obs_____]]=1
    end
    select!(df,Not(:_____obs_____))
    return done
end
pickone(df::DataFrame,groupvar::Symbol) = pickone(df, [groupvar])

"""
    duplicates(df::DataFrame, args::Symbol...; cmd::Symbol = :report)

Reports, tags, or deletes duplicates in a DataFrame for one or more columns. Use `cmd` to
request one of the actions:

- :report - for getting a frequency table for the number of duplicates
- :tag - for identifying rows with duplicate values
- :drop - for dropping all rows with duplicate values except for the first row.
"""
function duplicates(df::DataFrame, args::Union{Symbol,String}... ; cmd::Symbol = :report)

    if in(cmd,[:report,:drop,:tag]) == false
        error("`cmd = `", cmd, "` is not supported.")
    end

    # if args are not specified, use all variables
    nargs = length(args)
    if nargs == 0
        args = tuple(names(df)...)
    end

    gdf = groupby(df,collect(args))
    dfx = combine(gdf,nrow => :__dups__)

    if cmd == :report
        na = freqtable(dfx, :__dups__)
        setdimnames!(na,"copies",1)
        return na
    end

    # substract 1 from :x in dfx (we are reporting 0 for unique observations)
    dfx[:__dups__] .-= 1
    df = leftjoin(df, dfx, on = collect(args))

    if cmd == :tag
        return df[:__dups__]
    end

    # if cmd == :drop
    ba = [x == 1 ? true : false for x in pickone(df,collect(args))]
    return select(df[ba,:],Not(:__dups__))
end

"""
    sample(df::AbstractDataFrame,num::Real)

Creates an Int8 array that identifies a randomly selected sample (1 = sample; 0 otherwise) from the input dataframe. `num`
specifies the sample size. If `num` is an integer, it indicates the number of rows to be selected.
If it is a float, it indicates the percentage of the input dataframe to be randomly selected.
Use `Random.seed!()` to set a seed before using `sample()`.

##Example
To select 100 rows, use

```
julia> df[!,:sample] = sample(df,100)
```

To select 20% of the original dataframe, use

```
julia> df[!,:sample2] = sample(df,.2)
```

"""
function sample(df::AbstractDataFrame,n::Real)
    len = size(df,1)
    asamp = zeros(Int8,len)
    if isa(n,AbstractFloat)
        asamp[randsubseq(collect(1:len),n)] .= 1
    else
        asamp[randperm(len)[1:n]] .= 1
    end
    return asamp
end

"""
    dfsample(df::AbstractDataFrame,num::Real)

Creates a dataframe with a randomly selected sample from the input dataframe. `num`
specifies the sample size. If `num` is an integer, it indicates the number of rows to be selected.
If it is a float, it indicates the percentage of the input dataframe to be randomly selected.
Use `Random.seed!()` to set a seed before running `dfsample()`.

##Example
To select 100 rows, use

```
julia> df2 = dfsample(df,100)
```

To select 20% of the original dataframe, use

```
julia> df2 = dfsample(df,.2)
```

"""
function dfsample(df::AbstractDataFrame,num::Real)
    return df[sample(df,num) .== 1,:]
end

"""
    dfmerge(df1::DataFrame,df2::DataFrame,linkers::Union{Symbol,Vector}; kind::Symbol = :outer)

Produces a merged dataframe with a `:_merge` variable indicating how the merge was done.
This function is a wrapper for `join` in the DataFrames package and mimics `merge ... using ...`
command in Stata. Currently, it supports a merge of two sources only.

The default merge is a `outer` join that keeps records from both sources in the merged
dataframe. `:_merge` values indicates:

- `1`: in the first source alone
- `2`: in the second source alone
- `3`: in both sources. A successful merge was achieved.

It is the responsibility of the user to keep or discard records with `:_merge` values 1 or 2.

The third argument can be either a single symbol or a vector of symbols to be used as the
linkage variables.

Optionally, `kind` argument can be specified. Available options are
`:left`, `:right`, `:inner`,`:outer`,`:semi`, and `:anti`. Consult the documentation for
the `join` function in the DataFrames package.
"""
function dfmerge(df1::DataFrame,df2::DataFrame,linkers::Union{Symbol,Vector};kind::Symbol = :outer)
    if in(:_merge,propertynames(df1))
        error("`:_merge' exists in the first dataframe")
    end
    if in(:_merge,propertynames(df2))
        error("`:_merge' exists in the second dataframe")
    end

    df1[:___mergeleft___] = ones(Int8,size(df1,1))
    df2[:___mergeright___] = ones(Int8,size(df2,1))

    df_merged = join(df1,df2,on = linkers,kind=kind)
    df_merged[:_merge] = zeros(Int8,size(df_merged,1))
    for i = 1:size(df_merged,1)
        if ismissing(df_merged[i,:___mergeright___])
            df_merged[i,:_merge] = 1
        elseif ismissing(df_merged[i,:___mergeleft___])
            df_merged[i,:_merge] = 2
        elseif df_merged[i,:___mergeleft___] == 1 && df_merged[i,:___mergeright___] == 1
            df_merged[i,:_merge] = 3
        end
    end
    delete!(df_merged,[:___mergeleft___,:___mergeright___])
    return df_merged
end

function get_class(val::Real,thresholds::Vector,lower::Bool)
    if lower == false
        for i = 1:length(thresholds)
            if i == 1 && val < thresholds[1]
                return 1
            elseif i > 1 && thresholds[i-1] <= val < thresholds[i]
                return i
            end
        end
    else
        for i = 1:length(thresholds)
            if i == 1 && val <= thresholds[1]
                return 1
            elseif i > 1 && thresholds[i-1] < val <= thresholds[i]
                return i
            end
        end
    end
    return length(thresholds)+1
end

"""
    classify(da::AbstractVector,thresholds::Vector; lower::Bool = false)

Produces a categorical variable based on a data array and a vector of thresholds.
Use `lower = true` to classify the threshold value to the lower class.
"""
function classify(da::AbstractVector,thresholds::Vector; lower::Bool = false)
    da2 = Vector{Union{Missing,Int8}}(undef,length(da))
    if length(thresholds) > 100
        error("Cannot use more than 100 threshold values.")
    end

    for i = 1:size(da,1)
        if ismissing(da[i])
            da2[i] = missing
        else
            da2[i] = get_class(da[i],thresholds,lower)
        end
    end
    return da2
end

"""
    addterms(fm::Formula,v::Vector{Symbol})

Adds covariates to an existing formula. `fm` is an object of `Formula` type
created by `@formula` macro. `v` is a vector of Symbols. This function is useful to
build hierarchically nested models.

### Example

```
julia> using DataFrames

julia> fm = @formula(income ~ age + male)
Formula: income ~ age + male

julia> fm2 = addterms(fm,[:race, :educ])
Formula: income ~ age + male + race + educ
```

"""
function addterms(fmm::Formula,v::Vector{Symbol})

    # create a new Formula object
    fm = deepcopy(fmm)

    # if fm.rhs is a null modle with rhs == 1 or
    # a unadjusted model with only one rhs variable,
    # convert it to an Expr first
    if fm.rhs == 1
        fm.rhs = Expr(:call,:+,Tuple(v)...)
    elseif typeof(fm.rhs) == Symbol
        fm.rhs = Expr(:call,:+,fm.rhs,Tuple(v)...)
    else
        push!(fm.rhs.args,Tuple(v)...)
    end

    return fm
end

"""
    Stella.formula(o::Symbol,v::Vector{Symbol})

Creates a Formula from a vector of covariates. This function is not exported for potential conflict with other packages.

- o - A symbol for the dependent variable
- v - A vector of symbols for covariates

## Usage

julia> fm = Stella.formula(:readmit, [:age,:male,:race,:chd,:exercise])
Formula: readmit ~ age + male + race + chd + exercise
"""
function formula(o::Symbol,v::Vector{Symbol})
    if length(v) == 1
        rhs = v[1]
    else
        rhs = Expr(:call,:+,Tuple(v)...)
    end
    e = :($o ~ $rhs)
    return Formula(e,e,o,rhs)
end

"""
    renvars!(df::DataFrame; vars::Array{Symbol,1}, case = "lower")

Rename column names in `vars` to either upper or lower cases. The default is to convert
all columns to lower case names.
"""
function renvars!(df::DataFrame; vars=[], case="lower")
    numvar = length(vars)
    symnames = propertynames(df)

    if numvar == 0
        varnames = propertynames(df)
    else
        varnames = propertynames(df[vars])
    end

    for nm in varnames

        if case == "lower" || case == "LOWER"
            newname = lowercase(string(nm))
        elseif case == "upper" || case == "UPPER"
            newname = uppercase(string(nm))
        else
            error(option," ", case," is not allowed.")
        end
        if nm != Symbol(newname)
            rename!(df,nm => Symbol(newname))
        end
    end
end

vidx(df::DataFrame,varname::Symbol) = findfirst(x->x == varname, propertynames(df))

"""
    destring(da::AbstractArray;force=true)
    destring(df::DataFrmae,strvar::Symbol;force=true)
    destring!(df::DataFrame,strvars; newvars::Vector{Symbol} = [], force=true, replace=false)

Convert a string vector to a numeric vector. Use `force = true` to coerce conversion of alphanumeric strings to
missing values. If `force = false`, any vector containing non-numeric characters are not converted.
Use `replace = true` in `destring!` to replace the original string vector with a new converted numeric vector.
If `replace` option is specified, `newvars` array is ignored.
"""
function destring(da::AbstractArray; force=true)
    if length(da) == 0
        error("Input array is empty!")
    end
    if nonmissingtype(eltype(da)) <: Number
        error("Input array is a numeric vector. No conversion needed.")
    end

    # check if the values include any alphabetic characters or decimals
    isfloat = false
    alpha = false
    da_safe = collect(skipmissing(da))
    for i in length(da_safe)
        if sum([isdigit(x) == false for x in da_safe[i]]) > 0
            alpha = true
        end
        if occursin(r"[,0-9]*\.[0-9]+",da_safe[i])
            isfloat = true
        end
    end

    if alpha && force == false
        error("Input array contains alphabetic letters. Use 'force=true' option to coerce conversion.")
    end

    T = isfloat ? Float64 : Int64
    da2 = [ismissing(x) ? missing : parse(T,x) for x in da]

    return acompress(da2)
end
destring(df::AbstractDataFrame,strvar::Symbol; force=true) = destring(df[:,strvar],force=force)
function destring!(df::DataFrame,strvars; newvars::Vector{Symbol} = [], force=true, replace=false)

    if replace
        for v in strvars
            df[v] = destring(df[v],force=force)
        end
    else
        # check if there are same number of symbols in strvars
        if length(strvars) != length(newvars)
            error("The number of symbols in ", strvars, " and ", newvars, " are not the same.")
        end
        for i in 1:length(strvars)
            df[newvars[i]] = destring(df[strvars[i]], force=force)
        end
    end
end

# """
#     rowsum(df::DataFrame)
#
# Creates a DataArray that contains the row total of all values on the same row of `df`.
# If one of the columns contain an NA value on a row, an NA value will be returned for that
# row. This function emulates Stata's `egen rowsum = rowtotal(var1 - var3)`.
#
# ```
# julia>df[:rowsum] = df[[:var1,:var2,:var3]]
# ```
#
# If the position numbers for `:var1` (e.g., 4), `:var2` (5), `:var3` (6) are known and consecutive,
# you can specify them as follows:
#
# ```
# julia>df[:rowsum] = df[collect(4:6)]
# ```
# """
# function rowsum(df::DataFrame)
#
#     isfloat = false
#     for i in 1:size(df,2)
#         if eltype(df[i]) <: AbstractFloat
#             isfloat = true
#         end
#     end
#
#     if isfloat
#         da = zeros(Union{Missing,Float64},size(df,1))
#     else
#         da = zeros(Union{Missing,Int64},size(df,1))
#     end
#
#     ba = completecases(df)
#     for i = 1:size(df,1)
#         if ba[i] == false
#             da[i] = missing
#         else
#             for j = 1:size(df,2)
#                 da[i] += df[i,j]
#             end
#         end
#     end
#     return da
# end

"""
    rowstat(df::DataFrame,func::Function)

Creates a vector that contains the row statistic of all values on the same row of `df`
produced by the `func` function. If one of the columns contain a missing value on a row, the returned value will be mssing for that
row. This function emulates Stata's `egen` row functions such as `rowtotal`, `rowmean`, etc.

```
julia>df[:rowmean] = rowstat(df[[:var1,:var2,:var3]],mean)
```
"""
function rowstat(df::DataFrame,func::Function)

    # df size
    nrow, ncol = size(df)

    # array to be returned
    da = zeros(Union{Missing,Float64},size(df,1))

    # temporary array of row values
    ta = Vector{Union{Missing,Float64}}(undef,size(df,2))

    # iterate over all rows
    for i = 1:nrow

        # counter for non-missing values in each row
        k=0
        for j = 1:ncol
            if ismissing(df[i,j]) == false
                k += 1
                ta[k] = df[i,j]
            end
        end

        # if any of the columns has a missing value, return missing
        if k < ncol
            da[i] = missing
        else
            tmpfloat = func(ta[1:k])
            da[i] = isnan(tmpfloat) ? missing : tmpfloat
        end
    end
    return da
end

"""
    xtile(da::AbstractArray;nq::Int = 4, lower::Bool = false, cutoffs::Union{Nothing,AbstractVector} = nothing)
    xtile(df::DataFrame,varname::Symbol;nq::Int = 4, lower::Bool = false, cutoffs::Union{Nothing,AbstractVector} = nothing)

Create a DataArray{Int8,1} that identifies `nq` categories based on values in `da`.
The default `nq` is 4. `cutoffs` vector can be provided to make custom categories.
`cutoffs` vector is expected to contain `nq - 1` elements. The minimum and maximum values
will be computed. `lower = true` will make the threshold values to go in the lower category (default: false).

```
julia> df[:agecat] = xtile(df[:age], nq = 3)
```
"""
function xtile(da::AbstractArray ; nq::Int = 4, lower::Bool = false, cutoffs::Union{Nothing,AbstractVector} = nothing)

    if nq == 1
        error("`nq` must be greater than 1")
    end

    da2 = collect(skipmissing(da))
    if cutoffs == nothing
	    cutoffs = nquantile(da2,nq)
    elseif length(cutoffs) == nq - 1
        cutoffs = vcat(minimum(da2), cutoffs, maximum(da2))
    else
        error("`cutoffs` vector length is not consistent with `nq`. It must be 1 greater or 1 less than `nq`.")
    end

    if lower
        return [ismissing(x) ? missing : qval_low(x,cutoffs) for x in da2]
    else
        return [ismissing(x) ? missing : qval_high(x,cutoffs) for x in da2]
    end
end
xtile(df::DataFrame,arg::Symbol; nq::Int = 4, cutoffs::Union{Nothing,AbstractVector} = nothing) = xtile(df[!,arg], nq = nq, cutoffs = cutoffs)

function qval_high(val::Real,cut::Vector)
    cl = length(cut)

    for i in 2:cl-1
        if i == 2 && val < cut[i]
            return 1
        elseif i == cl-1 && cut[i] <= val
            return i
        elseif cut[i] <= val < cut[i+1]
            return i
        end
    end
    warn("Error - check qval function")
end

function qval_low(val::Real,cut::Vector)
    cl = length(cut)

    for i in 2:cl-1
        if i == 2 && val <= cut[i]
            return 1
        elseif i == cl-1 && cut[i] <= val
            return i
        elseif cut[i] < val <= cut[i+1]
            return i
        end
    end
    warn("Error - check qval function")
end

"""
    rapply(df::DataFrame,r::OrderedDict)

Creates a recoded vector of values according to the rule set specified as an ordered dictionary.
Keys in the rule set must be the recoded value and the values must be the rules.
The rule must be an expression or a string that return a boolean value (`true` or `false`).

## Example

julia> ruleset = OrderedDict(
            1 => :( df[:race] .== 1 && df[:hispanic] .== 0),
            2 => :( df[:race] .== 2 && df[:hispanic] .== 0),
            3 => :( df[:hispanic] .== 1),
            4 => :( df[:hispanic] .== 0 && in.(df[:race],[1,2]) == false)
       )
OrderedDict{Int64,Expr} with 4 entries:
  1 => :(df[:race] .== 1 && df[:hispanic] .== 0)
  2 => :(df[:race] .== 2 && df[:hispanic] .== 0)
  3 => :(df[:hispanic] .== 1)
  4 => :(df[:hispanic] .== 0 && in.(df[:race], [1, 2]) == false)

julia> recode(adf,ruleset)

"""
function rapply(df::DataFrame,r::OrderedDict)

    # values
    vals = sort(collect(keys(r)))

    # types
    vtyp = eltype(vals)

    # empty Vector
    vec = Vector{Union{vtyp,Missing}}(undef,size(df,1))

    # go through the rules and assign values
    for v in vals

        if typeof(r[v]) == String
            ba = eval(parse(r[v]))
        elseif typeof(r[v]) == Expr
            ba = eval(r[v])
        else
            error(typeof(r[v]), " is not an allowed type for rule ", string(r[v]))
        end

        for i in findall(x->x==true,ba)
            vec[i] = v
        end
    end

    return vec
end

"""
	uncategorical!(df::AbstractDataFrame,vv::Union{Symbol,Vector{Symbol}})

uncategorizes a CategoricalArray by replacing the refs with its original values.
"""
function uncategorical!(df::AbstractDataFrame,vv::Union{Symbol,Vector{Symbol}})
    for v in vcat(vv)
        if isa(df[:,v],CategoricalArray)
            lev = df[:,v].pool.levels
            df[!,v] = [x == 0 ? missing : lev[x] for x in df[:,v].refs]
        end
    end
end

"""
	uncategorize(v::AbstractArray)

produces an array by uncategorizing a CategoricalArray by replacing the refs with its original values.
"""
function uncategorize(v::CategoricalArray)
    return [x == 0 ? missing : v.pool.levels[x] for x in v.refs]
end

# convert the Arrow data types from the feather file
function convert_feather(df::DataFrame)
    for v in propertynames(df)
        if typeof(df[!,v]) <: DictEncoding
            df[!,v] = CategoricalArray(df[!,v])
	elseif eltype(df[!,v]) == Feather.Arrow.Datestamp
	    df[!,v] = convert.(Date,df[!,v])
        else
            df[!,v] = Array(df[!,v])
        end
    end

    return df
end


text:
#include("../src/ttv.jl")
#include("/Users/ericagol/Computer/Julia/regress.jl")

@testset "ttvbv_elements" begin

# This routine takes derivative of transit times, vsky & bsky^2 with respect
# to the initial orbital elements.
#n = 8
n = 3
H = [3,1,1]
n_body = n
t0 = 7257.93115525-7300.0
#h  = 0.12
h  = 0.04
#tmax = 600.0
#tmax = 800.0
#tmax = 1000.0
tmax = 10.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')
elements[:,3] .-= 7300.0
# Make masses of planets bigger
#elements[2,1] *= 10.0
#elements[3,1] *= 10.0

ntt = zeros(Int64,n)

# Make an array, ttbv,  to hold transit times, vsky & b^2:
# First, though, make sure it is large enough:
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
dtbvdq0 = zeros(3,n,maximum(ntt),7,n)
ttbv  = zeros(3,n,maximum(ntt))
ttbv1 = zeros(3,n,maximum(ntt))
ttbv2 = zeros(3,n,maximum(ntt))
ttbv3 = zeros(3,n,maximum(ntt))
ttbv4 = zeros(3,n,maximum(ntt))
ttbv8 = zeros(3,n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12
dq = ttvbv_elements!(H,t0,h,tmax,elements,ttbv1,count1,0.0,0,0,rstar)
# Now call with half the timestep:
count2 = zeros(Int64,n)
count3 = zeros(Int64,n)
dq = ttvbv_elements!(H,t0,h/10.,tmax,elements,ttbv2,count2,0.0,0,0,rstar)

mask = zeros(Bool, size(dtbvdq0))
for itbv=1:3, jq=1:n_body
  for iq=1:7
    if iq == 7; ivary = 1; else; ivary = iq+1; end  # Shift mass variation to end
    for i=2:n
      for k=1:count2[i]
        # Ignore inclination & longitude of nodes variations:
        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
          mask[itbv,i,k,iq,jq] = true
        end
      end
    end
  end
end

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtbvdelements0 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements0 = ttvbv_elements!(H,t0,h,tmax,elements,ttbv,count,dtbvdq0,rstar)
dtbvdq2 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements2 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements2 = ttvbv_elements!(H,t0,h/2.,tmax,elements,ttbv2,count,dtbvdq2,rstar)
dtbvdq4 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements4 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements4 = ttvbv_elements!(H,t0,h/4.,tmax,elements,ttbv4,count,dtbvdq4,rstar)
dtbvdq8 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements8 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements8 = ttvbv_elements!(H,t0,h/8.,tmax,elements,ttbv8,count,dtbvdq8,rstar)
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdelements0[mask])-asinh.(dtbvdelements2[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdq0)-asinh.(dtbvdq2))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdelements2[mask])-asinh.(dtbvdelements4[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdq2)-asinh.(dtbvdq4))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdelements4[mask])-asinh.(dtbvdelements8[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdq4)-asinh.(dtbvdq8))))

# Compute derivatives numerically:
#nq = 15
# This "summarizes" best numerical derivative:
dtbvdelements0_num = zeros(BigFloat,3,n,maximum(ntt),7,n)

# Compute derivatives with BigFloat for additional precision:
elements0 = copy(elements)
dq0 = big(1e-10)
t0big = big(t0); tmaxbig = big(tmax); hbig = big(h)
zilch = big(0.0)
# Compute the transit times in BigFloat precision:
ttbv_big = big.(ttbv); elementsbig = big.(elements0); rstarbig = big(rstar)
dqbig = ttvbv_elements!(H,t0big,hbig,tmaxbig,elementsbig,ttbv_big,count,big(0.0),0,0,rstarbig)

mask = zeros(Bool, size(dtbvdq0))
for itbv=1:3, jq=1:n_body
  for iq=1:7
    if iq == 7; ivary = 1; else; ivary = iq+1; end  # Shift mass variation to end
    for i=2:n
      for k=1:count2[i]
        # Ignore inclination & longitude of nodes variations:
        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
          mask[itbv,i,k,iq,jq] = true
        end
      end
    end
  end
end

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtbvdelements0 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements0 = ttvbv_elements!(H,t0,h,tmax,elements,ttbv,count,dtbvdq0,rstar)
dtbvdq2 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements2 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements2 = ttvbv_elements!(H,t0,h/2.,tmax,elements,ttbv2,count,dtbvdq2,rstar)
dtbvdq4 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements4 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements4 = ttvbv_elements!(H,t0,h/4.,tmax,elements,ttbv4,count,dtbvdq4,rstar)
dtbvdq8 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements8 = zeros(3,n,maximum(ntt),7,n)
dtbvdelements8 = ttvbv_elements!(H,t0,h/8.,tmax,elements,ttbv8,count,dtbvdq8,rstar)
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdelements0[mask])-asinh.(dtbvdelements2[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdq0)-asinh.(dtbvdq2))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdelements2[mask])-asinh.(dtbvdelements4[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdq2)-asinh.(dtbvdq4))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdelements4[mask])-asinh.(dtbvdelements8[mask]))))
println("Maximum error on derivative: ",maximum(abs.(asinh.(dtbvdq4)-asinh.(dtbvdq8))))

# Compute derivatives numerically:
#nq = 15
# This "summarizes" best numerical derivative:
dtbvdelements0_num = zeros(BigFloat,3,n,maximum(ntt),7,n)

# Compute derivatives with BigFloat for additional precision:
elements0 = copy(elements)
dq0 = big(1e-10)
ttbv2 = big.(ttbv2)
ttbv3 = big.(ttbv3)
t0big = big(t0); tmaxbig = big(tmax); hbig = big(h)
zilch = big(0.0)
# Compute the transit times in BigFloat precision:
ttbv_big = big.(ttbv); elementsbig = big.(elements0)
# Now, compute derivatives numerically:
for jq=1:n_body
  for iq=1:7
    elementsbig = big.(elements0)
    if iq == 7; ivary = 1; else; ivary = iq+1; end  # Shift mass variation to end
    elementsbig[jq,ivary] += dq0
    dq_plus = ttvbv_elements!(H,t0big,hbig,tmaxbig,elementsbig,ttbv2,count2,zilch,0,0,big(rstar))
    elementsbig[jq,ivary] -= 2dq0
    dq_minus = ttvbv_elements!(H,t0big,hbig,tmaxbig,elementsbig,ttbv3,count2,zilch,0,0,big(rstar))
    #xm,vm = init_nbody(elements,t0,n_body)
    for i=2:n
      for k=1:count2[i]
        # Compute double-sided derivative for more accuracy:
        for itbv=1:3
          dtbvdelements0_num[itbv,i,k,iq,jq] = (ttbv2[itbv,i,k]-ttbv3[itbv,i,k])/(2dq0)
        # Ignore inclination & longitude of nodes variations:
          if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
            mask[itbv,i,k,iq,jq] = true
          end
        end
      end
    end
  end
end

println("Max diff asinh(dtbvdelements): ",maximum(abs.(asinh.(dtbvdelements0[mask])-asinh.(dtbvdelements0_num[mask]))))

#=
using PyPlot

# Make a plot of the fractional errors:
for itbv=1:3, i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtbvdelements0_num[itbv,i,2:count1[i],k,l])) > 0
    loglog(ttbv[itbv,i,2:count1[i]].-ttbv[itbv,i,1],convert(Array{Float64,1},abs.(dtbvdelements0[itbv,i,2:count1[i],k,l]./dtbvdelements0_num[itbv,i,2:count1[i],k,l].-1.)))
  end
end

clf()
# Plot the difference in the TTVs:
for itbv=1:3, i=2:3
  diff1 = convert(Array{Float64,1},abs.(ttbv1[itbv,i,2:count1[i]]./ttbv_big[itbv,i,2:count1[i]].-1.0));
  loglog(ttbv[itbv,i,2:count1[i]].-ttbv[itbv,i,1],diff1);
end
for itbv=1:3, i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtbvdelements0_num[itbv,i,2:count1[i],k,l])) > 0 && k != 5 && k != 6
    diff1 = convert(Array{Float64,1},abs.(dtbvdelements0[itbv,i,2:count1[i],k,l]./dtbvdelements0_num[itbv,i,2:count1[i],k,l].-1));
    diff2 = convert(Array{Float64,1},abs.(asinh.(dtbvdelements0[itbv,i,2:count1[i],k,l])-asinh.(dtbvdelements0_num[itbv,i,2:count1[i],k,l])));
    loglog(ttbv[itbv,i,2:count1[i]].-ttbv[itbv,i,1],diff1);
    println(i," ",k," ",l," frac error: ",convert(Float64,maximum(diff1))," asinh error: ",convert(Float64,maximum(diff2))); #read(STDIN,Char);
  end
end

# Plot a line that scales as time^{3/2}:

loglog([1.0,1024.0],1e-09*[1,2^15],":")
loglog([1.0,1024.0],1e-12*[1,2^15],":")
loglog([1.0,1024.0],1e-15*[1,2^15],":")
=#

@test isapprox(asinh.(dtbvdelements0[mask]),asinh.(dtbvdelements0_num[mask]);norm=maxabs)
dtbvdelements0_num = convert(Array{Float64,5},dtbvdelements0_num)
#@save "ttvbv_elements.jld2" dtbvdelements0 dtbvdelements0_num ttbv
end


text:
using LowRankApprox

@testset "Truncation with standard SVD works correctly" begin
    D = 100
    Dcut = D - 1
    tol = 1E-8

    a = rand(D, D)

    U1, Σ1, V1 = svd(a)

    δ = min(Dcut, size(Σ1)...)
    U1 = U1[:, 1:δ]
    Σ1 = Σ1[1:δ]
    V1 = V1[:, 1:δ]

    U2, Σ2, V2 = svd(a)

    δ = min(Dcut, size(Σ2)...)
    U2 = U2[:, 1:δ]
    Σ2 = Σ2[1:δ]
    V2 = V2[:, 1:δ]

    r1 = U1 * Diagonal(Σ1) * V1'
    r2 = U2 * Diagonal(Σ2) * V2'

    @test norm(r1-r2) < tol
end

@testset "Truncation with with random SVD works correctly" begin

    D = 100
    Dcut = D - 1
    tol = 1E-8

    a = rand(D, D)

    U1, Σ1, V1 = psvd(a, rank=Dcut, atol=1E-16, rtol=1E-16)
    U2, Σ2, V2 = psvd(a, rank=Dcut, atol=1E-16, rtol=1E-16)

    r1 = U1 * Diagonal(Σ1) * V1'
    r2 = U2 * Diagonal(Σ2) * V2'

    println(norm(a-r2))
    println(norm(a-r1))

    @test norm(r1-r2) < tol
end


text:
function prob(n)
  n = BigFloat(n)
  two = BigFloat(2)
  x = n/2
  a = factorial(n)
  b = factorial(x)
  b *= b
  c = two^n
  return a / b / c
end

for n in [3030, 1000, 500, 200, 100, 50, 10]
  p = float64(round(prob(n)*100, 2))
  println(n, " : ", p, " %")
end

max_stars_repo_name: jgoizueta/binomial-techniques

max_issues_repo_name: jgoizueta/binomial-techniques

max_forks_repo_name: jgoizueta/binomial-techniques


text:
"""
    degrees_submatrices(exponents_matrix)

Returns the submatrices with the corresponding degree. This takes the degrees
from the bottom and assume a reverse lexicographic order of the columns.

## Example
```julia
julia> degrees_submatrices([1 2 2; 1 3 2; 1 1 2])
([1, 2], [[1 2; 1 3], [1 ; 1]])
"""
function degrees_submatrices(E)
    submatrices = Vector{typeof(@view E[1:end-1, 1:1])}()# Vector{Matrix{Int}}()
    degrees = Int[]
    j = 1
    last_d_change = 1
    d = E[end, 1]
    n = size(E, 2)
    while j < n
        if E[end,j + 1] != d
            push!(submatrices, @view E[1:end-1, last_d_change:j])
            push!(degrees, d)
            d = E[end,j + 1]
            last_d_change = j+1
        end
        j += 1
    end
    push!(submatrices, @view E[1:end-1, last_d_change:end])
    push!(degrees, E[end, end])
    degrees, submatrices
end

pow(x::AbstractFloat, k::Integer) = Base.FastMath.pow_fast(x, k)
# simplified from Base.power_by_squaring
function pow(x::Number, p::Integer)
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end

"""
    diffvars(variables, parameters)

Remove `parameters` from `variables`. Deals with the case that `parameters` is `nothing`.
"""
diffvars(variables, parameters) = setdiff(variables, parameters)
diffvars(variables, ::Nothing) = variables

# Implementation from Base.sort adapted to also reorder an associated vector
"""
    revlexicographic_cols_perm(A, v)

Sorts the columns of `A` in reverse lexicographic order and returns the permutation vector
to obtain this ordering.
"""
function revlexicographic_cols_perm(A::AbstractMatrix; kws...)
    inds = axes(A,2)
    cols = map(i -> (@view A[end:-1:1, i]), inds)
    sortperm(cols; kws...)
end

@generated function assemble_matrix(vs::SVector{M, SVector{N, T}}) where {T, N, M}
    quote
        SMatrix{$M, $N, $T, $(M*N)}(
            tuple($([:(vs[$i][$j]) for j=1:N for i=1:M]...))
        )
    end
end


text:
function Atanh()
    return Atanh(())
end

function value(obj::Atanh, arg0::Derivative)
    return jcall(obj, "value", Derivative, (Derivative,), arg0)
end

function value(obj::Atanh, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Atanh, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
module JULES

using Oceananigans
using Oceananigans: AbstractGrid

using Oceananigans.TurbulenceClosures:
    IsotropicDiffusivity, DiffusivityFields, with_tracers

export
    Entropy, Energy,
    DryEarth, DryEarth3,
    CompressibleModel,
    update_total_density!,
    time_step!,
    cfl, acoustic_cfl

include("Operators/Operators.jl")

include("lazy_fields.jl")
include("thermodynamics.jl")
include("pressure_gradients.jl")
include("time_steppers.jl")
include("compressible_model.jl")

include("source_terms.jl")
include("time_stepping_kernels.jl")
include("time_stepping.jl")
include("utils.jl")

end # module


text:
@symbol_func function P_E(cur_reactor::AbstractReactor)
  cur_P = P_T(cur_reactor)

  cur_P *= cur_reactor.eta_T

  cur_P
end


text:
# balance a dataset by oversampling:
# works for vectors, arrays (by last dim) and dataframes (by row).
# returns new x, y-Tuple:
#
function do_balance(x,y)

    return MLDataUtils.oversample((x, y))
end

function do_split(x,y; at=0.5)

    # return MLDataUtils.splitobs(MLDataUtils.shuffleobs((x, y)), at=at)
    return MLDataUtils.stratifiedobs((x, y), p=at)
end

function do_shuffle(x,y)

    return collect(MLDataUtils.shuffleobs((x, y)))
end

"""
    function crop_array(x, crop_sizes)

Crop a n-dimensional array to the given size. Cropping is always
centered (i.e. a margin is removed).

### Arguments:
+ `x`: n-dim AbstractArray
+ `crop_sizes`: Tuple of target sizes to which the array is cropped.
        Allowed values are Int or `:`. If `crop_sizes` defines less
        dims as x has, the remaining dims will not be cropped (assuming `:`).
        If a demanded crop size is bigger as the actual size of x,
        it is ignored.
"""
function crop_array(x, crop_size)

    x_size = size(x)
    ranges = []

    for (cs,s) in zip(crop_size, x_size)
        if cs isa Colon
            push!(ranges, :)
        elseif cs >= s
            push!(ranges, :)
        else
            margin = fld(s-cs, 2)
            push!(ranges, margin+1:margin + cs)
        end
    end
    while length(ranges) < length(x_size)
        push!(ranges, :)
    end
    return x[ranges...]
end

"""
    function init0(siz...)

Initialise a vector or array of size `siz` with zeros.
If a GPU is detected type of the returned value is `KnetArray{Float32}`,
otherwise `Array{Float32}`.

### Examples:
```
julia> init0(2,10)
2×10 Array{Float32,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

 julia> init0(0,10)
 0×10 Array{Float32,2}
```
"""
function init0(siz...)

    x = zeros(Float32, siz...)

    if CUDA.functional()
        return KnetArray(x)
    else
        return x
    end
end

"""
    function convert2KnetArray(x)

Convert an array `x` to a `KnetArray{Float32}` only in GPU context
(if `CUDA.functional()`) or to an `Array{Float32}` otherwise.
"""
function convert2KnetArray(x)

    if CUDA.functional()
        return Knet.KnetArray{Float32}(x)
    else
        return Array{Float32}(x)
    end
end

"""
function blowup_array(x, n)

Blow up an array `x` with an additional dimension
and repeat the content of the array `n` times.

### Arguments:
+ `x`: Array of any dimension
+ `n`: number of repeats. ´n == 1´ will return an
array with an addritional dimension of size 1.

### Examples:

```Julia
julia> x = [1,2,3,4]; blowup_array(x, 3)
4×3 Array{Int64,2}:
 1  1  1
 2  2  2
 3  3  3
 4  4  4

julia> x = [1 2; 3 4]; blowup_array(x, 3)
2×2×3 Array{Int64,3}:
[:, :, 1] =
 1  2
 3  4

[:, :, 2] =
 1  2
 3  4

[:, :, 3] =
 1  2
 3  4
```
"""
function blowup_array(x, n)

    siz = size(x)
    out = x
    i = 1
    while i < n
        out = cat1d(out,x)
        i += 1
    end
    return reshape(out, siz...,:)
end

"""
function recycle_array(x, n)

Recycle an array `x` along the last dimension
and repeat the content of the array `n` times.
The number of dims stays unchanged, but the array
valueas are repeated `n` times.

### Arguments:
+ `x`: Array of any dimension
+ `n`: number of repeats. ´n == 1´ will return an unchanged
        array.

### Examples:

```Julia
julia> recycle_array([1,2],3)
6-element Array{Int64,1}:
 1
 2
 1
 2
 1
 2

julia> x = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> recycle_array(x,3)
2×6 Array{Int64,2}:
 1  2  1  2  1  2
 3  4  3  4  3  4
```
"""
function recycle_array(x, n)

    dim = ndims(x)
    out = x
    i = 1
    while i < n
        out = cat(out,x, dims=dim)
        i += 1
    end
    return out
end

# split the iterator in 2 parts at at
# and return the first and the second part as separate iterators.
#
# function split_iterator(itr, at)
#
#     len = length(itr)
#     last_trn = Int(round(len * at))
#     return Iterators.take(itr, last_trn), Iterators.drop(itr, last_trn)
# end

"""
    function de_embed(x; remove_dim=true)

Replace the maximum of the first dimension of an n-dimensional array
by its index (aka argmax()).
If `remove_dim=false` the first dim is preserved with size=1; otherwise
the returned array has the first dimension removed.
"""
function de_embed(x; remove_dim=true)

    siz = size(x)
    depth = siz[1]
    siz = siz[2:end]

    x = reshape(x, depth, :)
    x = softmax(x, dims=1)
    x = [argmax(x[:,i]) for i in 1:size(x)[2]]
    if remove_dim
        return reshape(x, siz...)
    else
        return reshape(x, 1,siz...)
    end
end


text:
mutable struct Scope
    parent::Union{Scope,Nothing}
    expr::EXPR
    names::Dict{String,Binding}
    modules::Union{Nothing,Dict{Symbol,Any}}
    ismodule::Bool
end
Scope(expr) = Scope(nothing, expr, Dict{Symbol,Binding}(), nothing, typof(expr) === CSTParser.ModuleH || typof(expr) === CSTParser.BareModule)
function Base.show(io::IO, s::Scope)
    printstyled(io, typof(s.expr))
    printstyled(io, " ", join(keys(s.names), ","), color = :yellow)
    s.modules isa Dict && printstyled(io, " ", join(keys(s.modules), ","), color = :blue)
    println(io)
end

"""
scopehasmodule(s::Scope, mname::Symbol)::Bool

Checks whether the module `mname` has been `using`ed in `s`.
"""
scopehasmodule(s::Scope, mname::Symbol) = s.modules !== nothing && haskey(s.modules, mname)

"""
    addmoduletoscope!(s, m, [mname::Symbol])

Adds module `m` to the list of used modules in scope `s`.
"""
function addmoduletoscope!(s::Scope, m, mname::Symbol)
    if s.modules === nothing
        s.modules = Dict{Symbol,Any}()
    end
    s.modules[mname] = m
end
addmoduletoscope!(s::Scope, m::SymbolServer.ModuleStore) = addmoduletoscope!(s, m, m.name.name)
addmoduletoscope!(s::Scope, m::EXPR) =  CSTParser.defines_module(m) && addmoduletoscope!(s, scopeof(m), Symbol(valof(CSTParser.get_name(m))))
addmoduletoscope!(s::Scope, s1::Scope) = CSTParser.defines_module(s1.expr) && addmoduletoscope!(s, s1, Symbol(valof(CSTParser.get_name(s1.expr))))

getscopemodule(s::Scope, m::Symbol) = s.modules[m]

"""
    scopehasbinding(s::Scope, n::String)

Checks whether s has a binding for variable named `n`.
"""
scopehasbinding(s::Scope, n::String) = haskey(s.names, n)

"""
    introduces_scope(x::EXPR, state)

Does this expression introduce a new scope?
"""
function introduces_scope(x::EXPR, state)
    #TODO: remove unused 2nd argument.
    if typof(x) === CSTParser.BinaryOpCall
        if kindof(x[2]) === CSTParser.Tokens.EQ && CSTParser.is_func_call(x[1])
            return true
        elseif kindof(x[2]) === CSTParser.Tokens.EQ && typof(x[1]) === CSTParser.Curly
            return true
        elseif kindof(x[2]) === CSTParser.Tokens.ANON_FUNC
            return true
        else
            return false
        end
    elseif typof(x) === CSTParser.WhereOpCall
        # unless in func def signature
        return !_in_func_def(x)
    elseif typof(x) === CSTParser.TupleH && length(x) > 2 && typof(x[1]) === CSTParser.PUNCTUATION && is_assignment(x[2])
        return true
    elseif typof(x) === CSTParser.FunctionDef ||
            typof(x) === CSTParser.Macro ||
            typof(x) === CSTParser.For ||
            typof(x) === CSTParser.While ||
            typof(x) === CSTParser.Let ||
            typof(x) === CSTParser.Generator || # and Flatten?
            typof(x) === CSTParser.Try ||
            typof(x) === CSTParser.Do ||
            typof(x) === CSTParser.ModuleH ||
            typof(x) === CSTParser.BareModule ||
            typof(x) === CSTParser.Abstract ||
            typof(x) === CSTParser.Primitive ||
            typof(x) === CSTParser.Mutable ||
            typof(x) === CSTParser.Struct
        return true
    end
    return false
end

hasscope(x::EXPR) = hasmeta(x) && hasscope(x.meta)
scopeof(x) = nothing
scopeof(x::EXPR) = scopeof(x.meta)
CSTParser.parentof(s::Scope) = s.parent

function setscope!(x::EXPR, s)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.scope = s
end

"""
    scopes(x::EXPR, state)

Called when traversing the syntax tree and handles the association of
scopes with expressions. On the first pass this will add scopes as
necessary, on following passes it empties it.
"""
function scopes(x::EXPR, state)
    clear_scope(x)
    if scopeof(x) === nothing && introduces_scope(x, state)
        setscope!(x, Scope(x))
    end
    s0 = state.scope
    if typof(x) === FileH
        setscope!(x, state.scope)
    elseif scopeof(x) isa Scope
        if CSTParser.defines_function(x) || CSTParser.defines_macro(x)
            state.delayed = true # Allow delayed resolution
        end
        scopeof(x) != s0 && setparent!(scopeof(x), s0)
        state.scope = scopeof(x)
        if typof(x) === ModuleH # Add default modules to a new module
            state.scope.modules = Dict{Symbol,Any}()
            state.scope.modules[:Base] = getsymbolserver(state.server)[:Base]
            state.scope.modules[:Core] = getsymbolserver(state.server)[:Core]
        elseif typof(x) === BareModule
            state.scope.modules = Dict{String,Any}()
            state.scope.modules[:Core] = getsymbolserver(state.server)[:Core]
        end
        if (typof(x) === CSTParser.ModuleH || typof(x) === CSTParser.BareModule) && bindingof(x) !== nothing # Add reference to out of scope binding (i.e. itself)
            # state.scope.names[bindingof(x).name] = bindingof(x)
            # TODO: move this to the binding stage
            add_binding(x, state)
        elseif typof(x) === CSTParser.Flatten && typof(x[1]) === CSTParser.Generator && length(x[1]) > 0 && typof(x[1][1]) === CSTParser.Generator
            setscope!(x[1][1], nothing)
        end
    end
    return s0
end


text:
# packing and unpacking of Graph related types

using IncrementalInference
using DistributedFactorGraphs

using Test

dfg = GraphsDFG{SolverParams}(params=SolverParams())
setSerializationModule!(dfg, Main)

@testset "hard-coded test of PackedPrior to Prior" begin

pt = PackedPrior("KDE:100:[1.5015]:[-98.8276 -101.803 -98.1296 -96.1897 -99.3076 -99.1881 -101.721 -101.332 -100.431 -97.7293 -96.7652 -99.3806 -95.5593 -104.237 -94.9318 -101.691 -102.255 -98.9559 -99.3386 -99.2361 -102.483 -102.896 -97.0244 -98.9643 -99.4457 -101.298 -103.263 -2.75251 5.14065 0.327863 3.60042 -0.604114 -0.0564047 -0.804898 3.05469 1.4974 1.34657 2.22745 4.78117 1.89485 -0.48091 6.63068 0.63013 -3.11422 1.73705 5.22904 -1.73223 2.47453 1.10584 -0.0179944 3.65585 4.50016 -1.95042 98.2664 98.9983 103.748 100.789 98.4127 101.397 104.364 102.125 96.3685 103.59 99.0581 100.913 101.461 105.211 103.513 99.3325 101.201 98.05 103.508 99.9785 104.624 100.202 100.258 101.579 96.6931 95.4181 299.02 296.804 301.322 298.127 299.578 298.36 296.339 300.156 299.641 297.731 299.822 296.941 295.857 299.482 302.531 301.875 302.192 301.999 300.634 294.084 300.44]")

tt = convert(BallTreeDensity, pt.Z)
@test isa(tt, BallTreeDensity)

F = Prior

upt = convert(F, pt)

@test isa(upt, Prior)

  # TODO add more tests
end

fg = initfg()
N=100

doors = zeros(1,4)
doors[1,1:4] = [-100.0;0.0;100.0;300.0]
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
bws2 = Array{Float64,2}(undef, length(bws),1)
bws2[:,1] = bws[:]
doors2 = getPoints(pd);
v1 = addVariable!(fg,:x1, ContinuousScalar,N=N)
f1  = addFactor!(fg,[:x1], Prior(kde!(doors2, bws2[:])))

v2 = addVariable!(fg,:x2, ContinuousScalar, N=N)
lc = LinearConditional(  Normal(50.0, 2.0) )
f2 = addFactor!(fg, [:x1; :x2], lc)

@testset "Testing conversion to packed function node data structure and back" begin

topack = solverData(f1)
dd = convert(PackedFunctionNodeData{PackedPrior},topack)
upd = convert(FunctionNodeData{CommonConvWrapper{Prior}}, dd)

@test compare(topack, upd)

topack = solverData(f2)
dd =  convert(IncrementalInference.PackedFunctionNodeData{PackedLinearConditional},topack)
upd = convert(IncrementalInference.FunctionNodeData{CommonConvWrapper{LinearConditional}}, dd)

@test compare(topack, upd)

end

@testset "Testing conversion to packed variable node data structure and back" begin

dat = solverData(getVariable(fg,:x1))

# dat.BayesNetVertID

pd = pack(dfg, dat)
unpckd = unpack(dfg, pd)

@test compareFields(dat, unpckd, skip=[:softtype])
@test compareFields(dat.softtype, unpckd.softtype)
@test isa(dat.softtype, ContinuousScalar)
@test isa(unpckd.softtype, ContinuousScalar)

end

#


text:
using Documenter
using Pidfile

makedocs(
    modules = [Pidfile],
    format = :html,
    sitename = "Pidfile.jl",
    pages = [
        "Home" => "index.md",
    ],
)


text:
using JuMP, BARON, Base.Test

m = Model(solver=BaronSolver())
ub = [6,4]
@variable(m, 0 ≤ x[i=1:2] ≤ ub[i])

@NLconstraint(m, x[1]*x[2] ≤ 4)

@NLobjective(m, Min, -x[1] - x[2])

solve(m)

@test isapprox(getvalue(x[1]), 6, rtol=1e-6)
@test isapprox(getvalue(x[2]), 2/3, rtol=1e-6)
@test isapprox(getobjectivevalue(m), -20/3, rtol=1e-6)


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using SparseArrays
Random.seed!(233)

function mpi_create_matrix(indices,values,ilower,iupper)
    mpi_create_matrix_ = load_op_and_grad("./build/libMPITensor","mpi_create_matrix", multiple=true)
    indices,values,ilower,iupper = convert_to_tensor(Any[indices,values,ilower,iupper], [Int64,Float64,Int64,Int64])
    mpi_create_matrix_(indices,values,ilower,iupper)
end

function mpi_get_matrix(rows,ncols,cols,ilower,iupper,values, N)
    mpi_get_matrix_ = load_op_and_grad("./build/libMPITensor","mpi_get_matrix", multiple=true)
    rows,ncols,cols,ilower_,iupper_,values = convert_to_tensor(Any[rows,ncols,cols,ilower,iupper,values], [Int32,Int32,Int32,Int64,Int64,Float64])
    indices, vals = mpi_get_matrix_(rows,ncols,cols,ilower_,iupper_,values)
    SparseTensor(tf.SparseTensor(indices, vals, (iupper-ilower+1, N)), false)
end

function mpi_tensor_solve(rows,ncols,cols,values,rhs,ilower,iupper,solver = "BoomerAMG",printlevel = 2)
    mpi_tensor_solve_ = load_op_and_grad("./build/libMPITensor","mpi_tensor_solve")
    rows,ncols,cols,values,rhs,ilower,iupper,printlevel = convert_to_tensor(Any[rows,ncols,cols,values,rhs,ilower,iupper,printlevel], [Int32,Int32,Int32,Float64,Float64,Int64,Int64,Int64])
    mpi_tensor_solve_(rows,ncols,cols,values,rhs,ilower,iupper,solver,printlevel)
end


text:
function stairs3(x)
    s1 = 0.125
    s2 = 0.375
    s3 = 0.625
    s4 = 0.875
    a0 = 0.00
    a1 = 0.25
    a2 = 0.50
    a3 = 0.75
    y = IfElse.ifelse(x[1] < s1, a0,
            IfElse.ifelse(x[1] < s2, a1,
                IfElse.ifelse(x[1] < s3, a2,
                    IfElse.ifelse(x[1] < s4, a3, a0))))
    return y
end

function d_stairs3(x)
    return 0.0
end

stairs3_2D_lc = Environment{R2, LinearizedCone}(stairs3, d_stairs3)

# X = Vector(0:0.01:1.5)
# plot(X, stairs3steps.(X))

# Kernel
function f(x, center, radius, sharpness)
    return sharpness*(-((x-center)/radius)^2 + 1.0)
end

function smoothed_stairs(x)
    c = [0.0, 0.5, 1.0, 1.5, 2.0]
    a = [0.00, 0.25, 0.50, 0.75, 0.00]
    r = 0.25
    s = 3
    # Kernel vector
    v = [f(x, c[i], r, s) for i=1:5]
    # softmax
    ev = exp.(v)
    w = ev/sum(ev)
    # altitude selection
    return w'*a
end

# X = Vector(0:0.01:3)
# plot(X, smoothed_stairs.(X))