text:
import Base: copy!, eltype, getindex, isvalid, length, ndims, setindex!, show, size

#### Interpolation of evenly-spaced data ####

# This implements "generalized interpolation." See, for example,
# P. Thevenaz, T. Blu, and M. Unser (2000). "Interpolation Revisited."
# IEEE Transactions on Medical Imaging, 19: 739-758.

# This type manages temporary storage needed for efficient
# interpolation on a grid, so that once created it is possible to
# calculate many interpolated values without additional
# temporaries. It is designed to support both single- and multi-valued
# interpolation, and computation of the gradient with respect to the
# position of evaluation.
# An example of multi-valued interpolation is sub-pixel interpolation
# within an RGB image, and getting the 3 color values out.
type InterpGridCoefs{T, IT<:InterpType}
    coord1d::Vector{Vector{Int}}  # for 1-d positions
    coef1d::Vector{Vector{T}}     # for 1-d coefficients
    gcoef1d::Vector{Vector{T}}    # for 1-d coefficients of the gradient
    hcoef1d::Vector{Vector{T}}    # for 1-d coefficients of the hessian
    c1d::Vector{Vector{T}}        # temp space for switch btw value & gradient
    dims::Vector{Int}
    strides::Vector{Int}
    offset::Vector{Int}  # the default N-d offset when there is no wrapping
    offset_base::Int
    index::Vector{Int}   # the alternative N-d index, needed when wrapping
    wrap::Bool           # which of these two indices to use
    valid::Bool          # some BCs return NaN when wrapping
    coef::Vector{T}      # the N-d coefficient of each neighbor (val & gradient)
end

#### InterpGrid ####
# This is the high-level interface. Use this for scalar-valued
# interpolation, and for calculation of the gradient with respect to
# the location of the evaluation point.
# For efficient multi-valued interpolation, use the low-level
# interface below. (getindex and valgrad provide examples of how to use the
# low-level interface.)
abstract AbstractInterpGrid{T, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractArray{T,N}

type InterpGrid{T<:AbstractFloat, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractInterpGrid{T,N,BC,IT}
    coefs::Array{T,N}
    ic::InterpGridCoefs{T, IT}
    x::Vector{T}
    fillval::T  # used only for BCfill (if ever)
end
function InterpGrid{T<:AbstractFloat, N, BC<:BoundaryCondition, IT<:InterpType}(A::Array{T, N}, ::Type{BC}, ::Type{IT})
    if BC == BCfill
        error("Construct BCfill InterpGrids by supplying the fill value")
    end
    coefs = copy(A)
    interp_invert!(coefs, BC, IT, 1:N)
    ic = InterpGridCoefs(coefs, IT)
    x = zeros(T, N)
    InterpGrid{T, N, BC, IT}(coefs, ic, x, convert(T, NaN))
end
function InterpGrid{T<:AbstractFloat, N, BC<:Union{BCnil,BCnan,BCna}}(A::Array{T, N}, ::Type{BC}, ::Type{InterpCubic})
    # Cubic interpolation requires padding
    coefs=pad1(copy(A),0,1:N)
    interp_invert!(coefs, BC, InterpCubic, 1:N)
    ic = InterpGridCoefs(coefs, InterpCubic)
    x = zeros(T, N)
    InterpGrid{T, N, BC, InterpCubic}(coefs, ic, x, convert(T, NaN))
end
function InterpGrid{T<:AbstractFloat, N, IT<:InterpType}(A::Array{T, N}, f::Number, ::Type{IT})
    coefs = pad1(A, f, 1:N)
    interp_invert!(coefs, BCnearest, IT, 1:N)
    ic = InterpGridCoefs(coefs, IT)
    x = zeros(T, N)
    InterpGrid{T, N, BCfill, IT}(coefs, ic, x, convert(T, f))
end

setx{T}(G::InterpGrid{T,1}, x::Real) = @inbounds G.x[1] = x
function setx{T}(G::InterpGrid{T,2}, x::Real, y::Real)
    xG = G.x
    @inbounds xG[1] = x
    @inbounds xG[2] = y
end
function setx{T}(G::InterpGrid{T,3}, x::Real, y::Real, z::Real)
    xG = G.x
    @inbounds xG[1] = x
    @inbounds xG[2] = y
    @inbounds xG[3] = z
end
function setx{T,N}(G::InterpGrid{T,N}, x::Real...)
    if length(x) != N
        error("Incorrect number of dimensions supplied")
    end
    xG = G.x
    @inbounds for idim = 1:N
        xG[idim] = x[idim]
    end
end
# a version that corrects for the padding of BCfill types
setx{T}(G::InterpGrid{T,1,BCfill}, x::Real) = @inbounds G.x[1] = x+1
function setx{T}(G::InterpGrid{T,2,BCfill}, x::Real, y::Real)
    xG = G.x
    @inbounds xG[1] = x+1
    @inbounds xG[2] = y+1
end
function setx{T}(G::InterpGrid{T,3,BCfill}, x::Real, y::Real, z::Real)
    xG = G.x
    @inbounds xG[1] = x+1
    @inbounds xG[2] = y+1
    @inbounds xG[3] = z+1
end
function setx{T,N}(G::InterpGrid{T,N,BCfill}, x::Real...)
    if length(x) != N
        error("Incorrect number of dimensions supplied")
    end
    xG = G.x
    @inbounds for idim = 1:N
        xG[idim] = x[idim]+1
    end
end
# InterpCubic also needs to compensate for padding
setx{T,BC<:Union{BCna,BCnan,BCnil}}(G::InterpGrid{T,1,BC,InterpCubic}, x::Real) = @inbounds G.x[1] = x+1
function setx{T,BC<:Union{BCna,BCnan,BCnil}}(G::InterpGrid{T,2,BC,InterpCubic}, x::Real, y::Real)
    xG = G.x
    @inbounds xG[1] = x+1
    @inbounds xG[2] = y+1
end
function setx{T,BC<:Union{BCna,BCnan,BCnil}}(G::InterpGrid{T,3,BC,InterpCubic}, x::Real, y::Real, z::Real)
    xG = G.x
    @inbounds xG[1] = x+1
    @inbounds xG[2] = y+1
    @inbounds xG[3] = z+1
end
function setx{T,N,BC<:Union{BCna,BCnan,BCnil}}(G::InterpGrid{T,N,BC,InterpCubic}, x::Real...)
    if length(x) != N
        error("Incorrect number of dimensions supplied")
    end
    xG = G.x
    @inbounds for idim = 1:N
        xG[idim] = x[idim]+1
    end
end

## Evaluation at single points
function _getindex{T}(G::InterpGrid{T})
    set_position(G.ic, boundarycondition(G), false, false, G.x)
    interp(G.ic, G.coefs)
end
getindex{T}(G::InterpGrid{T,1}, x::Real) = (setx(G, x); _getindex(G))
function getindex{T}(G::InterpGrid{T,1}, x::Real, y::Real)
    if y != 1
        throw(BoundsError())
    end
    setx(G, x)
    _getindex(G)
end
getindex{T}(G::InterpGrid{T,2}, x::Real, y::Real) = (setx(G, x, y); _getindex(G))
getindex{T}(G::InterpGrid{T,3}, x::Real, y::Real, z::Real) = (setx(G, x, y, z); _getindex(G))
getindex{T,N}(G::InterpGrid{T,N}, x::Real...) = (setx(G, x...); _getindex(G))

function _valgrad{T,N}(g::AbstractVector{T}, G::InterpGrid{T,N})
    if length(g) != N
        error("Wrong number of components for the gradient")
    end
    ic = G.ic
    coefs = G.coefs
    set_position(ic, boundarycondition(G), true, false, G.x)
    val = interp(ic, coefs)
    for idim = 1:N
        set_gradient_coordinate(ic, idim)
        g[idim] = interp(ic, coefs)
    end
    return val
end
function _valgrad{T}(G::InterpGrid{T,1})
    ic = G.ic
    coefs = G.coefs
    set_position(ic, boundarycondition(G), true, false, G.x)
    val = interp(ic, coefs)
    set_gradient_coordinate(ic, 1)
    g = interp(ic, coefs)
    return val, g
end
function valgrad{T}(G::AbstractInterpGrid{T,1}, x::Real)
    setx(G, x)
    _valgrad(G)
end
function valgrad{T}(G::AbstractInterpGrid{T,2}, x::Real, y::Real)
    g = Array(T, 2)
    val = valgrad(g, G, x, y)
    return val, g
end
function valgrad{T}(G::AbstractInterpGrid{T,3}, x::Real, y::Real, z::Real)
    g = Array(T, 3)
    val = valgrad(g, G, x, y, z)
    return val, g
end
function valgrad{T}(G::AbstractInterpGrid{T}, x::Real...)
    g = Array(T, length(x))
    val = valgrad(g, G, x...)
    return val, g
end
valgrad{T}(g::AbstractVector{T}, G::InterpGrid{T}, x::Real) = (setx(G, x); _valgrad(g, G))
valgrad{T}(g::AbstractVector{T}, G::InterpGrid{T}, x::Real, y::Real) = (setx(G, x, y); _valgrad(g, G))
valgrad{T}(g::AbstractVector{T}, G::InterpGrid{T}, x::Real, y::Real, z::Real) = (setx(G, x, y, z); _valgrad(g, G))
valgrad{T}(g::AbstractVector{T}, G::InterpGrid{T}, x::Real...) = (setx(G, x...); _valgrad(g, G))

function _valgradhess{T,N}(g::AbstractVector{T}, h::AbstractMatrix{T}, G::InterpGrid{T,N})
    if length(g) != N
        error("Wrong number of components for the gradient")
    end
    if size(h,1) != N || size(h,2) != N
        error("Wrong number of components for the hessian")
    end
    ic = G.ic
    coefs = G.coefs
    set_position(ic, boundarycondition(G), true, true, G.x)
    val = interp(ic, coefs)
    for idim = 1:N
        set_gradient_coordinate(ic, idim)
        g[idim] = interp(ic, coefs)
        # diagonal elements of the hessian
        set_hessian_coordinate(ic, idim, idim)
        h[idim,idim] = interp(ic, coefs)
        # off-diagonal element. exploit the symmetry!
        for jdim = idim+1:N
            set_hessian_coordinate(ic, idim, jdim)
            h[idim,jdim] = interp(ic, coefs)
            h[jdim,idim] = h[idim,jdim]
        end
    end
    return val
end
function _valgradhess{T}(G::InterpGrid{T,1})
    ic = G.ic
    coefs = G.coefs
    set_position(ic, boundarycondition(G), true, true, G.x)
    val = interp(ic, coefs)
    set_gradient_coordinate(ic, 1)
    g = interp(ic, coefs)
    set_hessian_coordinate(ic, 1, 1)
    h = interp(ic, coefs)
    return val, g, h
end
function valgradhess{T}(G::AbstractInterpGrid{T,1}, x::Real)
    setx(G, x)
    _valgradhess(G)
end
function valgradhess{T}(G::AbstractInterpGrid{T,2}, x::Real, y::Real)
    g = Array(T, 2)
    h = Array(T, 2, 2)
    val = valgradhess(g, h, G, x, y)
    return val, g, h
end
function valgradhess{T}(G::AbstractInterpGrid{T,3}, x::Real, y::Real, z::Real)
    g = Array(T, 3)
    h = Array(T, 3, 3)
    val = valgradhess(g, h, G, x, y, z)
    return val, g, h
end
function valgradhess{T}(G::AbstractInterpGrid{T}, x::Real...)
    g = Array(T, length(x))
    h = Array(T, length(x), length(x))
    val = valgradhess(g, h, G, x...)
    return val, g, h
end
valgradhess{T}(g::Vector{T}, h::Matrix{T}, G::InterpGrid{T, 1}, x::Real) = (setx(G, x); _valgradhess(g, h, G))
valgradhess{T}(g::Vector{T}, h::Matrix{T}, G::InterpGrid{T, 2}, x::Real, y::Real) = (setx(G, x, y); _valgradhess(g, h, G))
valgradhess{T}(g::Vector{T}, h::Matrix{T}, G::InterpGrid{T, 3}, x::Real, y::Real, z::Real) = (setx(G, x, y, z); _valgradhess(g, h, G))
function valgradhess{T}(g::Vector{T}, h::Matrix{T}, G::InterpGrid{T}, x::Real...)
    setx(G, x...)
    _valgradhess(g, h, G)
end

## Vectorized evaluation at multiple points
function getindex{T,R<:Number}(G::AbstractInterpGrid{T,1}, x::AbstractVector{R})
    n = length(x)
    v = Array(T, n)
    for i = 1:n
        v[i] = getindex(G,x[i])
    end
    v
end
getindex{T,N,R<:Number}(G::AbstractInterpGrid{T,N}, x::AbstractVector{R}) = error("Linear indexing not supported")
function getindex{T,R<:Number}(G::AbstractInterpGrid{T,2}, x::AbstractVector{R}, y::AbstractVector{R})
    nx, ny = length(x), length(y)
    v = Array(T, nx, ny)
    for i = 1:nx
        for j = 1:ny
            v[i,j] = getindex(G, x[i], y[j])
        end
    end
    v
end
function getindex{T,N,R<:Number}(G::AbstractInterpGrid{T,N}, x::AbstractVector{R}, xrest::AbstractVector{R}...)
    if length(xrest) != N-1
        error("Dimensionality mismatch")
    end
    nx = length(x)
    nrest = [length(y) for y in xrest]
    v = Array(T, nx, nrest...)
    for c in Counter(nrest)
        for i = 1:nx
            setx(G, x[i], ntuple(j->xrest[j][c[j]], N-1)...)  # FIXME performance?? May not matter...
            v[i,c...] = _getindex(G)
        end
    end
    v
end

# This should work for the convert in the InterpIrregular constructor below,
# but for some reason it does not
#convert{T,S}(::Type{Array{T,1}}, r::Ranges{S}) = [convert(T, x) for x in r]

#### Non-uniform grid interpolation ####

# Currently supports only 1d, nearest-neighbor or linear
# Consequently, the internal representation may change in the future
# BCperiodic and BCreflect not supported
type InterpIrregular{T<:Number, S, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractInterpGrid{S,N,BC,IT}
    grid::Vector{Vector{T}}
    coefs::Array{S,N}
    x::Vector{T}
    fillval::S  # used only for BCfill (if ever)
end
InterpIrregular{T<:Number, BC<:BoundaryCondition, IT<:Union{InterpForward,InterpBackward,InterpNearest,InterpLinear}}(grid::Vector{T}, A::AbstractArray, ::Type{BC}, ::Type{IT}) =
    InterpIrregular(Vector{T}[grid], A, BC, IT) # special 1d syntax
InterpIrregular{T<:Number, BC<:BoundaryCondition, IT<:Union{InterpForward,InterpBackward,InterpNearest,InterpLinear}}(grid::(@compat Tuple{Vararg{Vector{T}}}), A::AbstractArray, ::Type{BC}, ::Type{IT}) =
    InterpIrregular(Vector{T}[grid...], A, BC, IT)
function InterpIrregular{T<:Number, S, N, BC<:BoundaryCondition, IT<:Union{InterpForward,InterpBackward,InterpNearest,InterpLinear}}(grid::Vector{Vector{T}}, A::AbstractArray{S, N}, ::Type{BC}, ::Type{IT})
    if length(grid) != 1
        error("Sorry, for now only 1d is supported")
    end
    if BC == BCreflect || BC == BCperiodic
        error("Sorry, reflecting or periodic boundary conditions not yet supported")
    end
    for i = 1:length(grid)
        if !issorted(grid[i])
            error("Coordinates must be supplied in increasing order")
        end
    end
    grid = copy(grid)
    coefs = Array(S, size(A))
    copy!(coefs, A)
    x = zeros(T, N)
    InterpIrregular{T, S, N, BC, IT}(grid, coefs, x, convert(S, NaN))
end
function InterpIrregular{IT<:InterpType}(grid, A::Array, f::Number, ::Type{IT})
    iu = InterpIrregular(grid, A, BCfill, IT)
    iu.fillval = f
    iu
end

function _getindexii{T,S,BC<:Union{BCfill,BCna,BCnan}}(G::InterpIrregular{T,S,1,BC}, x::Number)
    g = G.grid[1]
    i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
    (i == 1 || i == length(g)+1) ? G.fillval : _interpu(x, g, i, G.coefs, interptype(G))
end
function _getindexii{T,S}(G::InterpIrregular{T,S,1,BCnil}, x::Number)
    g = G.grid[1]
    i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
    (i == 1 || i == length(g)+1) ? throw(BoundsError()) : _interpu(x, g, i, G.coefs, interptype(G))
end
function _getindexii{T,S}(G::InterpIrregular{T,S,1,BCnearest}, x::Number)
    g = G.grid[1]
    i = (x == g[1]) ? 2 : searchsortedfirst(g, x)
    i == 1 ? G.coefs[1] : i == length(g)+1 ? G.coefs[end] : _interpu(x, g, i, G.coefs, interptype(G))
end
# This next is necessary for precedence
getindex(G::InterpIrregular, x::Real) = _getindexii(G, x)
getindex(G::InterpIrregular, x::Number) = _getindexii(G, x)

_interpu(x, g, i, coefs, ::Type{InterpForward}) = coefs[i]
_interpu(x, g, i, coefs, ::Type{InterpBackward}) = coefs[i-1]
_interpu(x, g, i, coefs, ::Type{InterpNearest}) = (x-g[i-1] < g[i]-x) ? coefs[i-1] : coefs[i]
function _interpu(x, g, i, coefs, ::Type{InterpLinear})
    f = (x-g[i-1])/(g[i]-g[i-1])
    (1-f)*coefs[i-1] + f*coefs[i]
end

# Low-level interpolation interface

# Low-level interpolation. The process is split into two parts:
#   1. Identify the neighbors that will contribute and calculate the
#      interpolation coefficients
#   2. Use the neighbor indices and coefficients to produce
#      interpolated values
# Splitting it this way helps support efficient multi-valued
# interpolation and gradient computation.
#
# set_position handles step 1 of this process.
function set_position{T,BC<:BoundaryCondition,IT<:InterpType}(ic::InterpGridCoefs{T,IT}, ::Type{BC}, calc_grad::Bool, calc_hess::Bool, x::Vector{T})
    N = ndims(ic)
    if length(x) != N
        error("Dimensionality mismatch")
    end
    valid = true
    for idim = 1:N
        if !isvalid(BC, IT, x[idim], ic.dims[idim])
            valid = false
            break
        end
    end
    if !valid && BC == BCnil
        error("Invalid location")
    end
    ic.valid = valid
    wrap = false
    if valid
        ib::Int = 0
        for idim = 1:N
            ix::Int, dx::T, twrap::Bool = interp_coords_1d(ic.coord1d[idim], BC, IT, x[idim], ic.dims[idim])
            wrap = wrap | twrap
            ib += (ix-1)*ic.strides[idim]
            interp_coefs_1d(ic.coef1d[idim], BC, IT, dx)
            if calc_grad
                interp_gcoefs_1d(ic.gcoef1d[idim], BC, IT, dx)
            end
            if calc_hess
                interp_hcoefs_1d(ic.hcoef1d[idim], BC, IT, dx)
            end
        end
        if wrap
            interp_index_coef(ic.index, ic.coef, ic.coord1d, ic.coef1d, ic.strides)
        else
            ic.offset_base = ib
            interp_coef(ic.coef, ic.coef1d)
        end
    end
    ic.wrap = wrap
end
set_position{T,BC<:BoundaryCondition,IT<:InterpType}(ic::InterpGridCoefs{T,IT}, ::Type{BC}, calc_grad::Bool, x::Vector{T}) = set_position(ic, BC, calc_grad, false, x)

# Once you're done calculating interpolated _values_, call this function
# for each component of the _gradient_.
# It effectively re-does step #1 above for the gradient component,
# although it skips over the calculation of the neighbor indices. You
# then call "interp" below to evaluate the gradient, just as you to
# evaluate the value. Put each component in a loop to get all the
# gradient components.
# Naturally, if you're interpolating a multi-valued function, it's far
# better to evaluate all values, then all first components of the
# gradient, then all second components of the gradient, etc, than it
# is to recalculate the coefficients for each of the multiple values.
function set_gradient_coordinate(ic::InterpGridCoefs, i::Int)
    if !ic.valid
        return
    end
    N = ndims(ic)
    if i < 1 || i > N
        error("Wrong dimension index")
    end
    for idim = 1:N
        if idim == i
            ic.c1d[idim] = ic.gcoef1d[idim]
        else
            ic.c1d[idim] = ic.coef1d[idim]
        end
    end
    interp_coef(ic.coef, ic.c1d)
end
function set_hessian_coordinate(ic::InterpGridCoefs, i1::Int, i2::Int)
    if !ic.valid
        return
    end
    N = ndims(ic)
    if !(1 <= i1 <= N && 1 <= i2 <= N)
        error("Wrong dimension index")
    end
    for idim = 1:N
        if i1 == i2
            if idim == i1
                ic.c1d[idim] = ic.hcoef1d[idim]
            else
                ic.c1d[idim] = ic.coef1d[idim]
            end
        else
            if idim == i1 || idim == i2
                ic.c1d[idim] = ic.gcoef1d[idim]
            else
                ic.c1d[idim] = ic.coef1d[idim]
            end
        end
    end
    interp_coef(ic.coef, ic.c1d)
end

# "interp" evaluates the interpolation. Call "set_position" first. If
# you want to also evaluate the gradient, call
# "set_gradient_coordinate" followed by "interp" for each component.
function interp{T}(ic::InterpGridCoefs{T}, A::AbstractArray, index::Int)
    if !ic.valid
        return convert(T, NaN)
    end
    coef = ic.coef
    if ic.wrap
        offset = ic.index
        index -= 1
    else
        offset = ic.offset
        index += ic.offset_base
    end
    val = zero(T)
    for i = 1:length(coef)
        @inbounds c = coef[i]
        @inbounds val += c == zero(T) ? zero(T) : c*convert(T,A[offset[i]+index])
    end
    return convert(T, val)
end
interp(ic::InterpGridCoefs, A::AbstractArray) = interp(ic, A, 1)

# Change the array dimensions
function set_size(ic::InterpGridCoefs, dims, strides)
    N = ndims(ic)
    if length(strides) != N
        error("Strides do not have the correct dimensionality")
    end
    if length(dims) != N
        error("Dimensions do not have the correct dimensionality")
    end
    d = [dims...]
    s = [strides...]
    for idim = 1:N
        ic.dims[idim] = d[idim]
        ic.strides[idim] = s[idim]
        interp_coords_1d(ic.coord1d[idim], interptype(ic))
    end
    interp_index(ic.offset, ic.coord1d, s)
end
function set_size(ic::InterpGridCoefs, dims)
    N = length(dims)
    s = Array(Int, N)
    s[1] = 1
    for i = 1:N-1
        s[i+1] = s[i]*dims[i]
    end
    set_size(ic, dims, s)
end

# strides is supplied separately to allow working with subarrays,
# and/or specific dimensions. For example, an RGB image might be
# interpolated with respect to the spatial dimensions but not the
# dimension that holds color.
function InterpGridCoefs{T,IT<:InterpType}(::Type{T}, ::Type{IT}, dims::Union{Dims,Vector{Int}}, strides::Union{Dims,Vector{Int}})
    N = length(strides)
    if length(dims) != N
        error("Length of dims and strides must match")
    end
    coord1d = Array(Vector{Int}, N)
    coef1d = Array(Vector{T}, N)
    gcoef1d = Array(Vector{T}, N)
    hcoef1d = Array(Vector{T}, N)
    c1d = Array(Vector{T}, N)
    l = npoints(IT)
    for idim = 1:N
        coord1d[idim] = Array(Int, l)
        coef1d[idim] = Array(T, l)
        gcoef1d[idim] = Array(T, l)
        hcoef1d[idim] = Array(T, l)
        # do not allocate entries for c1d
    end
    n_coef = l^N
    offset = Array(Int, n_coef)
    index = Array(Int, n_coef)
    coef = Array(T, n_coef)
    # Pre-calculate the default offset from the strides
    for idim = 1:N
        interp_coords_1d(coord1d[idim], IT)
    end
    interp_index(offset, coord1d, strides)
    InterpGridCoefs{T,IT}(coord1d,coef1d,gcoef1d,hcoef1d,c1d,[dims...],[strides...],offset,0,index,false,false,coef)
end
InterpGridCoefs{IT<:InterpType}(A::Array, ::Type{IT}) = InterpGridCoefs(eltype(A), IT, [size(A)...], [strides(A)...])

# Interpolation support routines

npoints(::Type{InterpForward}) = 1
npoints(::Type{InterpBackward}) = 1
npoints(::Type{InterpNearest}) = 1
npoints(::Type{InterpLinear}) = 2
npoints(::Type{InterpQuadratic}) = 3
npoints(::Type{InterpCubic}) = 4

# Test whether a given coordinate will yield an interpolated result
isvalid{BC<:BoundaryCondition,IT<:InterpType}(::Type{BC},::Type{IT}, x, len::Int) = true
isvalid{BC<:Union{BCnil,BCnan,BCna}}(::Type{BC}, ::Type{InterpForward}, x, len::Int) = x >= 0.0 && x <= len
isvalid{BC<:Union{BCnil,BCnan,BCna}}(::Type{BC}, ::Type{InterpBackward}, x, len::Int) = x >= 0.0 && x <= len
isvalid{BC<:Union{BCnil,BCnan,BCna}}(::Type{BC}, ::Type{InterpNearest}, x, len::Int) = x >= 0.5 && x-0.5 <= len
isvalid{BC<:Union{BCnil,BCnan,BCna}}(::Type{BC}, ::Type{InterpLinear}, x, len::Int) = x >= 1 && x <= len
isvalid{BC<:Union{BCnil,BCnan,BCna}}(::Type{BC}, ::Type{InterpQuadratic}, x, len::Int) = x >= 1 && x <= len
isvalid{BC<:Union{BCnil,BCnan,BCna}}(::Type{BC}, ::Type{InterpCubic}, x, len::Int) = 1 <= x <= len -1

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpForward})
    coord1d[1] = 0
end
# version for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpForward}, x, len::Int)
    ix = wrap(BC, round(Int, real(x)), len)
    coord1d[1] = ix
    return ix, zero(typeof(x)), false
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpBackward})
    coord1d[1] = 0
end
# version for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpBackward}, x, len::Int)
    ix = wrap(BC, round(Int, real(x)), len)
    coord1d[1] = ix
    return ix, zero(typeof(x)), false
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpNearest})
    coord1d[1] = 0
end
# version for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpNearest}, x, len::Int)
    ix = wrap(BC, round(Int, real(x)), len)
    coord1d[1] = ix
    return ix, zero(typeof(x)), false
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpLinear})
    coord1d[1] = 0
    coord1d[2] = 1
end
# version for indices
function interp_coords_1d{T,BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpLinear}, x::T, len::Int)
    ifx = floor(Int, x)
    dx = x-ifx
    ifx -= x < convert(T, ifx)
    ix = wrap(BC, ifx, len)
    @inbounds coord1d[1] = ix
    iswrap = (ix == len && dx > 0) || (ix == 1 && ifx < 1)
    @inbounds coord1d[2] = iswrap ? ix : wrap(BC, ix+1, len)
    return ix, dx, iswrap
end
function interp_coords_1d{T}(coord1d::Vector{Int}, ::Type{BCreflect}, ::Type{InterpLinear}, x::T, len::Int)
    ix = mod(floor(Int, x)-1, 2*len)
    dx = x-floor(Int, x)
    if ix < len
        ix += 1
        ixp = ix+1
        if ixp > len
            ixp = len
        end
    else
        ix = 2*len-ix-1
        dx = one(T)-dx
        ixp = ix+1
        if ix == 0
            ix = 1
        end
    end
    @inbounds coord1d[1] = ix
    @inbounds coord1d[2] = wrap(BCreflect, ixp, len)
    iswrap = (ix == ixp && dx > 0)
    return ix, dx, iswrap
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpQuadratic})
    coord1d[1] = -1
    coord1d[2] = 0
    coord1d[3] = 1
end
# versions for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpQuadratic}, x, len::Int)
    ix = round(Int, real(x))
    dx = x-ix
    ix = wrap(BC, ix, len)
    coord1d[2] = ix
    iswrap = (ix == 1 || ix == len)
    if iswrap
        coord1d[1] = wrap(BC, ix-1, len)
        coord1d[3] = wrap(BC, ix+1, len)
    else
        coord1d[1] = ix-1
        coord1d[3] = ix+1
    end
    return ix, dx, iswrap
end
# for InterpQuadratic, several require special handling
# BCnil, Bnan, BCna: for 1 <= x <= 1.5, continue the quadratic centered at x=2
function interp_coords_1d{BC<:Union{BCnil,BCnan,BCna}}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpQuadratic}, x, len::Int)
    if x > 1.5 && x+0.5 < len
        ix = round(Int, real(x))
    elseif x <= 1.5
        ix = 2
    else
        ix = len-1
    end
    dx = x-ix
    coord1d[1] = ix-1
    coord1d[2] = ix
    coord1d[3] = ix+1
    return ix, dx, false
end
# BCnearest & BCfill: for 1 <= x <= 1.5, ensure the slope tends to 0 at x=1
function interp_coords_1d{T,BC<:Union{BCnearest,BCfill}}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpQuadratic}, x::T, len::Int)
    if x > 1.5 && x+0.5 < len
        ix = round(Int, real(x))
        coord1d[1] = ix-1
        coord1d[2] = ix
        coord1d[3] = ix+1
        iswrap = false
        dx = x-convert(T,ix)
    elseif x <= 1.5
        ix = 1
        coord1d[1] = 2
        coord1d[2] = 1
        coord1d[3] = 2
        iswrap = true
        if x < 1
            dx = zero(T)
        else
            dx = x-convert(T,ix)
        end
    else
        ix = len
        coord1d[1] = len-1
        coord1d[2] = len
        coord1d[3] = len-1
        iswrap = true
        if x > len
            dx = zero(T)
        else
            dx = x-convert(T,ix)
        end
    end
    return ix, dx, iswrap
end
function interp_coords_1d(coord1d::Vector{Int}, ::Type{BCreflect}, ::Type{InterpQuadratic}, x, len::Int)
    ix = round(Int, real(x))
    dx = x-ix
    ix = mod(ix-1, 2*len)
    if ix < len
        ix += 1
    else
        dx = -dx
        ix = 2*len-ix
    end
    coord1d[2] = ix
    iswrap = (ix == 1 || ix == len)
    if iswrap
        coord1d[1] = wrap(BCreflect, ix-1, len)
        coord1d[3] = wrap(BCreflect, ix+1, len)
    else
        coord1d[1] = ix-1
        coord1d[3] = ix+1
    end
    return ix, dx, iswrap
end

# version for offsets
function interp_coords_1d(coord1d::Vector{Int}, ::Type{InterpCubic})
    coord1d[1] = -1
    coord1d[2] = 0
    coord1d[3] = 1
    coord1d[4] = 2
end
#versions for indices
function interp_coords_1d{BC<:BoundaryCondition}(coord1d::Vector{Int}, ::Type{BC}, ::Type{InterpCubic}, x, len::Int)
    ix = trunc(Int, real(x))
    dx = x - ix
    coord1d[2] = ix
    # the outermost non-wrapping indices are 2 and len-2
    iswrap = isvalid(BC, ix, 2, len-1)
    if iswrap
        coord1d[1] = wrap(BC, ix-1, len)
        coord1d[3] = wrap(BC, ix+1, len)
        coord1d[4] = wrap(BC, ix+2, len)
    else
        coord1d[1] = ix-1
        coord1d[3] = ix+1
        coord1d[4] = ix+2
    end
    return ix, dx, iswrap
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpForward}, dx::T)
    coef1d[1] = one(T)
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpForward}, dx::T)
    coef1d[1] = zero(T)
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpBackward}, dx::T)
    coef1d[1] = one(T)
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpBackward}, dx::T)
    coef1d[1] = zero(T)
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpNearest}, dx::T)
    coef1d[1] = one(T)
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpNearest}, dx::T)
    coef1d[1] = zero(T)
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpLinear}, dx::T)
    coef1d[1] = 1.0-dx
    coef1d[2] = dx
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpLinear}, dx::T)
    coef1d[1] = -one(T)
    coef1d[2] = one(T)
end

function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = (dx-0.5)^2/2
    coef1d[2] = 0.75-dx^2
    coef1d[3] = (dx+0.5)^2/2
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = dx-0.5
    coef1d[2] = -2dx
    coef1d[3] = dx+0.5
end
function interp_hcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpQuadratic}, dx::T)
    coef1d[1] = one(T)
    coef1d[2] = -2*one(T)
    coef1d[3] = one(T)
end
function interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpCubic}, dx::T)
    coef1d[1] = (1-dx)^3/6
    coef1d[2] = 2/3-dx^2*(2-dx)/2
    coef1d[3] = 2/3-(dx-1)^2*(1+dx)/2
    coef1d[4] = dx^3/6
end
function interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpCubic}, dx::T)
    coef1d[1] = -(1-dx)^2/2
    coef1d[2] = 3dx^2/2-2dx
    coef1d[3] = 1/2+dx*(1-(3/2)*dx)
    coef1d[4] = dx^2/2
end
function interp_hcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Vector{T}, ::Type{BC}, ::Type{InterpCubic}, dx::T)
    coef1d[1] = 1-dx
    coef1d[2] = 3dx-2
    coef1d[3] = 1-3dx
    coef1d[4] = dx
end

eltype{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = T
ndims{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = N
boundarycondition{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = BC
interptype{T, N, BC, IT}(G::InterpGrid{T, N, BC, IT}) = IT
size{T,N}(G::InterpGrid{T,N,BCfill}) = ntuple(i->size(G.coefs,i)-2, N)
size{T,N}(G::InterpGrid{T,N,BCfill}, i::Integer) = size(G.coefs, i)-2
size{T,N}(G::InterpGrid{T,N}) = size(G.coefs)
size{T,N}(G::InterpGrid{T,N}, i::Integer) = size(G.coefs, i)

eltype{T, S, N, BC, IT}(G::InterpIrregular{T, S, N, BC, IT}) = S
ndims{T, S, N, BC, IT}(G::InterpIrregular{T, S, N, BC, IT}) = N
boundarycondition{T, S, N, BC, IT}(G::InterpIrregular{T, S, N, BC, IT}) = BC
interptype{T, S, N, BC, IT}(G::InterpIrregular{T, S, N, BC, IT}) = IT
size(G::InterpIrregular) = size(G.coefs)
size(G::InterpIrregular, i::Integer) = size(G.coefs, i)

eltype{T,IT<:InterpType}(ic::InterpGridCoefs{T,IT}) = T
interptype{IT<:InterpType,T}(ic::InterpGridCoefs{T,IT}) = IT
ndims(ic::InterpGridCoefs) = length(ic.coord1d)
show(io::IO, ic::InterpGridCoefs) = print(io, "InterpGridCoefs{", eltype(ic), ",", interptype(ic), "}")

# Generalized interpolation of higher order than InterpLinear requires
# inversion of the interpolation operator. See, for example, the
# Thevanez citation above.

# We need efficient 1d slices:
immutable Slice1D{T,N} <: DenseArray{T,1}
    A::Array{T,N}
    index::StepRange{Int,Int}
end
size(S::Slice1D) = (length(S.index),)
size(S::Slice1D, d) = d==1 ? length(S.index) : 1
length(S::Slice1D) = length(S.index)

# Note: unsafe, no bounds-checking implemented
getindex{T}(S::Slice1D{T}, i::Real) = S.A[first(S.index)+(i-1)*step(S.index)]
getindex{T}(S::Slice1D{T}, i::Real, j::Real) = S[i]  # assume j=1, because a bounds-check will make it too slow
setindex!{T}(S::Slice1D{T}, x, i::Real) = (setindex!(S.A, x, first(S.index)+(i-1)*step(S.index)); S)
setindex!{T}(S::Slice1D{T}, x, i::Real, j::Real) = setindex!(S, x, i)

function copy!(dst::Vector, src::Slice1D)
    n = length(src)
    length(dst) == n || throw(DimensionMismatch("Length $(length(dst)) of destination does not equal length $n of source"))
    p = src.A
    rng = src.index
    j = first(rng)
    s = step(rng)
    for i = 1:n
        @inbounds dst[i] = p[j]
        j += s
    end
    dst
end

const Q3inv = [7/8 1/4 -1/8; -1/8 5/4 -1/8; -1/8 1/4 7/8] # for handling "snippets" of size 3 along each dimension (InterpQuadratic)
# This works in place. If instead it allocated the output for you, then
# calling multiple times (e.g., to apply inversions for different
# interptypes along different dimensions) would result in unnecessary
# allocations.
interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpForward}, dimlist) = A
interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpBackward}, dimlist) = A
interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpNearest}, dimlist) = A
interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpLinear}, dimlist) = A
function interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpQuadratic}, dimlist)
    sizeA = [size(A)...]
    stridesA = [strides(A)...]
    for idim = dimlist
        n = size(A, idim)
        # Set up the tridiagonal system
        du = fill(convert(eltype(A), 1/8), n-1)
        d = fill(convert(eltype(A), 3/4), n)
        dl = copy(du)
        M = _interp_invert_matrix(BC, InterpQuadratic, dl, d, du)
        sizeA[idim] = 1  # don't iterate over the dimension we're solving on
        for cc in Counter(sizeA)
            rng = range(coords2lin(cc, stridesA), stridesA[idim], n)
            A_ldiv_B!(M, Slice1D(A, rng))
        end
        sizeA[idim] = size(A, idim)
    end
    A
end
function interp_invert!{BC<:BoundaryCondition}(A::Array, ::Type{BC}, ::Type{InterpCubic}, dimlist)
    sizeA = [size(A)...]
    stridesA = [strides(A)...]

    for idim = dimlist
        n = size(A,idim)
        # Set up tridiagonal system: (c_i-1 + 4 c_i + c_i+1) / 6 = f_i
        du = fill(convert(eltype(A), 1/6), n-1)
        d = fill(convert(eltype(A), 4/6), n)
        dl = copy(du)
        M = _interp_invert_matrix(BC, InterpCubic, dl, d, du)
        sizeA[idim] = 1
        for cc in Counter(sizeA)
            rng = range(coords2lin(cc, stridesA), stridesA[idim], n)
            A_ldiv_B!(M, Slice1D(A, rng))
        end
        sizeA[idim] = size(A,idim)
    end
    A
end
interp_invert!{BC<:Union{BoundaryCondition,Number}}(A::Array, ::Type{BC}, IT) = interp_invert!(A, BC, IT, 1:ndims(A))
interp_invert!{BC<:Union{BoundaryCondition,Number}}(A::Array, ::Type{BC}, IT, dimlist...) = interp_invert!(A, BC, IT, dimlist)

function _interp_invert_matrix{BC<:Union{BCnil,BCnan,BCna},T}(::Type{BC}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    # For these, the quadratic centered on x=2 is continued down to
    # 1 rather than terminating at 1.5
    n = length(d)
    d[1] = d[n] = 9/8
    dl[n-1] = du[1] = -1/4
    MT = lufact!(Tridiagonal(dl, d, du))
    # Woodbury correction to add 1/8 for row 1, col 3 and row n, col n-2
    U = zeros(T, n, 2)
    V = zeros(T, 2, n)
    C = zeros(T, 2, 2)
    C[1,1] = C[2,2] = 1/8
    U[1,1] = U[n,2] = 1
    V[1,3] = V[2,n-2] = 1
    M = Woodbury(MT, U, C, V)
end
function _interp_invert_matrix{T}(::Type{BCreflect}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    n = length(d)
    d[1] += 1/8
    d[n] += 1/8
    M = lufact!(Tridiagonal(dl, d, du))
end
function _interp_invert_matrix{T}(::Type{BCperiodic}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    n = length(d)
    MT = lufact!(Tridiagonal(dl, d, du))
    # Woodbury correction to wrap around
    U = zeros(T, n, 2)
    V = zeros(T, 2, n)
    C = zeros(T, 2, 2)
    C[1,1] = C[2,2] = 1/8
    U[1,1] = U[n,2] = 1
    V[1,n] = V[2,1] = 1
    M = Woodbury(MT, U, C, V)
end
function _interp_invert_matrix{T,BC<:Union{BCnearest,BCfill}}(::Type{BC}, ::Type{InterpQuadratic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    n = length(d)
    du[1] += 1/8
    dl[n-1] += 1/8
    M = lufact!(Tridiagonal(dl, d, du))
end

function _interp_invert_matrix{BC<:Union{BCnil,BCnan,BCna},T}(::Type{BC}, ::Type{InterpCubic}, dl::Vector{T}, d::Vector{T}, du::Vector{T})
    # The most common way to terminate Cubic B-spline interpolations at the edges is by setting the second derivative to 0
    # Doing this yields c_1 = f_1 and c_N = f_N at the low and high end respectively
    # The BC also adds the equations a_0 = 2a_1 - a_2 for a ghost point a_0 outside the lower end of the domain,
    # and similarly for the higher end. These equations are represented in the top and bottom rows.
    d[1:2] = d[end-1:end] = 1
    du[1] = dl[end] = convert(T,-2)
    du[2] = dl[end-1] = du[end] = dl[1] = zero(T)
    MT = lufact!(Tridiagonal(dl, d, du))
    # Since the edge equations adds off-diagonal elemetns, we need Woodbury correction
    n = length(d)
    U = zeros(T, n, 2)
    V = zeros(T, 2, n)
    C = zeros(T, 2, 2)
    C[1,1] = C[2,2] = one(T)
    U[1,1] = U[n,2] = one(T)
    V[1,3] = V[2,n-2] = one(T)
    Woodbury(MT, U, C, V)
end

# We have 3 functions for computing indices and coefficients:
#   interp_coef: fast, computes just the coefficients
#   interp_index_coef: fast, computes both indices and coefficients
#   interp_index: slow, computes just the indices
# The latter is rarely used (presumably, only upon InterpGridCoefs
# allocation), and hence should not be performance-critical, so the
# code is kept simple.
# The motivation for this organization is that in the "generic" case
# (see below), the overhead of looping is such that it's better to
# compute the indices and coefficients at the same time. However, when
# you don't need the indices, for fixed dimensions it's a significant
# savings (roughly 2x, at least for 3d) to just compute the
# coefficients.
function interp_index(index::Vector{Int}, coord1d::Vector{Vector{Int}}, s)
    N = length(coord1d)
    l = map(length, coord1d)
    i = 1
    for c in Counter(l)
        indx = s[1]*coord1d[1][c[1]]
        for idim = 2:N
            indx += s[idim]*coord1d[idim][c[idim]]
        end
        index[i] = indx
        i += 1
    end
end
# This is the main bottleneck, so it pays handsomely to specialize on
# dimension.
function interp_index_coef{T}(indices::Array{Int}, coefs::Array{T}, coord1d::Vector{Vector{Int}}, coef1d::Vector{Vector{T}}, strides)
    n_dims = length(coef1d)
    l = length(coef1d[1])
    if n_dims == 1
        copy!(coefs, coef1d[1])
        s = strides[1]
        c = coord1d[1]
        for i = 1:l
            indices[i] = (c[i]-1)*s + 1
        end
    elseif n_dims == 2
        s1 = strides[1]
        s2 = strides[2]
        cd1 = coord1d[1]
        cd2 = coord1d[2]
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        ind = 1
        for i2 = 1:l
            offset = (cd2[i2]-1)*s2
            p = cf2[i2]
            for i1 = 1:l
                indices[ind] = offset + (cd1[i1]-1)*s1 + 1
                coefs[ind] = p*cf1[i1]
                ind += 1
            end
        end
    elseif n_dims == 3
        s1 = strides[1]
        s2 = strides[2]
        s3 = strides[3]
        cd1 = coord1d[1]
        cd2 = coord1d[2]
        cd3 = coord1d[3]
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        ind = 1
        for i3 = 1:l
            offset3 = (cd3[i3]-1)*s3
            p3 = cf3[i3]
            for i2 = 1:l
                offset2 = offset3 + (cd2[i2]-1)*s2
                p2 = p3*cf2[i2]
                for i1 = 1:l
                    indices[ind] = offset2 + (cd1[i1]-1)*s1 + 1
                    coefs[ind] = p2*cf1[i1]
                    ind += 1
                end
            end
        end
    elseif n_dims == 4
        s1 = strides[1]
        s2 = strides[2]
        s3 = strides[3]
        s4 = strides[4]
        cd1 = coord1d[1]
        cd2 = coord1d[2]
        cd3 = coord1d[3]
        cd4 = coord1d[4]
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        cf4 = coef1d[4]
        ind = 1
        for i4 = 1:l
            offset4 = (cd4[i4]-1)*s4
            p4 = cf4[i4]
            for i3 = 1:l
                offset3 = offset4 + (cd3[i3]-1)*s3
                p3 = p4*cf3[i3]
                for i2 = 1:l
                    offset2 = offset3 + (cd2[i2]-1)*s2
                    p2 = p3*cf2[i2]
                    for i1 = 1:l
                        indices[ind] = offset2 + (cd1[i1]-1)*s1 + 1
                        coefs[ind] = p2*cf1[i1]
                        ind += 1
                    end
                end
            end
        end
    else
        interp_index_coef_generic(indices, coefs, coord1d, coef1d, strides)
    end
end
function interp_coef{T}(coefs::Array{T}, coef1d::Vector{Vector{T}})
    n_dims = length(coef1d)
    l = length(coef1d[1])
    if n_dims == 1
        copy!(coefs, coef1d[1])
    elseif n_dims == 2
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        ind = 1
        for i2 = 1:l
            p = cf2[i2]
            for i1 = 1:l
                coefs[ind] = p*cf1[i1]
                ind += 1
            end
        end
    elseif n_dims == 3
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        ind = 1
        for i3 = 1:l
            p3 = cf3[i3]
            for i2 = 1:l
                p2 = p3*cf2[i2]
                for i1 = 1:l
                    coefs[ind] = p2*cf1[i1]
                    ind += 1
                end
            end
        end
    elseif n_dims == 4
        cf1 = coef1d[1]
        cf2 = coef1d[2]
        cf3 = coef1d[3]
        cf4 = coef1d[4]
        ind = 1
        for i4 = 1:l
            p4 = cf4[i4]
            for i3 = 1:l
                p3 = p4*cf3[i3]
                for i2 = 1:l
                    p2 = p3*cf2[i2]
                    for i1 = 1:l
                        coefs[ind] = p2*cf1[i1]
                        ind += 1
                    end
                end
            end
        end
    else
        interp_coef_generic(coefs, coef1d)
    end
end
# Compared to the dimension-specialized versions, these next two will
# be slow due to having many more arrayref operations. But at least it
# avoids an extra loop compared to the naive implementation using
# Counter, so it could be a lot worse.
# Consider making something like make_arrayind_loop_nest?
function interp_index_coef_generic{T}(indices::Array{Int}, coefs::Array{T}, coord1d::Vector{Vector{Int}}, coef1d::Vector{Vector{T}}, strides)
    n_dims = length(coord1d)
    l = length(coord1d[1])
    c = fill(1, n_dims)
    offset = Array(Int, n_dims)
    p = Array(T, n_dims)
    offset[end] = (coord1d[end][1]-1)*strides[end]
    p[end] = coef1d[end][1]
    for idim = n_dims-1:-1:1
        offset[idim] = offset[idim+1] + (coord1d[idim][1]-1)*strides[idim]
        p[idim] = p[idim+1] * coef1d[idim][1]
    end
    ind = 1
    while c[n_dims] <= l
        indices[ind] = offset[1]+1
        coefs[ind] = p[1]
        ind += 1
        c[1] += 1
        thisc = c[1]
        if thisc <= l
            if n_dims > 1
                p[1] = p[2] * coef1d[1][thisc]
                offset[1] = offset[2] + (coord1d[1][thisc]-1)*strides[1]
            else
                p[1] = coef1d[1][thisc]
                offset[1] = (coord1d[1][thisc]-1)*strides[1]
            end
        else
            idim = 1
            while c[idim] > l && idim < n_dims
                c[idim] = 1
                c[idim+1] += 1
                idim += 1
            end
            if idim == n_dims
                thisc = c[idim]
                if thisc <= l
                    p[idim] = coef1d[idim][thisc]
                    offset[idim] = (coord1d[idim][thisc]-1)*strides[idim]
                    idim -= 1
                end
            end
            while idim > 0 && idim < n_dims
                thisc = c[idim]
                offset[idim] = offset[idim+1] + (coord1d[idim][thisc]-1)*strides[idim]
                p[idim] = p[idim+1] * coef1d[idim][thisc]
                idim -= 1
            end
        end
    end
end
function interp_coef_generic{T}(coefs::Array{T}, coef1d::Vector{Vector{T}})
    n_dims = length(coef1d)
    l = length(coef1d[1])
    c = ones(Int, n_dims)
    offset = Array(Int, n_dims)
    p = Array(T, n_dims)
    p[end] = coef1d[end][1]
    for idim = n_dims-1:-1:1
        p[idim] = p[idim+1] * coef1d[idim][1]
    end
    ind = 1
    while c[n_dims] <= l
        coefs[ind] = p[1]
        ind += 1
        c[1] += 1
        thisc = c[1]
        if thisc <= l
            if n_dims > 1
                p[1] = p[2] * coef1d[1][thisc]
            else
                p[1] = coef1d[1][thisc]
            end
        else
            idim = 1
            while c[idim] > l && idim < n_dims
                c[idim] = 1
                c[idim+1] += 1
                idim += 1
            end
            if idim == n_dims
                thisc = c[idim]
                if thisc <= l
                    p[idim] = coef1d[idim][thisc]
                    idim -= 1
                end
            end
            while idim > 0 && idim < n_dims
                thisc = c[idim]
                p[idim] = p[idim+1] * coef1d[idim][thisc]
                idim -= 1
            end
        end
    end
end

# Padding an array by 1 element on each boundary of chosen dimensions
# This is needed to support BCfill
function pad1_index(szt::Tuple, dimpad)
    szv = [szt...]
    for i in dimpad
        szv[i] = szv[i].+2
    end
    N = length(szt)
    ind = Array{Any}(N)
    for i in 1:N
        if szv[i] > szt[i]
            ind[i] = 2:szv[i]-1
        else
            ind[i] = 1:szv[i]
        end
    end
    return szv, ind
end
pad1{BC<:BoundaryCondition}(Ain::AbstractArray, ::Type{BC}, dimpad...) = copy(Ain)  # no padding needed for most types
pad1(Ain::AbstractArray, ::Type{BCfill}, dimpad...) = error("For BCfill, specify the fill value instead of the type BCfill")
function pad1(Ain::AbstractArray, f::Number, dimpad...)
    sz, ind = pad1_index(size(Ain), dimpad)
    A = fill(convert(eltype(Ain), f), sz...)
    A[ind...] = Ain
    return A
end


text:
#Signal sources
#-------------------------------------------------------------------------------

#==
===============================================================================#

#==
===============================================================================#
function VSource(sig::Sine; name)
	𝜔 = 2π*sig.f

	@parameters A f
	defaults = Dict(A=>sig.A, f=>sig.f)

	@named(p = Pin()); @named(n = Pin())
	@variables v(t) i(t)
	eqs = [
		v ~ p.v - n.v
		0 ~ p.i + n.i #in = -out
		i ~ -p.i #Positive current flows *out* of p terminal
		v ~ sig.A*sin(𝜔*t)
	]
	ODESystem(eqs, t, [v, i], keys(defaults), systems=[p, n]; defaults, name)
end

function VSource(V::Real; name)
	#Rename before re-using variables:
	_V = V
	@parameters V
	defaults = Dict(V=>_V)

	@named(p = Pin()); @named(n = Pin())
	@variables v(t) i(t)
	eqs = [
		v ~ p.v - n.v
		0 ~ p.i + n.i #in = -out
		i ~ -p.i #Positive current flows *out* of p terminal
		v ~ V
	]
	ODESystem(eqs, t, [v, i], keys(defaults), systems=[p, n]; defaults, name)
end

Ideal

#Last line

ext: jl


text:
#list of added packages

#routines for adding the packages
add Atom Juno
add DataFrames Revise Distributions StatsBase GLM Dates NLopt ForwardDiff Formatting
add DataStructures LinearAlgebra Gadfly Measures Formatting CodecZlib CodecLz4
add CSV NLopt CategoricalArrays Serialization Distributed Random
add StaticArrays VegaLite CuArrays
add https://github.com/clintonTE/Finometrics

Thank you for the fresh take on portfolio choice. Since smoothing rules usually imply substantial dissaving during periods of high marginal utility, how strong would this effect be for modern endowments?

Are colelges underwater? Real estate is down, college admissions down

285797
dxnes
adaptive_de_rand_1_bin_radiuslimited
xnes
de_rand_1_bin_radiuslimited
adaptive_de_rand_1_bin
generating_set_search
de_rand_1_bin
separable_nes
resampling_inheritance_memetic_search
probabilistic_descent
resampling_memetic_search
de_rand_2_bin_radiuslimited
de_rand_2_bin
random_search
simultaneous_perturbation_stochastic_approximation

---
September 24 Meeting
September 1 is deadline for JMP

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
include("utils.jl")

run_test(parse(Int, ARGS[1]), parse(Int, ARGS[2]))

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
#!/usr/bin/env julia

using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, Pose2D, PoseStamped, Vector3
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan
@rosimport mavs_control.msg: Control

@rosimport gazebo_msgs.msg: ModelState
@rosimport gazebo_msgs.srv: SetModelState, GetModelState, GetWorldProperties

rostypegen()
using geometry_msgs.msg
using std_srvs.srv
using nav_msgs.srv.GetPlan
using gazebo_msgs.msg
using gazebo_msgs.srv
using mavs_control.msg
using PyCall
@pyimport tf.transformations as tf
using VehicleModels # for Linear_Spline()

# TODO
# 1) check to see if model is paused
# 2) make a while loop that publishes the current message until there is a new one ready or the current one ran out -> in this case stop the vehicle
# 3) interpolate message in time, is there a while! new message

function callback(ctr_msg::Control, get_state::ServiceProxy{GetModelState}, set_state::ServiceProxy{SetModelState})
    # interpolate optimal controls
    sp_X=Linear_Spline(ctr_msg.t,ctr_msg.x);
    sp_Y=Linear_Spline(ctr_msg.t,ctr_msg.y);
    sp_PSI=Linear_Spline(ctr_msg.t,ctr_msg.psi);

    modelName = RobotOS.get_param("robotName")
    println(" service...")

    t=to_sec(get_rostime())
    t0=to_sec(get_rostime())
    while((t-t0) < 0.5)
        println(t-t0)
        # Get the current position of the Gazebo model
        gs = GetModelStateRequest()
        gs.model_name = modelName
        gs_r = get_state(gs)

        if !gs_r.success
            error(string(" calling /gazebo/get_model_state service: ", gs_r.status_message))
        end

        # Define the current time
        t = to_sec(get_rostime())

        # Define position to move robot
        vehPose = gs_r.pose  # use the current position
        vehPose.position.x = sp_X[t]
        vehPose.position.y = sp_Y[t]
        roll = 0; pitch = 0; yaw = sp_PSI[t];
        Q = tf.quaternion_from_euler(roll, pitch, yaw)
        vehPose.orientation.x = Q[1]
        vehPose.orientation.y = Q[2]
        vehPose.orientation.z = Q[3]
        vehPose.orientation.w = Q[4]

        # Define the robot state
        ms = ModelState()
        ms.model_name = modelName
        ms.pose = vehPose

        # Set the state of the Gazebo model
        ss = SetModelStateRequest()
        ss.model_state = ms
    #    println("Calling 'gazebo/set_model_state' service...")
        ss_r = set_state(ss)

        if !ss_r.success
            error(string(" calling /gazebo/set_model_state service: ", ss_r.status_message))
        end
    end
end

function main()
    init_node("rosjl_move_hmmwv")

    # Set up service to get Gazebo model state
    const get_state = ServiceProxy("/gazebo/get_model_state", GetModelState)
    println("Waiting for 'gazebo/get_model_state' service...")
    wait_for_service("gazebo/get_model_state")

    # Set up service to set Gazebo model state
    const set_state = ServiceProxy("/gazebo/set_model_state", SetModelState)
    println("Waiting for 'gazebo/set_model_state' service...")
    wait_for_service("gazebo/set_model_state")

    sub_control = Subscriber{Control}("/mavs/optimal_control", callback, (get_state, set_state), queue_size = 2)

    spin()
end

if ! isinteractive()
    main()
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestOntologyJSON
using Test

using Catlab
using SemanticFlowGraphs

# Concepts
##########

TestPres = Presentation(String)
A, B, C, A0, B0 = Ob(Monocl, "A", "B", "C", "A0", "B0")
I = munit(Monocl.Ob)
f, f0 = Hom("f", A, B), Hom("f0", A0, B0)
add_generators!(TestPres, [A, B, C, A0, B0])
add_generators!(TestPres, [SubOb(A0, A), SubOb(B0, B)])
add_generators!(TestPres, [f, f0])
add_generator!(TestPres, Hom("g", I, otimes(A,B)))
add_generator!(TestPres, SubHom(f0, f, SubOb(A0, A), SubOb(B0, B)))

concept(pairs...) = Dict("schema" => "concept", pairs...)

const docs = [
  concept(
    "kind" => "type",
    "id" => "A",
  ),
  concept(
    "kind" => "type",
    "id" => "B",
  ),
  concept(
    "kind" => "type",
    "id" => "C",
  ),
  concept(
    "kind" => "type",
    "id" => "A0",
    "is-a" => "A",
  ),
  concept(
    "kind" => "type",
    "id" => "B0",
    "is-a" => "B",
  ),
  concept(
    "kind" => "function",
    "id" => "f",
    "inputs" => [
      Dict("type" => "A"),
    ],
    "outputs" => [
      Dict("type" => "B"),
    ]
  ),
  concept(
    "kind" => "function",
    "id" => "f0",
    "is-a" => "f",
    "inputs" => [
      Dict("type" => "A0"),
    ],
    "outputs" => [
      Dict("type" => "B0"),
    ]
  ),
  concept(
    "kind" => "function",
    "id" => "g",
    "inputs" => [],
    "outputs" => [
      Dict("type" => "A"),
      Dict("type" => "B"),
    ]
  )
]

pres = presentation_from_json(docs)
@test generators(pres, Monocl.Ob) == generators(TestPres, Monocl.Ob)
@test generators(pres, Monocl.Hom) == generators(TestPres, Monocl.Hom)
@test generators(pres, Monocl.SubOb) == generators(TestPres, Monocl.SubOb)
@test generators(pres, Monocl.SubHom) == generators(TestPres, Monocl.SubHom)

# Annotations
#############

pres = Presentation(String)
A, B, C, D = Ob(Monocl, "A", "B", "C", "D")
f = Hom("f", A, B)
g = Hom("g", B, C)
h = Hom("h", D, D)
add_generators!(pres, [A, B, C, D, f, g])

note = annotation_from_json(Dict(
  "schema" => "annotation",
  "kind" => "type",
  "language" => "python",
  "package" => "mypkg",
  "id" => "a",
  "class" => "ClassA",
  "definition" => "A"
), pres)
@test note.name == AnnotationID("python", "mypkg", "a")
@test note.language == Dict(:class => "ClassA")
@test note.definition == A

note = annotation_from_json(Dict(
  "schema" => "annotation",
  "kind" => "function",
  "language" => "python",
  "package" => "mypkg",
  "id" => "a-do-composition",
  "class" => "ClassA",
  "method" => "do_composition",
  "definition" => [
    "otimes",
    ["compose", "f", "g" ],
    ["Hom", "h", "D", "D" ],
  ]
), pres)
@test note.name == AnnotationID("python", "mypkg", "a-do-composition")
@test note.language == Dict(:class => "ClassA", :method => "do_composition")
@test note.definition == otimes(compose(f, g), h)

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
# > julia -L test/test.jl

using UnitFloat

a = UFloat(.24)
b = UFloat(.12)

v = UFloat.(rand((10000, 1)))
w = UFloat.(rand((10000, 1)))


text:
export Decoder, add!, decode, decode!, get_source, get_source!, RankDeficiencyException

"""
    Decoder{Tv,Ti<:Integer,Tm<:AbstractMatrix{Tv}}

Inactivation decoder compatible with Raptor10 (rfc5053) and RaptorQ (rfc6330) codes.
"""
mutable struct Decoder{Tv,Ti<:Integer,Tm<:AbstractMatrix{Tv}}
    # Constraint matrix
    columns::Vector{Vector{Ti}} # columns[cpi] contains the indices of rows (rpi) neighboring column cpi
    dense::Tm # dense (inactivated) symbols are stored separately
    # Permutation vectors
    colperm::Vector{Ti} # colperm[ri] gives an index for a vector in sparse
    colperminv::Vector{Ti} # inverse of rowperm
    rowperm::Vector{Ti} # sparse[rowperm[ri]] is the SparseVector of the ri-th row
    rowperminv::Vector{Ti} # inverse of rowperm
    uperm::Vector{Ti} # maps ui to upi
    uperminv::Vector{Ti} # maps upi to ui
    # Scalar values
    num_symbols::Ti # number of source symbols
    num_decoded::Ti # denoted by i in the R10 spec.
    num_inactivated::Ti # denoted by u in the R10 spec.
    num_rowops::Ti # for performance tracking (not used by the decoder)
    # Row schedule
    rowpq::PriorityQueue{Ti,Tuple{Ti,Ti},Base.Order.ForwardOrdering}
    componentpq::PriorityQueue{Ti,Ti,Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}}
    components::IntDisjointSets{Ti}
end

"""

Create a decoder object from a matrix `A`, where each column of `A` corresponds to a constraint.
"""
function Decoder(A::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}; initial_inactivation_storage::Integer=32) where {Tv,Ti}
    k, n = size(A)
    count(iszero, nonzeros(A)) == 0 || throw(ArgumentError("Structural zeros are not supported, i.e., there may be no explicitly stored zeros in A"))

    # permutation vectors
    rowperm = collect(Ti, 1:n)
    rowperminv = collect(Ti, 1:n)
    colperm = collect(Ti, 1:k)
    colperminv = collect(Ti, 1:k)
    uperm = collect(Ti, 1:initial_inactivation_storage)
    uperminv = collect(Ti, 1:initial_inactivation_storage)

    # Index each constraint by the source symbols it neighbors
    rows = rowvals(A)
    columns = [Vector{Ti}() for _ in 1:k]
    for j in 1:n
        for i in nzrange(A, j)
            cpi = rows[i]
            push!(columns[cpi], j)
        end
    end

    # The degree of a constraint is equal to the number of non-zero entries it has, whereas the
    # vdegree is the number of non-zero entries corresponding to source symbols that are neither
    # decoded nor inactivated. Rows are prioritized first by vdegree and second by degree. Note
    # that the vdegree and degree is upper-bounded by k.
    rowpq = PriorityQueue{Ti, Tuple{Ti,Ti}}(Base.Order.Forward)
    componentpq = PriorityQueue{Ti,Ti}(Base.Order.Reverse)
    components = IntDisjointSets{Ti}(k)
    for j in 1:n
        Is = nzrange(A, j)
        degree = length(Is)
        vdegree = degree # degree consisting of only symbols that are neither decoded nor inactivated
        priority = (vdegree, degree)
        enqueue!(rowpq, j, priority)

        # Two constraints are refered to as neighbors if both constraints have non-zero entries
        # corresponding to the same source symbol. Constraints with vdegree 2 are prioritized by
        # the number of neighboring that also have vdegree 2.
        if degree == 2
            cpi, cpj = rows[Is]
            merge_components!(componentpq, components, cpi, cpj)
        end
    end

    # Storage for inactivated symbols
    dense = zeros(Tv, initial_inactivation_storage, n)

    Decoder{Tv,Ti,typeof(dense)}(
        columns, dense,
        colperm, colperminv, rowperm, rowperminv, uperm, uperminv,
        k, 0, 0, 0,
        rowpq, componentpq, components
    )
end

struct RankDeficiencyException <: Exception end

"""

Print the decoder state to stdout (used for debugging).
"""
function print_state(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC)
    n = size(A, 2)
    k = d.num_symbols
    println("### Sparse constraint matrix ###")
    for ri in 1:n
        rpi = d.rowperm[ri]
        s = "($ri, $rpi)"
        s *= reduce(*, [" " for _ in 1:30-length(s)])
        s *= "["
        print(s)
        for ci in 1:k
            if ci == d.num_decoded + 1
                print(" | ")
            end
            cpi = d.colperm[ci]
            coef = A[cpi, rpi]
            print(" $(Int(coef)) ")
            # if iszero(coef)
            #     print(" ")
            # else
            #     print(".")
            # end
            if ci == n - d.num_inactivated
                print(" | ")
            end
        end
        println("]")
    end
    println()
    println("### Dense constraint matrix ###")
    for ri in 1:n
        rpi = d.rowperm[ri]
        # setinactive!(d, rpi)
        s = "($ri, $rpi)"
        s *= reduce(*, [" " for _ in 1:30-length(s)])
        s *= "["
        print(s)
        for ci in 1:k
            cpi = d.colperm[ci]
            coef = getdense(d, rpi, cpi)
            print(" $(Int(coef)) ")
            # if iszero(coef)
            #     print(" ")
            # else
            #     print(".")
            # end
        end
        println("]")
    end
    println()
    # println("### Permutation vectors ###")
    # println(d.rowperm)
    # println(d.rowperminv)
    # println(d.colperm)
    # println(d.colperminv)
    # println(d.uperm)
    # println(d.uperminv)
    # println()
    return
end

## getting, setting, indexing for the dense submatrix ###

# the dense matrix stores inactivated entries and should be seen as
# covering the num_inactivated rightmost columns of the constraint
# matrix. whenever a column is inactivated it covers one additional
# column. the first column of the dense matrix corresponds to the
# first column to be inactivated, i.e., the rightmost column of the
# constraint matrix, and the second column of the dense matrix
# corresponds to second column to be inactivated, i.e., the
# (num_symbols-1)-th column, and so on. the below methods convert
# column indices of the constraint matrix (ci) to/from column indices
# of the dense matrix (ui).
@inline _ci2ui(d::Decoder, ci::Integer) = d.num_symbols-ci+1
@inline _ui2ci(d::Decoder, ui::Integer) = d.num_symbols-ui+1

"""

Set the coefficient of the dense matrix corresponding to element `[rpi, cpi]` to `v`.
"""
function setdense!(d::Decoder, rpi::Integer, cpi::Integer, v)
    ci = d.colperminv[cpi]
    ui = _ci2ui(d, ci)
    upi = d.uperm[ui]
    d.dense[upi, rpi] = v
end

"""

Return the element of the dense matrix corresponding to permuted row and column indices
`[rpi, cpi]`.
"""
function getdense(d::Decoder, rpi::Integer, cpi::Integer)
    ci = d.colperminv[cpi]
    ui = _ci2ui(d, ci)
    if ui > length(d.uperm)
        return zero(eltype(d.dense))
    end
    upi = d.uperm[ui]
    d.dense[upi, rpi]
end

"""

Expand the matrix storing inactivated symbols.
"""
function expand_dense!(d::Decoder)
    m = size(d.dense, 1)
    while m < d.num_inactivated
        m *= 2
    end
    n = size(d.dense, 2)
    if m == size(d.dense, 1) && n == size(d.dense, 2)
        return
    end
    d.dense = expand_dense!(d.dense, m, n)
    return
end

function expand_dense!(dense::AbstractMatrix, m::Integer, n::Integer)
    mp, np = size(dense)
    if m < mp
        throw(ArgumentError("can't expand a matrix with $mp rows into $m rows"))
    end
    if n < np
        throw(ArgumentError("can't expand a matrix with $np columns into $n columns"))
    end
    rv = zeros(eltype(dense), m, n)
    rv[1:mp, 1:np] .= dense
    return rv
end

function expand_dense!(dense::QMatrix, m::Integer, n::Integer)
    return resize!(dense, m, n)
end

"""check if an intermediate symbol is covered."""
@inline function iscovered(d::Decoder, i::Integer)::Bool
    return length(d.columns[i]) > 0
end

## column/row permutation functions ###
"""swap columns `ci` and `cj` of the constraint matrix."""
@inline function swap_cols!(d::Decoder, ci::Integer, cj::Integer)
    d.colperm[ci], d.colperm[cj] = d.colperm[cj], d.colperm[ci]
    d.colperminv[d.colperm[ci]] = ci
    d.colperminv[d.colperm[cj]] = cj
end

"""swap columns `ui` and `uj` of the dense submatrix, which is denoted by `u` in the spec."""
@inline function swap_dense_cols!(d::Decoder, ui::Integer, uj::Integer)
    d.uperm[ui], d.uperm[uj] = d.uperm[uj], d.uperm[ui]
    d.uperminv[d.uperm[ui]] = ui
    d.uperminv[d.uperm[uj]] = uj
end

"""swap rows `ri` and `rj` of the constraint matrix."""
@inline function swap_rows!(d::Decoder, ri::Integer, rj::Integer)
    d.rowperm[ri], d.rowperm[rj] = d.rowperm[rj], d.rowperm[ri]
    d.rowperminv[d.rowperm[ri]] = ri
    d.rowperminv[d.rowperm[rj]] = rj
end

## functions for subtracting one row from another ##

# isnonzero(v::AbstractFloat) = abs(v) > eps(typeof(v))
isnonzero(v::Real) = !iszero(v)

function subtract!(Vs::AbstractVector{<:Number}; coef, rpi_src::Integer, rpi_dst::Integer)
    Vs[rpi_dst] -= coef*Vs[rpi_src]
    return
end

function subtract!(Vs::AbstractVector{<:AbstractArray}; coef, rpi_src::Integer, rpi_dst::Integer)
    Vs[rpi_dst] .-= coef.*Vs[rpi_src]
    return
end

"""

Subtract `coef_dst / coef_src` multiplied by the `rpi_src`-th constraint from the `rpi_dst`-th
constraint. Updates `Vs` and the dense sub-matrix in-place.
"""
function subtract!(d::Decoder, Vs; rpi_src::Integer, rpi_dst::Integer, coef_src::Tv, coef_dst::Tv, compute_dense::Bool=true) where Tv
    if !isnonzero(coef_dst)
        return
    end
    coef::Tv = coef_dst / coef_src
    if !isnonzero(coef)
        return
    end

    # dense submatrix is stored explicitly
    if compute_dense
        for ui in 1:d.num_inactivated
            d.dense[ui, rpi_dst] -= coef * d.dense[ui, rpi_src]
        end
    end

    # subtract the values
    subtract!(Vs; coef, rpi_src, rpi_dst)
    d.num_rowops += 1
    return
end

"""

Set the dense elements of row `rpi` based on which columns have been inactivated.
"""
function setinactive!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, rpi::Integer)
    rows = rowvals(A)
    vals = nonzeros(A)
    for i in nzrange(A, rpi)
        cpi = rows[i]
        ci = d.colperminv[cpi]
        if ci_is_inactivated(d, ci)
            setdense!(d, rpi, cpi, vals[i])
        end
    end
end

"""

Mark the column with permuted index cpi as decoded. Permutes the
columns such that the decoded column is the rightmost column in I.
"""
function mark_decoded!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, cpi::Integer)
    ci = d.colperminv[cpi]
    ci_is_active(d, ci) || throw(ArgumentError("expected $ci to be active"))
    d.num_decoded += 1
    swap_cols!(d, ci, d.num_decoded)
    update_schedule!(d, A, cpi)
    return
end

"""

Mark the column with permuted index cpi as inactivated. Permutes the
columns such that the inactivated column is the leftmost column of
u. Expands the dense submatrix and permutation vectors as needed.
"""
function mark_inactive!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, cpi::Integer)
    ci = d.colperminv[cpi]
    ci_is_active(d, ci) || throw(ArgumentError("expected $ci to be active"))
    rightmost_active_col = d.num_symbols - d.num_inactivated
    swap_cols!(d, ci, rightmost_active_col)
    d.num_inactivated += 1

    # extend the dense matrix permutation vectors when necessary
    @assert length(d.uperm) == length(d.uperminv)
    @assert d.num_inactivated <= length(d.uperminv) + 1
    if d.num_inactivated > length(d.uperm)
        push!(d.uperm, d.num_inactivated)
        push!(d.uperminv, d.num_inactivated)
    end

    # Update the priority of adjacent rows
    update_schedule!(d, A, cpi)
    return
end

"""

Merge the components containing the constraints indexed by `cpi` and `cpj`, and set the priority
of the resulting component equal to the size of the resulting combined component.
"""
function merge_components!(componentpq, components, cpi::Integer, cpj::Integer)
    x = find_root!(components, cpi)
    y = find_root!(components, cpj)
    if x == y
        return
    end
    root = union!(components, x, y)
    hasx = haskey(componentpq, x)
    hasy = haskey(componentpq, y)
    sizex = hasx ? componentpq[x] : 1
    sizey = hasy ? componentpq[y] : 1
    size = sizex + sizey
    if root == x && hasy
        delete!(componentpq, y)
    elseif root == y && hasx
        delete!(componentpq, x)
    end
    componentpq[root] = size
end

merge_components!(d::Decoder, cpi::Integer, cpj::Integer) = merge_components!(d.componentpq, d.components, cpi, cpj)

"""

Return the indices of non-zero coefficients in the `rpi`-th constraint corresponding to active
symbols. Assumes that there exactly two non-zero coefficients.
"""
function component_indices(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, rpi::Integer)
    Is = nzrange(A, rpi)
    rows = rowvals(A)
    f = (x) -> cpi_is_active(d, x)::Bool
    i::Int = findnext(f, rows, first(Is))::Int
    j::Int = findnext(f, rows, i+1)::Int
    @assert i < j <= last(Is)
    rows[i], rows[j]
end

"""

Update the row schedule after decoding/inactivating column cpi.
"""
function update_schedule!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, cpi::Integer)

    # Decrease the vdegree of neighboring rows by 1
    for rpi in d.columns[cpi]
        if !haskey(d.rowpq, rpi)
            continue
        end

        # Reduce the vdegree by 1, since a neighboring column was decoded or inactivated
        vdegree, degree = d.rowpq[rpi]
        vdegree -= 1
        d.rowpq[rpi] = (vdegree, degree)

        # Update the connected components if the vdegree becomes 2
        if vdegree == 2
            cpi, cpj = component_indices(d, A, rpi)
            merge_components!(d, cpi, cpj)
        end
    end
end

"""Return `true` if the `ci`-th symbol has been decoded."""
ci_is_decoded(d::Decoder, ci::Integer) = ci <= d.num_decoded
cpi_is_decoded(d::Decoder, cpi::Integer) = ci_is_decoded(d, d.colperminv[cpi])

"""Return `true` if the `ci`-th symbol has been inactivated."""
ci_is_inactivated(d::Decoder, ci::Integer) = ci > d.num_symbols - d.num_inactivated
cpi_is_inactivated(d::Decoder, cpi::Integer) = cpi_is_inactivateded(d, d.colperminv[cpi])

"""Return `true` if the `ci`-th symbol is in the active portion of the matrix (denoted by V in the spec.)."""
ci_is_active(d::Decoder, ci::Integer) = !ci_is_decoded(d, ci) && !ci_is_inactivated(d, ci)
cpi_is_active(d::Decoder, cpi::Integer) = ci_is_active(d, d.colperminv[cpi])

"""

Select the next row to process in the diagonalization phase.
"""
function select_row(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC)

    # Drop rows without elements in V, i.e., that don't have non-zero elements in columns
    # corresponding to symbols that aren't decoded or inactivated
    while length(d.rowpq) > 0 && iszero(peek(d.rowpq)[2][1])
        dequeue!(d.rowpq)
    end

    # If the minimum vdegree is 2, inactivate a column part of the largest component
    while length(d.rowpq) > 0 && peek(d.rowpq)[2][1] == 2

        # Drop components corresponding to already decoded/inactivated columns
        while length(d.componentpq) > 0 && !cpi_is_active(d, peek(d.componentpq)[1])
            dequeue!(d.componentpq)
        end

        # Get a column part of the largest component and inactivate it
        cpi = dequeue!(d.componentpq)
        mark_inactive!(d, A, cpi)
    end

    # The matrix is rank deficit if there are no remaining rows
    if length(d.rowpq) == 0
        throw(RankDeficiencyException())
    end

    # Return the row with lowest original degree out of the rows with minimal vdegree.
    rpi = dequeue!(d.rowpq)
    return rpi
end

"""

Zero out any elements of the `rpi`-th constraint below the diagonal.
"""
function zerodiag!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, Vs, rpi::Integer)
    rows = rowvals(A)
    vals = nonzeros(A)
    for i in nzrange(A, rpi)
        cpi = rows[i]
        ci = d.colperminv[cpi]
        coef_dst = vals[i]
        if ci_is_decoded(d, ci) && !ci_is_inactivated(d, ci)
            rpi_src = d.rowperm[ci]
            coef_src = A[cpi, rpi_src]
            subtract!(d, Vs; rpi_src, rpi_dst=rpi, coef_src, coef_dst)
        end
    end
    return
end

"""

Peel away previously decoded symbols from a row. Has to be carried out each time
a row is selected.
"""
function peel_row!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, Vs, rpi::Integer)
    expand_dense!(d)
    setinactive!(d, A, rpi)
    zerodiag!(d, A, Vs, rpi)
    return
end

"""

Perform row and column operations to put the submatrix consisting of the first L-u columns into
diagonal form. Referred to as the first phase in rfc6330.
"""
function diagonalize!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, Vs)
    f = (x) -> cpi_is_active(d, x)
    rows = rowvals(A)
    vals = nonzeros(A)
    while d.num_decoded + d.num_inactivated < d.num_symbols

        # select a constraint to operate on and move it into position
        rpi = select_row(d, A)
        ri = d.rowperminv[rpi]
        swap_rows!(d, ri, d.num_decoded+1)

        # previously decoded symbols are are zeroed out lazily
        peel_row!(d, A, Vs, rpi)

        # swap any non-zero active symbol (not decoded or inactivated) into the first column of V
        Is = nzrange(A, rpi)
        i::Int = findnext(f, rows, first(Is))::Int
        @assert i <= last(Is)
        cpi = rows[i]
        ci = d.colperminv[cpi]
        mark_decoded!(d, A, cpi)

        # inactivate the remaining neighboring symbols
        for j in i+1:last(Is)
            cpi = rows[j]
            ci = d.colperminv[cpi]
            if ci_is_active(d, ci)
                mark_inactive!(d, A, cpi)
                expand_dense!(d)
                coef = vals[j]
                setdense!(d, rpi, cpi, coef)
            end
        end
    end
    return d
end

"""

Zero out all elements of the dense matrix below the diagonal of the `rpj`-th constraint.
"""
function peel_dense_left!(d::Decoder, Vs, rpj::Integer)
    for cj in (d.num_symbols-d.num_inactivated+1):d.num_decoded
        cpj = d.colperm[cj]
        coef = getdense(d, rpj, cpj)
        rpk = d.rowperm[cj]
        coef_src = getdense(d, rpk, cpj)
        subtract!(d, Vs; rpi_src=rpk, rpi_dst=rpj, coef_src, coef_dst=coef)
    end
    return
end

"""

Zero out all elements of u_lower above the diagonal.
"""
function peel_dense_above!(d::Decoder, Vs)
    for rj in (d.num_symbols-d.num_inactivated+1):d.num_decoded
        rpj = d.rowperm[rj]
        cpj = d.colperm[d.num_decoded+1]
        coef_dst = getdense(d, rpj, cpj)
        rpk = d.rowperm[d.num_decoded+1]
        coef_src = getdense(d, rpk, cpj)
        subtract!(d, Vs; rpi_src=rpk, rpi_dst=rpj, coef_src, coef_dst)
    end
    return
end

"""

Solve for the inactivated symbols from the system of equations consisting of the inactivated
columns using Gaussian Elimination.
"""
function solve_dense!(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, Vs)
    nconstraints = size(A, 2)

    # prioritize the remaining constraints by number of non-zero entries
    for ri in (d.num_decoded+1):nconstraints
        rpi = d.rowperm[ri]
        if !haskey(d.rowpq, rpi)
            enqueue!(d.rowpq, rpi, (0, length(nzrange(A, rpi))))
        end
    end

    while length(d.rowpq) > 0 && d.num_decoded < d.num_symbols

        # select a constraint with at least 1 non-zero entry in the section of the matrix
        # consisting of inactivated symbols
        upi = 0
        rpi = 0
        while length(d.rowpq) > 0 && iszero(rpi)
            rpj = dequeue!(d.rowpq)

            # zero out elements below the diagonal
            peel_row!(d, A, Vs, rpj)
            peel_dense_left!(d, Vs, rpj)

            # check if there are any non-zero elements remaining
            for upj in 1:d.num_inactivated
                v = d.dense[upj, rpj]
                uj = d.uperminv[upj]
                cj = _ui2ci(d, uj)
                isdecoded = cj <= d.num_decoded # necessary when operating over floating-point coefficients due to inexact arithmetic
                if !isdecoded && isnonzero(v)
                    rpi = rpj
                    upi = upj
                    break
                end
            end
        end
        if iszero(rpi)
            break
        end
        ri = d.rowperminv[rpi]
        swap_rows!(d, d.num_decoded+1, ri)

        # swap columns such that the non-zero entry found is on the diagonal
        ci = d.num_decoded + 1
        ui = d.uperminv[upi]
        cj = _ui2ci(d, ui)
        uj = _ci2ui(d, ci)
        swap_dense_cols!(d, ui, uj)
        swap_cols!(d, ci, cj)

        # zero out elements above the entry just swapped onto the diagonal
        peel_dense_above!(d, Vs)
        d.num_decoded += 1
    end
    if d.num_decoded < d.num_symbols
        throw(RankDeficiencyException())
    end
    return d
end

"""

Subtract the symbols decoded in solve_dense from the above rows of the constraint matrix.
"""
function backsolve!(d::Decoder, Vs)
    for ri in 1:(d.num_symbols-d.num_inactivated)
        rpi = d.rowperm[ri]
        for upi in 1:d.num_inactivated
            coef = d.dense[upi, rpi]
            @assert !isnan(coef)
            if isnonzero(coef)
                ui = d.uperminv[upi]
                ci = _ui2ci(d, ui)
                cpi = d.colperm[ci]
                rpj = d.rowperm[ci]

                rpi_dst = rpi
                rpi_src = rpj
                coef_dst = coef
                coef_src = getdense(d, rpi_src, cpi)
                subtract!(d, Vs; rpi_src, rpi_dst, coef_src, coef_dst, compute_dense=false)
            end
        end
    end
    return d
end

"""

Return the vector of decoded symbols.
"""
function get_source(d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, Vs::AbstractVector{Tv}) where {Tv}
    get_source!(Vector{Tv}(undef, d.num_symbols), d, A, Vs)
end

"""

In-place version of get_source()
"""
function get_source!(dec::AbstractVector, d::Decoder, A::SparseArrays.AbstractSparseMatrixCSC, Vs)
    length(dec) == d.num_symbols || throw(DimensionMismatch("dec has dimension $(length(dec)), but num_symbols is $(d.num_symbols))"))
    for ri in 1:(d.num_symbols-d.num_inactivated)
        rpi = d.rowperm[ri]
        cpi = d.colperm[ri]
        coef = A[cpi, rpi]
        dec[cpi] = isone(coef) ? Vs[rpi] : Vs[rpi] / coef
    end
    for ri in (d.num_symbols-d.num_inactivated+1):d.num_symbols
        rpi = d.rowperm[ri]
        cpi = d.colperm[ri]
        coef = getdense(d, rpi, cpi)
        if isnonzero(coef)
            dec[cpi] = isone(coef) ? Vs[rpi] : Vs[rpi] / coef
        end
    end
    dec
end

"""

Return `true` if each source symbols is part of at least one constraint, and `false` otherwise.
"""
function all_symbols_covered(d::Decoder)
    for v in d.columns
        if length(v) == 0
            return false
        end
    end
    return true
end

"""
    decode(A::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}, Vs::AbstractVector; decoder::Decoder{Tv,Ti}=Decoder(A)) where {Tv,Ti}

Solve the system of equations `transpose(A)*x = Vs` for `x` using the inactivation decoding
algorithm. All algorithm state is stored in the `Decoder` struct, i.e., neither `A` or `Vs` is
mutated.
"""
function decode(A::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}, Vs::AbstractVector; decoder::Decoder{Tv,Ti}=Decoder(A)) where {Tv,Ti}
    k, n = size(A)
    length(Vs) == n || throw(DimensionMismatch("A has dimensions $(size(A)), but Vs has dimension $(length(Vs))"))
    n >= k || throw(RankDeficiencyException())
    all_symbols_covered(decoder) || throw(RankDeficiencyException())
    diagonalize!(decoder, A, Vs)
    solve_dense!(decoder, A, Vs)
    backsolve!(decoder, Vs)
    get_source(decoder, A, Vs)
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
# ---------------------------------------------------------------------------- #
#
#   basisFuncTriangleLeg.jl
#
#   Get Legendre basis functions for Triangle
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

function basisFuncTriangleLeg( P::P1, s::Vector{Float64}, t::Vector{Float64} )

  dim = 2
  nx  = length( s )

  # Linear

  nphi = 3
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  # value
  phi[:,1] = 1
  phi[:,2] = sqrt(6).*(s - t)
  phi[:,3] = 3.*sqrt(2).*(s + t - 2./3.)

  # derivative in x
  dphi[:,1,1] =  0
  dphi[:,2,1] =  sqrt(6)
  dphi[:,3,1] =  3.*sqrt(2)

  # derivative in y
  dphi[:,1,1] =  0
  dphi[:,2,1] = -sqrt(6)
  dphi[:,3,1] =  3.*sqrt(2)

  # transpose
  phi   = phi'

  dphi2 = fill( 0.0, size(dphi,2), size(dphi,1), 2 )
  dphi2[:,:,1] = dphi[:,:,1]'
  dphi2[:,:,2] = dphi[:,:,2]'

  return phi, dphi2

end

function basisFuncTriangleLeg( P::P2, s::Vector{Float64}, t::Vector{Float64} )

  dim = 2
  nx  = length( s )

  # Quadratric

  nphi = 6
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  # value
  phi[:,1] = 1
  phi[:,2] = sqrt(6).*(s - t)
  phi[:,3] = 3.*sqrt(2).*(s + t - 2./3.)
  phi[:,4] = 3.*(s - t).*(4 - 5.*(s + t))
  phi[:,5] = sqrt(3).*(3 - 12.*(s + t) + 10.*(s + t).*(s + t))
  phi[:,6] = sqrt(15).*((s + t).*(s + t) - 3.*(s - t).*(s - t))/2.

  # derivative in x
  dphi[:,1,1] =  0
  dphi[:,2,1] =  sqrt(6)
  dphi[:,3,1] =  3.*sqrt(2)
  dphi[:,4,1] =  6.*(2 - 5.*s)
  dphi[:,5,1] = -4.*sqrt(3).*(3 - 5.*(s + t))
  dphi[:,6,1] = -2.*sqrt(15).*(s - 2.*t)

  # derivative in y
  dphi[:,1,1] =  0
  dphi[:,2,1] = -sqrt(6)
  dphi[:,3,1] =  3.*sqrt(2)
  dphi[:,4,1] = -6.*(2 - 5.*t)
  dphi[:,5,1] = -4.*sqrt(3).*(3 - 5.*(s + t))
  dphi[:,6,1] = -2.*sqrt(15).*(t - 2.*s)

  # transpose
  phi   = phi'

  dphi2 = fill( 0.0, size(dphi,2), size(dphi,1), 2 )
  dphi2[:,:,1] = dphi[:,:,1]'
  dphi2[:,:,2] = dphi[:,:,2]'

  return phi, dphi2

end

function basisFuncTriangleLeg( P::P3, s::Vector{Float64}, t::Vector{Float64} )

  dim = 2
  nx  = length( s )

  nphi = 10
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  # value
  phi[:, 1] = 1
  phi[:, 2] = sqrt(6).*(s - t)
  phi[:, 3] = 3.*sqrt(2).*(s + t - 2./3.)
  phi[:, 4] = 3.*(s - t).*(4 - 5.*(s + t))
  phi[:, 5] = sqrt(3).*(3 - 12.*(s + t) + 10.*(s + t).*(s + t))
  phi[:, 6] = sqrt(15).*((s + t).*(s + t) - 3.*(s - t).*(s - t))/2.
  phi[:, 7] = 2.*sqrt(3).*(s - t).*(10 - 30.*(s + t) + 21.*(s + t).*(s + t))
  phi[:, 8] =   sqrt(7).*(s - t).*(5.*(s - t).*(s - t) - 3.*(s + t).*(s + t))
  phi[:, 9] = 2.*(5.*((7.*(s + t) - 12).*(s + t) + 6).*(s + t) - 4)
  phi[:,10] = sqrt(5).*((s + t).*(s + t) - 3.*(s - t).*(s - t)).*(6 - 7.*(s + t))

  # derivative in x
  dphi[:, 1,1] =  0
  dphi[:, 2,1] =  sqrt(6)
  dphi[:, 3,1] =  3.*sqrt(2)
  dphi[:, 4,1] =  6.*(2 - 5.*s)
  dphi[:, 5,1] = -4.*sqrt(3).*(3 - 5.*(s + t))
  dphi[:, 6,1] = -2.*sqrt(15).*(s - 2.*t)
  dphi[:, 7,1] =  2.*sqrt(3).*(10.*(1 - 6.*s) + 21.*(s + t).*(3.*s - t))
  dphi[:, 8,1] =  6.*sqrt(7).*(3.*(s - t).*(s - t) - 2.*s.*s)
  dphi[:, 9,1] =  30.*(2 - 8.*(s + t) + 7.*(s + t).*(s + t))
  dphi[:,10,1] =  6.*sqrt(5).*(7.*((s - t).*(s - t) - 2.*t.*t) - 4.*(s - 2.*t))

  dphi[:, 1,1] =  0
  dphi[:, 2,1] = -sqrt(6)
  dphi[:, 3,1] =  3.*sqrt(2)
  dphi[:, 4,1] = -6.*(2 - 5.*t)
  dphi[:, 5,1] = -4.*sqrt(3).*(3 - 5.*(s + t))
  dphi[:, 6,1] = -2.*sqrt(15).*(t - 2.*s)
  dphi[:, 7,1] = -2.*sqrt(3).*(10.*(1 - 6.*t) - 21.*(s + t).*(s - 3.*t))
  dphi[:, 8,1] = -6.*sqrt(7).*(3.*(s - t).*(s - t) - 2.*t.*t)
  dphi[:, 9,1] =  30.*(2 - 8.*(s + t) + 7.*(s + t).*(s + t))
  dphi[:,10,2] =  6.*sqrt(5).*(7.*((s - t).*(s - t) - 2.*s.*s) - 4.*(t - 2.*s))

  # transpose
  phi   = phi'

  dphi2 = fill( 0.0, size(dphi,2), size(dphi,1), 2 )
  dphi2[:,:,1] = dphi[:,:,1]'
  dphi2[:,:,2] = dphi[:,:,2]'

  return phi, dphi2

end

function basisFuncTriangleLeg( P::P4, s::Vector{Float64}, t::Vector{Float64} )

  dim = 2
  nx  = length( s )

  nphi = 15
  phi  = Array{Float64}(nx, nphi)
  dphi = Array{Float64}(nx, nphi, dim)

  u = s + t
  v = s - t

  # value
  phi[:, 1] = 1

  phi[:, 2] = sqrt(6).*v
  phi[:, 3] = 3.*sqrt(2).*(u - 2./3.)

  phi[:, 4] = 3.*v.*(4 - 5.*u)
  phi[:, 5] = sqrt(3).*(3 - 12.*u + 10.*u.*u)
  phi[:, 6] = sqrt(15).*(u.*u - 3.*v.*v)/2.

  phi[:, 7] = 2.*sqrt(3).*v.*(10 - 30.*u + 21.*u.*u)
  phi[:, 8] =   sqrt(7).*v.*(5.*v.*v - 3.*u.*u)
  phi[:, 9] = 2.*(5.*((7.*u - 12).*u + 6).*u - 4)
  phi[:,10] = sqrt(5).*(u.*u - 3.*v.*v).*(6 - 7.*u)

  phi[:,11] = sqrt(15).*v.*(20 - 105.*u + 168.*u.*u - 84.*u.*u.*u)
  phi[:,12] = 0.5.*sqrt(35).*v.*(5.*v.*v - 3.*u.*u).*(8 - 9.*u)
  phi[:,13] = sqrt(5).*(5 - 60.*u + 210.*u.*u - 280.*u.*u.*u + 126.*u.*u.*u.*u)
  phi[:,14] = 2.5.*(u.*u - 3.*v.*v).*(21 - 4.*u.*(14 - 9.*u))
  phi[:,15] = 0.375.*sqrt(5).*(3.*u.*u.*u.*u - 30.*u.*u.*v.*v + 35.*v.*v.*v.*v)

  # derivative in x
  dphi[:, 1,1] =  0

  dphi[:, 2,1] =  sqrt(6)
  dphi[:, 3,1] =  3.*sqrt(2)

  dphi[:, 4,1] =  3.*(4 - 5.*(u + v))
  dphi[:, 5,1] = -4.*sqrt(3).*(3 - 5.*u)
  dphi[:, 6,1] =  sqrt(15).*(u - 3.*v)

  dphi[:, 7,1] =  2.*sqrt(3).*(10 - 30.*v - 3.*u.*(10 - 7.*(u + 2.*v)))
  dphi[:, 8,1] = -3.*sqrt(7).*(u.*u + 2.*u.*v - 5.*v.*v)
  dphi[:, 9,1] =  30.*(2 - 8.*u + 7.*u.*u)
  dphi[:,10,1] = -3.*sqrt(5).*(7.*u.*u + v.*(12 - 7.*v) - 2.*u.*(2 + 7.*v))

  dphi[:,11,1] =  sqrt(15).*(5.*(4 - 21.*v) - 21.*u.*(5 - 16.*v - 4.*u.*(2 - u - 3.*v)))
  dphi[:,12,1] = -1.5.*sqrt(35).*(u.*u.*(8 - 9.*u - 27.*v) + u.*v.*(16 + 45.*v) - 5.*v.*v.*(8 - 3.*v))
  dphi[:,13,1] = -12.*sqrt(5).*(5 - 7.*u.*(5 - 2.*u.*(5 - 3.*u)))
  dphi[:,14,1] = -15.*(4.*u.*u.*(7 - 6.*u + 9.*v) - u.*(7 + 4.*v.*(14 - 9.*v)) + 7.*v.*(3 - 4.*v))
  dphi[:,15,1] =  1.5.*sqrt(5).*(3.*u.*u.*(u - 5.*v) - 5.*v.*v.*(3.*u - 7.*v))

  # derivative in y
  dphi[:, 1,2] =  0

  dphi[:, 2,2] = -sqrt(6)
  dphi[:, 3,2] =  3.*sqrt(2)

  dphi[:, 4,2] = -3.*(4 - 5.*(u - v))
  dphi[:, 5,2] = -4.*sqrt(3).*(3 - 5.*u)
  dphi[:, 6,2] =  sqrt(15).*(u + 3.*v)

  dphi[:, 7,2] = -2.*sqrt(3).*(10 + 30.*v - 3.*u.*(10 - 7.*(u - 2.*v)))
  dphi[:, 8,2] =  3.*sqrt(7).*(u.*u - 2.*u.*v - 5.*v.*v)
  dphi[:, 9,2] =  30.*(2 - 8.*u + 7.*u.*u)
  dphi[:,10,2] = -3.*sqrt(5).*(7.*u.*u - v.*(12 + 7.*v) - 2.*u.*(2 - 7.*v))

  dphi[:,11,2] = -sqrt(15).*(5.*(4 + 21.*v) - 21.*u.*(5 + 16.*v - 4.*u.*(2 - u + 3.*v)))
  dphi[:,12,2] =  1.5.*sqrt(35).*(u.*u.*(8 - 9.*u + 27.*v) - u.*v.*(16 - 45.*v) - 5.*v.*v.*(8 + 3.*v))
  dphi[:,13,2] = -12.*sqrt(5).*(5 - 7.*u.*(5 - 2.*u.*(5 - 3.*u)))
  dphi[:,14,2] = -15.*(4.*u.*u.*(7 - 6.*u - 9.*v) - u.*(7 - 4.*v.*(14 + 9.*v)) - 7.*v.*(3 + 4.*v))
  dphi[:,15,2] =  1.5.*sqrt(5).*(3.*u.*u.*(u + 5.*v) - 5.*v.*v.*(3.*u + 7.*v))

  # transpose
  phi   = phi'

  dphi2 = fill( 0.0, size(dphi,2), size(dphi,1), 2 )
  dphi2[:,:,1] = dphi[:,:,1]'
  dphi2[:,:,2] = dphi[:,:,2]'

  return phi, dphi2

end


text:
"""
    importicd10amcodes(fname)

Import the ICD-10-AM electronic code lists into the package.

Electronic code lists can be obtained from the [Independent Hospital Pricing
Authority](https://www.ihpa.gov.au/what-we-do/icd-10-am-achi-acs-classification).

This function just imports the `disease.txt` table. Provide the path to the `disease.txt` file.
"""
function importicd10amcodes(fname)
  icdcodes = CSV.File(fname)
  Level = Int8.(icdcodes.Level)
  code_id = ICD10AM.(icdcodes.code_id)
  dagger = Bool.(icdcodes.dagger)
  asterisk = Bool.(icdcodes.asterisk)
  valid = Bool.(icdcodes.valid)
  aust_code = Bool.(icdcodes.aust_code)
  ascii_desc = String.(icdcodes.ascii_desc)
  ascii_short_desc = String.(icdcodes.ascii_short_desc)
  dtfmt = DateFormat("d/m/y")
  effective_from = passmissing(Date).(icdcodes.effective_from, dtfmt)
  inactive = passmissing(Date).(icdcodes.inactive, dtfmt)
  reactivated = passmissing(Date).(icdcodes.reactivated, dtfmt)
  sex = recode(
    categorical(icdcodes.sex; levels = 1:2, compress = true),
    1 => "male",
    2 => "female",
  )
  Stype = recode(
    categorical(icdcodes.Stype; levels = 1:2, compress = true),
    1 => "fatal",
    2 => "warning",
  )
  AgeL = passmissing(ICD10AMAge).(icdcodes.AgeL)
  AgeH = passmissing(ICD10AMAge).(icdcodes.AgeH)
  Atype = recode(
    categorical(icdcodes.Atype; levels = 1:2, compress = true),
    1 => "fatal",
    2 => "warning",
  )
  RDiag = Bool.(icdcodes.RDiag)
  Morph_Code = Bool.(icdcodes.Morph_Code)
  concept_change = passmissing(Date).(icdcodes.concept_change, dtfmt)
  UnacceptPDx = Bool.(icdcodes.UnacceptPDx)
  global ICD10AMcodes = Table(;
    level = Level,
    icdcode = code_id,
    dagger = dagger,
    asterisk = asterisk,
    validforcoding = valid,
    australiancode = aust_code,
    description = ascii_desc,
    shortdescription = ascii_short_desc,
    effectivefrom = effective_from,
    inactive = inactive,
    reactivated = reactivated,
    sex = sex,
    stype = Stype,
    agel = AgeL,
    ageh = AgeH,
    atype = Atype,
    rarenotifiabledx = RDiag,
    morphcoderequired = Morph_Code,
    conceptchange = concept_change,
    unacceptpdx = UnacceptPDx,
  )
  save(
    normpath(@__DIR__, "..", "data", "icd10amcodes.jld2"),
    Dict("icd10amcodes" => Table(ICD10AMcodes)),
  )
end

"""
    importachicodes(fname)

Import the ACHI electronic code lists into the package.

Electronic code lists can be obtained from the [Independent Hospital Pricing
Authority](https://www.ihpa.gov.au/what-we-do/icd-10-am-achi-acs-classification).

This function just imports the `interven.txt` table. Provide the path to the `interven.txt` file.
"""
function importachicodes(fname)
  achicodes = CSV.File(fname)
  Code_id = ACHI.(achicodes.Code_id)
  Code_id = ACHI.(achicodes.Code_id)
  Block = Int16.(achicodes.Block)
  dtfmt = DateFormat("d/m/y")
  ascii_desc = String.(achicodes.ascii_desc)
  ascii_short_desc = String.(achicodes.ascii_short_desc)
  effective_from = passmissing(Date).(achicodes.effective_from, dtfmt)
  inactive = passmissing(Date).(achicodes.inactive, dtfmt)
  sex = recode(
    categorical(achicodes.sex; levels = 1:2, compress = true),
    1 => "male",
    2 => "female",
  )
  Stype = recode(
    categorical(achicodes.Stype; levels = 1:2, compress = true),
    1 => "fatal",
    2 => "warning",
  )
  AgeL = passmissing(ICD10AMAge).(achicodes.AgeL)
  AgeL = passmissing(ICD10AMAge).(achicodes.AgeL)
  AgeH = passmissing(ICD10AMAge).(achicodes.AgeH)
  Atype = recode(
    categorical(achicodes.Atype; levels = 1:2, compress = true),
    1 => "fatal",
    2 => "warning",
  )
  global ACHIcodes = Table(;
    achicode = Code_id,
    block = Block,
    description = ascii_desc,
    shortdescription = ascii_short_desc,
    effectivefrom = effective_from,
    inactive = inactive,
    sex = sex,
    stype = Stype,
    agel = AgeL,
    ageh = AgeH,
    atype = Atype,
  )
  save(
    normpath(@__DIR__, "..", "data", "achicodes.jld2"),
    Dict("achicodes" => Table(ACHIcodes)),
  )
end


text:
# Autogenerated wrapper script for ParseGen_jll for armv6l-linux-musleabihf-cxx03
export libparsegen, parsegen_calc

JLLWrappers.@generate_wrapper_header("ParseGen")
JLLWrappers.@declare_library_product(libparsegen, "libparsegen.so")
JLLWrappers.@declare_executable_product(parsegen_calc)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libparsegen,
        "lib/libparsegen.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        parsegen_calc,
        "bin/parsegen-calc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()