Add developer docs and a note about abandonment (#365)

Author: timholy

README.md

@@ -10,6 +10,12 @@ Documentation:
 
 **NEWS** v0.9 was a breaking release. See the [news](NEWS.md) for details on how to update.
 
+**NEWS** This package is officially looking for a maintainer.
+The original authors are only fixing bugs that affect them personally.
+However, to facilitate transition to a long-term maintainer,
+for a period of time they will make a concerted attempt to review pull requests.
+Step forward soon to avoid the risk of not being able to take advantage of such offers of support.
+
 This package implements a variety of interpolation schemes for the
 Julia language.  It has the goals of ease-of-use, broad algorithmic
 support, and exceptional performance.

docs/Project.toml

@@ -2,4 +2,4 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 
 [compat]
-Documenter = "~0.21"
+Documenter = "0.24"

docs/make.jl

@@ -8,7 +8,9 @@ pages=["Home" => "index.md",
         "Interpolation algorithms" => "control.md",
         "Extrapolation" => "extrapolation.md",
         "Convenience Constructors" => "convenience-construction.md",
+        "Developer documentation" => "devdocs.md",
         "Library" => "api.md"]
 )
 
-deploydocs(repo="github.com/JuliaMath/Interpolations.jl")
+deploydocs(repo="github.com/JuliaMath/Interpolations.jl",
+            push_preview=true)

docs/src/control.md

@@ -63,41 +63,44 @@ sitp(5.6, 7.1) # approximately log(5.6 + 7.1)
 These use a very similar syntax to BSplines, with the major exception
 being that one does not get to choose the grid representation (they
 are all `OnGrid`). As such one must specify a set of coordinate arrays
-defining the knots of the array.
+defining the nodes of the array.
 
 In 1D
 ```julia
 A = rand(20)
-A_x = collect(1.0:2.0:40.0)
-knots = (A_x,)
-itp = interpolate(knots, A, Gridded(Linear()))
+A_x = 1.0:2.0:40.0
+nodes = (A_x,)
+itp = interpolate(nodes, A, Gridded(Linear()))
 itp(2.0)
 ```
 
-The spacing between adjacent samples need not be constant, you can use the syntax
+The spacing between adjacent samples need not be constant; indeed, if they
+are constant, you'll get better performance with `scaled`.
+
+The general syntax is
 ```julia
-itp = interpolate(knots, A, options...)
+itp = interpolate(nodes, A, options...)
 ```
-where `knots = (xknots, yknots, ...)` to specify the positions along
+where `nodes = (xnodes, ynodes, ...)` specifies the positions along
 each axis at which the array `A` is sampled for arbitrary ("rectangular") samplings.
 
 For example:
 ```julia
 A = rand(8,20)
-knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
-itp = interpolate(knots, A, Gridded(Linear()))
+nodes = ([x^2 for x = 1:8], [0.2y for y = 1:20])
+itp = interpolate(nodes, A, Gridded(Linear()))
 itp(4,1.2)  # approximately A[2,6]
 ```
 One may also mix modes, by specifying a mode vector in the form of an explicit tuple:
 ```julia
-itp = interpolate(knots, A, (Gridded(Linear()),Gridded(Constant())))
+itp = interpolate(nodes, A, (Gridded(Linear()),Gridded(Constant())))
 ```
 
 Presently there are only three modes for gridded:
 ```julia
 Gridded(Linear())
 ```
-whereby a linear interpolation is applied between knots,
+whereby a linear interpolation is applied between nodes,
 ```julia
 Gridded(Constant())
 ```

docs/src/devdocs.md

@@ -0,0 +1,170 @@
+# Developer documentation
+
+Conceptually, `Interpolations.jl` supports two operations: construction and usage of interpolants.
+
+## Interpolant construction
+
+Construction creates the interpolant object. In some situations this is relatively trivial:
+for example, when using only `NoInterp`, `Constant`, or `Linear` interpolation schemes,
+construction essentially corresponds to recording the array of values and the "settings" (the interpolation scheme) specified at the time of
+construction. This case is simple because interpolated values may be efficiently computed directly from the on-grid values supplied at construction time: ``(1-\Delta x) a_i  + \Delta x a_{i+1}`` reconstructs ``a_i`` when ``\Delta x = 0``.
+
+For `Quadratic` and higher orders, efficient computation requires that the array of values be *prefiltered*.
+This essentially corresponds to "inverting" the computation that will be performed during interpolation, so as to approximately reconstruct the original values at on-grid points. Generally speaking this corresponds to solving a nearly-tridiagonal system of equations, inverting an underlying interpolation
+scheme such as ``p(\Delta x) \tilde a_{i-1} + q(\Delta x) \tilde a_i + p(1-\Delta x) \tilde a_{i+1}`` for some functions ``p`` and ``q`` (see [`Quadratic`](@ref) for further details).
+Here ``\tilde a`` is the pre-filtered version of `a`, designed so that
+substituting ``\Delta x = 0`` (for which one may *not* get 0 and 1 for the ``p`` and ``q`` calls, respectively) approximately recapitulates ``a_i``.
+
+The exact system of equations to be solved depends on the interpolation order and boundary conditions.
+Boundary conditions often introduce deviations from perfect tridiagonality; these "extras" are handled efficiently by the `WoodburyMatrices` package.
+These computations are implemented independently along each axis using the `AxisAlgorithms` package.
+
+In the `doc` directory there are some old files that give some of the mathematical details.
+A useful reference is:
+
+Thévenaz, Philippe, Thierry Blu, and Michael Unser. "Interpolation revisited." IEEE Transactions on Medical Imaging 19.7 (2000): 739-758.
+
+!!! note
+    As an application of these concepts, note that supporting quadratic or cubic interpolation for `Gridded` would
+    only require that someone implement prefiltering schemes for non-uniform grids;
+    it's just a question of working out a little bit of math.
+
+## Interpolant usage
+
+Usage occurs when evaluating `itp(x, y...)`, or `Interpolations.gradient(itp, x, y...)`, etc.
+Usage itself involves two sub-steps: computation of the *weights* and then performing the *interpolation*.
+
+### Weight computation
+
+Weights depend on the interpolation scheme and the location `x, y...` but *not* the coefficients of the array we are interpolating.
+Consequently there are many circumstances where one might want to reuse previously-computed weights, and `Interpolations.jl` has been carefully designed
+with that kind of reuse in mind.
+
+The key concept here is the [`Interpolations.WeightedIndex`](@ref), and there is
+no point repeating its detailed docstring here.
+It suffices to add that `WeightedIndex` is actually an abstract type, with two
+concrete subtypes:
+
+- `WeightedAdjIndex` is for indexes that will address *adjacent* points of the coefficient array (ones where the index increments by 1 along the corresponding dimension). These are used when prefiltering produces padding that can be used even at the edges, or for schemes like `Linear` interpolation which require no padding.
+- `WeightedArbIndex` stores both the weight and index associated with each accessed grid point, and can therefore encode grid access patterns. These are used in specific circumstances--a prime example being periodic boundary conditions--where the coefficients array may be accessed at something other than adjacent locations.
+
+`WeightedIndex` computation reflects the interpolation scheme (e.g., `Linear` or `Quadratic`) and also whether one is computing values, `gradient`s, or `hessian`s. The handling of derivatives will be described further below.
+
+### Interpolation
+
+General `AbstractArray`s may be indexed with `WeightedIndex` indices,
+and the result produces the interpolated value. In other words, the end result
+is just `itp.coefs[wis...]`, where `wis` is a tuple of `WeightedIndex` indices.
+
+Derivatives along a particular axis can be computed just by substituting a component of `wis` for one that has been designed to compute derivatives rather than values.
+
+As a demonstration, let's see how the following computation occurs:
+
+```jldoctest derivs; setup = :(using Interpolations)
+julia> A = reshape(1:27, 3, 3, 3)
+3×3×3 reshape(::UnitRange{Int64}, 3, 3, 3) with eltype Int64:
+[:, :, 1] =
+ 1  4  7
+ 2  5  8
+ 3  6  9
+
+[:, :, 2] =
+ 10  13  16
+ 11  14  17
+ 12  15  18
+
+[:, :, 3] =
+ 19  22  25
+ 20  23  26
+ 21  24  27
+
+julia> itp = interpolate(A, BSpline(Linear()));
+
+julia> x = (1.2, 1.4, 1.7)
+(1.2, 1.4, 1.7)
+
+julia> itp(x...)
+8.7
+```
+
+!!! note
+    By using the debugging facilities of an IDE like Juno or VSCode,
+    or using Debugger.jl from the REPL, you can easily step in to
+    the call above and follow along with the description below.
+
+Aside from details such as bounds-checking, the key call is to [`Interpolations.weightedindexes`](@ref):
+
+```jldoctest derivs
+julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp)..., x)
+(Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8, 0.19999999999999996)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.6000000000000001, 0.3999999999999999)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.30000000000000004, 0.7)))
+
+julia> wis[1]
+Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8, 0.19999999999999996))
+
+julia> wis[2]
+Interpolations.WeightedAdjIndex{2,Float64}(1, (0.6000000000000001, 0.3999999999999999))
+
+julia> wis[3]
+Interpolations.WeightedAdjIndex{2,Float64}(1, (0.30000000000000004, 0.7))
+
+julia> A[wis...]
+8.7
+```
+
+You can see that each of `wis` corresponds to a specific position: 1.2, 1.4, and 1.7 respectively. We can index `A` at `wis`, and it returns the value of `itp(x...)`, which here is just
+
+      0.8 * A[1, wis[2], wis[3]] + 0.2 * A[2, wis[2], wis[3]]
+    = 0.6 * (0.8 * A[1, 1, wis[3]] + 0.2 * A[2, 1, wis[3]]) +
+      0.4 * (0.8 * A[1, 2, wis[3]] + 0.2 * A[2, 2, wis[3]])
+    = 0.3 * (0.6 * (0.8 * A[1, 1, 1] + 0.2 * A[2, 1, 1]) +
+             0.4 * (0.8 * A[1, 2, 1] + 0.2 * A[2, 2, 1])  ) +
+      0.7 * (0.6 * (0.8 * A[1, 1, 2] + 0.2 * A[2, 1, 2]) +
+             0.4 * (0.8 * A[1, 2, 2] + 0.2 * A[2, 2, 2])  )
+
+This computed the value of `itp` at `x...` because we called `weightedindexes` with just a single function, [`Interpolations.value_weights`](@ref) (meaning, "the weights needed to compute the value").
+
+!!! note
+    Remember that prefiltering is not used for `Linear` interpolation.
+    In a case where prefiltering is used, we would substitute `itp.coefs[wis...]` for `A[wis...]` above.
+
+To compute derivatives, we *also* pass additional functions like [`Interpolations.gradient_weights`](@ref):
+
+```
+julia> wis = Interpolations.weightedindexes((Interpolations.value_weights, Interpolations.gradient_weights), Interpolations.itpinfo(itp)..., x)
+((Interpolations.WeightedAdjIndex{2,Float64}(1, (-1.0, 1.0)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.6000000000000001, 0.3999999999999999)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.30000000000000004, 0.7))), (Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8, 0.19999999999999996)), Interpolations.WeightedAdjIndex{2,Float64}(1, (-1.0, 1.0)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.30000000000000004, 0.7))), (Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8, 0.19999999999999996)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.6000000000000001, 0.3999999999999999)), Interpolations.WeightedAdjIndex{2,Float64}(1, (-1.0, 1.0))))
+
+julia> wis[1]
+(Interpolations.WeightedAdjIndex{2,Float64}(1, (-1.0, 1.0)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.6000000000000001, 0.3999999999999999)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.30000000000000004, 0.7)))
+
+julia> wis[2]
+(Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8, 0.19999999999999996)), Interpolations.WeightedAdjIndex{2,Float64}(1, (-1.0, 1.0)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.30000000000000004, 0.7)))
+
+julia> wis[3]
+(Interpolations.WeightedAdjIndex{2,Float64}(1, (0.8, 0.19999999999999996)), Interpolations.WeightedAdjIndex{2,Float64}(1, (0.6000000000000001, 0.3999999999999999)), Interpolations.WeightedAdjIndex{2,Float64}(1, (-1.0, 1.0)))
+
+julia> A[wis[1]...]
+1.0
+
+julia> A[wis[2]...]
+3.000000000000001
+
+julia> A[wis[3]...]
+9.0
+```
+In this case you can see that `wis` is a 3-tuple-of-3-tuples. `A[wis[i]...]` can be used to compute the `i`th component of the gradient.
+
+If you look carefully at each of the entries in `wis`, you'll see that
+each "inner" 3-tuple copies two of the three elements in the `wis` we
+obtained when we called `weightedindexes` with just `value_weights` above.
+`wis[1]` replaces the first entry with
+a weighted index having weights `(-1.0, 1.0)`, which corresponds to computing the *slope* along this dimension.
+Likewise `wis[2]` and `wis[3]` replace the second and third value-index, respectively, with the same slope computation.
+
+Hessian computation is quite similar, with the difference that one sometimes needs to replace two different indices or the same index with a set of weights corresponding to a second derivative.
+
+Consequently derivatives along particular directions are computed simply by "weight replacement" along the corresponding dimensions.
+
+The code to do this replacement is a bit complicated due to the need to support arbitrary dimensionality in a manner that allows Julia's type-inference to succeed.
+It makes good use of *tuple* manipulations, sometimes called "lispy tuple programming."
+You can search Julia's discourse forum for more tips about how to program this way.
+It could alternatively be done using generated functions, but this would increase compile time considerably and can lead to world-age problems.

src/Interpolations.jl

@@ -40,9 +40,12 @@ import Base: convert, size, axes, promote_rule, ndims, eltype, checkbounds, axes
 
 abstract type Flag end
 abstract type InterpolationType <: Flag end
+"`NoInterp()` indicates that the corresponding axis must use integer indexing (no interpolation is to be performed)"
 struct NoInterp <: InterpolationType end
 abstract type GridType <: Flag end
+"`OnGrid()` indicates that the boundary condition applies at the first & last nodes"
 struct OnGrid <: GridType end
+"`OnCell()` indicates that the boundary condition applies a half-gridspacing beyond the first & last nodes"
 struct OnCell <: GridType end
 
 const DimSpec{T} = Union{T,Tuple{Vararg{Union{T,NoInterp}}},NoInterp}
@@ -73,23 +76,27 @@ An abstract type with one of the following values (see the help for each for det
 - `InPlace(gt)`
 - `InPlaceQ(gt)`
 
-where `gt` is the grid type, e.g., `OnGrid()` or `OnCell()`. `OnGrid` means that the boundary
-condition "activates" at the first and/or last integer location within the interpolation region,
-`OnCell` means the interpolation extends a half-integer beyond the edge before
-activating the boundary condition.
+where `gt` is the grid type, e.g., [`OnGrid()`](@ref) or [`OnCell()`](@ref).
 """
 abstract type BoundaryCondition <: Flag end
 # Put the gridtype into the boundary condition, since that's all it affects (see issue #228)
 # Nothing is used for extrapolation
+"`Throw(gt)` causes beyond-the-edge extrapolation to throw an error"
 struct Throw{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+"`Flat(gt)` sets the extrapolation slope to zero"
 struct Flat{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+"`Line(gt)` uses a constant slope for extrapolation"
 struct Line{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
 struct Free{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+"`Periodic(gt)` applies periodic boundary conditions"
 struct Periodic{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+"`Reflect(gt)` applies reflective boundary conditions"
 struct Reflect{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+"`InPlace(gt)` is a boundary condition that allows prefiltering to occur in-place (it typically requires padding)"
 struct InPlace{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
-# InPlaceQ is exact for an underlying quadratic. This is nice for ground-truth testing
-# of in-place (unpadded) interpolation.
+"""`InPlaceQ(gt)` is similar to `InPlace(gt)`, but is exact when the values being interpolated
+arise from an underlying quadratic. It is primarily useful for testing purposes,
+allowing near-exact (to machine precision) comparisons against ground truth."""
 struct InPlaceQ{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
 const Natural = Line
 
@@ -114,7 +121,7 @@ twotuple(x, y) = (x, y)
 
 Return the `bounds` of the domain of `itp` as a tuple of `(min, max)` pairs for each coordinate. This is best explained by example:
 
-```jldoctest
+```jldoctest; setup = :(using Interpolations)
 julia> itp = interpolate([1 2 3; 4 5 6], BSpline(Linear()));
 
 julia> bounds(itp)
@@ -337,7 +344,7 @@ Compute the weights for interpolation of the value at an offset `δx` from the "
 
 # Example
 
-```jldoctest
+```jldoctest; setup = :(using Interpolations)
 julia> Interpolations.value_weights(Linear(), 0.2)
 (0.8, 0.2)
 ```
@@ -354,7 +361,7 @@ Compute the weights for interpolation of the gradient at an offset `δx` from th
 
 # Example
 
-```jldoctest
+```jldoctest; setup = :(using Interpolations)
 julia> Interpolations.gradient_weights(Linear(), 0.2)
 (-1.0, 1.0)
 ```
@@ -371,7 +378,7 @@ Compute the weights for interpolation of the hessian at an offset `δx` from the
 
 # Example
 
-```jldoctest
+```jldoctest; setup = :(using Interpolations)
 julia> Interpolations.hessian_weights(Linear(), 0.2)
 (0.0, 0.0)
 ```