Improve docs

src/b-splines/constant.jl

@@ -1,20 +1,41 @@
 immutable Constant <: Degree{0} end
 
+"""
+Constant b-splines are *nearest-neighbor* interpolations, and effectively
+return `A[round(Int,x)]` when interpolating
+"""
+Constant
+
+"""
+`define_indices_d` for a constant b-spline calculates `ix_d = round(Int,x_d)`
+"""
 function define_indices_d(::Type{BSpline{Constant}}, d, pad)
     symix, symx = symbol("ix_",d), symbol("x_",d)
     :($symix = clamp(round(Int, $symx), 1, size(itp, $d)))
 end
 
+"""
+`coefficients` for a constant b-spline simply sets `c_d = 1` for compatibility
+with the general b-spline framework
+"""
 function coefficients(::Type{BSpline{Constant}}, N, d)
     sym, symx = symbol("c_",d), symbol("x_",d)
     :($sym = 1)
 end
 
+"""
+`gradient_coefficients` for a constant b-spline simply sets `c_d = 0` for
+compatibility with the general b-spline framework
+"""
 function gradient_coefficients(::Type{BSpline{Constant}}, d)
     sym, symx = symbol("c_",d), symbol("x_",d)
     :($sym = 0)
 end
 
+"""
+`hessian_coefficients` for a constant b-spline simply sets `c_d = 0` for
+compatibility with the general b-spline framework
+"""
 function hessian_coefficients(::Type{BSpline{Constant}}, d)
     sym = symbol("c_",d)
     :($sym = 0)

src/b-splines/cubic.jl

@@ -5,14 +5,14 @@ Cubic{BC<:Flag}(::BC) = Cubic{BC}()
 Assuming uniform knots with spacing 1, the `i`th piece of cubic spline
 implemented here is defined as follows.
 
-    y_i(x) = ϕm p(x-i) + ϕ q(x-i) + ϕp q(1- (x-i)) + ϕpp p(1 - (x-i))
+    y_i(x) = cm p(x-i) + c q(x-i) + cp q(1- (x-i)) + cpp p(1 - (x-i))
 
 where
 
-    p(☆) = 1/6 * (1-☆)^3
-    q(☆) = 2/3 - ☆^2 + 1/2 ☆^3
+    p(δx) = 1/6 * (1-δx)^3
+    q(δx) = 2/3 - δx^2 + 1/2 δx^3
 
-and the values `ϕX` for `X ⋹ {m, _, p, pp}` are the pre-filtered coefficients.
+and the values `cX` for `X ∈ {m, _, p, pp}` are the pre-filtered coefficients.
 
 For future reference, this expands out to the following polynomial:
 
@@ -24,6 +24,11 @@ When we derive boundary conditions we will use derivatives `y_0'(x)` and
 """
 Cubic
 
+"""
+`define_indices_d` for a cubic b-spline calculates `ix_d = floor(x_d)` and
+`fx_d = x_d - ix_d` (corresponding to `i` `and `δx` in the docstring for
+`Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`
+"""
 function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
     symix, symixm, symixp = symbol("ix_",d), symbol("ixm_",d), symbol("ixp_",d)
     symixpp, symx, symfx = symbol("ixpp_",d), symbol("x_",d), symbol("fx_",d)
@@ -39,6 +44,14 @@ function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
     end
 end
 
+"""
+`define_indices_d`  for a cubic, periodic b-spline calculates `ix_d = floor(x_d)`
+and `fx_d = x_d - ix_d` (corresponding to `i` and `δx` in the docstring entry
+for `Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`.
+
+If any `ixX_d` for `x ∈ {m, p, pp}` (note: not `c_d`) should fall outside of
+the data interval, they wrap around.
+"""
 function define_indices_d(::Type{BSpline{Cubic{Periodic}}}, d, pad)
     symix, symixm, symixp = symbol("ix_",d), symbol("ixm_",d), symbol("ixp_",d)
     symixpp, symx, symfx = symbol("ixpp_",d), symbol("x_",d), symbol("fx_",d)
@@ -55,15 +68,16 @@ padding{BC<:Flag}(::Type{BSpline{Cubic{BC}}}) = Val{1}()
 padding(::Type{BSpline{Cubic{Periodic}}}) = Val{0}()
 
 """
-In this function we assume that `fx_d = x-ix_d` and we produce `cX_d` for
-`X ⋹ {m, _, p, pp}` such that
+In `coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
+and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
 
     cm_d  = p(fx_d)
     c_d   = q(fx_d)
     cp_d  = q(1-fx_d)
     cpp_d = p(1-fx_d)
 
-where `p` and `q` are defined in the docstring entry for `Cubic`.
+where `p` and `q` are defined in the docstring entry for `Cubic`, and
+`fx_d` in the docstring entry for `define_indices_d`.
 """
 function coefficients{C<:Cubic}(::Type{BSpline{C}}, N, d)
     symm, sym =  symbol("cm_",d), symbol("c_",d)
@@ -82,15 +96,16 @@ function coefficients{C<:Cubic}(::Type{BSpline{C}}, N, d)
 end
 
 """
-In this function we assume that `fx_d = x-ix_d` and we produce `cX_d` for
-`X ⋹ {m, _, p, pp}` such that
+In `gradient_coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
+and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
 
     cm_d  = p'(fx_d)
     c_d   = q'(fx_d)
     cp_d  = q'(1-fx_d)
     cpp_d = p'(1-fx_d)
 
-where `p` and `q` are defined in the docstring for `Cubic`.
+where `p` and `q` are defined in the docstring entry for `Cubic`, and
+`fx_d` in the docstring entry for `define_indices_d`.
 """
 function gradient_coefficients{C<:Cubic}(::Type{BSpline{C}}, d)
     symm, sym, symp, sympp = symbol("cm_",d), symbol("c_",d), symbol("cp_",d), symbol("cpp_",d)
@@ -147,6 +162,14 @@ end
 # Prefiltering #
 # ------------ #
 
+"""
+`Cubic`: continuity in function value, first and second derivatives yields
+
+    2/3 1/6
+    1/6 2/3 1/6
+        1/6 2/3 1/6
+           ⋱  ⋱   ⋱
+"""
 function inner_system_diags{T,C<:Cubic}(::Type{T}, n::Int, ::Type{C})
     du = fill(convert(T, SimpleRatio(1, 6)), n-1)
     d = fill(convert(T, SimpleRatio(2, 3)), n)
@@ -155,8 +178,8 @@ function inner_system_diags{T,C<:Cubic}(::Type{T}, n::Int, ::Type{C})
 end
 
 """
-`Flat` `OnGrid` amounts to setting `y_0'(x) = 0` at `x = 0`. Applying this
-condition gives:
+`Cubic{Flat}` `OnGrid` amounts to setting `y_1'(x) = 0` at `x = 1`.
+Applying this condition yields
 
     -cm + cp = 0
 """
@@ -173,10 +196,18 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
 end
 
 """
-`Flat` `OnCell` amounts to setting `y_0'(x) = 0` at `x = -1/2`. Applying this
-condition gives:
+`Cubic{Flat}`, `OnCell` amounts to setting `y_1'(x) = 0` at `x = 1/2`.
+Applying this condition yields
+
+    -9/8 cm + 11/8 c - 3/8 cp + 1/8 cpp = 0
+
+or, equivalently,
 
     -9 cm + 11 c -3 cp + 1 cpp = 0
+
+(Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
+were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
+close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
 """
 function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                    ::Type{Cubic{Flat}}, ::Type{OnCell})
@@ -198,10 +229,14 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
 end
 
 """
-`Line` `OnCell` amounts to setting `y_0''(x) = 0` at `x = -1/2`. Applying this
-condition gives:
+`Cubic{Line}` `OnCell` amounts to setting `y_1''(x) = 0` at `x = 1/2`.
+Applying this condition yields
 
     3 cm -7 c + 5 cp -1 cpp = 0
+
+(Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
+were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
+close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
 """
 function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                    ::Type{Cubic{Line}}, ::Type{OnCell})
@@ -223,13 +258,10 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
 end
 
 """
-`Line` `OnGrid` amounts to setting `y_0''(x) = 0` at `x = 0`. Applying this
+`Cubic{Line}` `OnGrid` amounts to setting `y_1''(x) = 0` at `x = 1`. Applying this
 condition gives:
 
-    1 cm -2 c + 1 cp = 0 = 0
-
-This is the same system as `Quadratic{Line}`, `OnGrid` so we reuse the
-implementation
+    1 cm -2 c + 1 cp = 0
 """
 function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                    ::Type{Cubic{Line}}, ::Type{OnGrid})
@@ -247,6 +279,14 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
     Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
 end
 
+"""
+`Cubic{Periodic}` `OnGrid` closes the system by looking at the coefficients themselves
+as periodic, yielding
+
+    c0 = c(N+1)
+
+where `N` is the number of data points.
+"""
 function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                                 ::Type{Cubic{Periodic}}, ::Type{GT})
     dl, d, du = inner_system_diags(T,n,Cubic{Periodic})
@@ -260,14 +300,11 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
 end
 
 """
-The free boundary condition for either `OnGrid` or `OnCell` makes sure the
-interpoland has a continuous third derivative at the second-to-outermost cell
-boundary: `y_0'''(1) = y_1'''(1)` and `y_{n-1}'''(n) = y_n'''(n)`. Applying this
-condition gives:
+`Cubic{Free}` `OnGrid` and `Cubic{Free}` `OnCell` amount to requiring an extra
+continuous derivative at the second-to-last cell boundary; this means
+`y_1'''(2) = y_2'''(2)`, yielding
 
     1 cm -3 c + 3 cp -1 cpp = 0
-
-This is the same system as `Quadratic{Free}` so we reuse the implementation
 """
 function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                                 ::Type{Cubic{Free}}, ::Type{GT})

src/b-splines/linear.jl

@@ -1,5 +1,31 @@
 immutable Linear <: Degree{1} end
 
+"""
+Assuming uniform knots with spacing 1, the `i`th peice of linear b-spline
+implemented here is defined as follows.
+
+    y_i(x) = c p(x) + cp p(1-x)
+
+where
+
+    p(δx) = x
+
+and the values `cX` for `X ∈ {_, p}` are the coefficients.
+
+Linear b-splines are naturally interpolating, and require no prefiltering;
+there is therefore no need for boundary conditions to be provided.
+
+Also, although the implementation is slightly different in order to re-use
+the framework built for general b-splines, the resulting interpolant is just
+a piecewise linear function connecting each pair of neighboring data points.
+"""
+Linear
+
+"""
+`define_indices_d` for a linear b-spline calculates `ix_d = floor(x_d)` and
+`fx_d = x_d - ix_d` (corresponding to `i` and `δx` in the docstring for
+`Linear`), as well as the auxiliary quantity `ixp_d`
+"""
 function define_indices_d(::Type{BSpline{Linear}}, d, pad)
     symix, symixp, symfx, symx = symbol("ix_",d), symbol("ixp_",d), symbol("fx_",d), symbol("x_",d)
     quote
@@ -9,6 +35,16 @@ function define_indices_d(::Type{BSpline{Linear}}, d, pad)
     end
 end
 
+"""
+In `coefficients` for a linear b-spline we assume that `fx_d = x-ix-d` and
+we define `cX_d` for `X ∈ {_, p}` such that
+
+    c_d = p(fx_d)
+    cp_d = p(1-fx_d)
+
+where `p` is defined in the docstring entry for `Linear` and `fx_d` in the
+docstring entry for `define_indices_d`.
+"""
 function coefficients(::Type{BSpline{Linear}}, N, d)
     sym, symp, symfx = symbol("c_",d), symbol("cp_",d), symbol("fx_",d)
     quote
@@ -17,6 +53,16 @@ function coefficients(::Type{BSpline{Linear}}, N, d)
     end
 end
 
+"""
+In `gradient_coefficients` for a linear b-spline we assume that `fx_d = x-ix_d`
+and we define `cX_d` for `X ⋹ {_, p}` such that
+
+    c_d  = p'(fx_d)
+    cp_d = p'(1-fx_d)
+
+where `p` is defined in the docstring entry for `Linear`, and `fx_d` in the
+docstring entry for `define_indices_d`.
+"""
 function gradient_coefficients(::Type{BSpline{Linear}}, d)
     sym, symp, symfx = symbol("c_",d), symbol("cp_",d), symbol("fx_",d)
     quote
@@ -25,6 +71,16 @@ function gradient_coefficients(::Type{BSpline{Linear}}, d)
     end
 end
 
+"""
+In `hessian_coefficients` for a linear b-spline we assume that `fx_d = x-ix_d`
+and we define `cX_d` for `X ⋹ {_, p}` such that
+
+    c_d  = p''(fx_d)
+    cp_d = p''(1-fx_d)
+
+where `p` is defined in the docstring entry for `Linear`, and `fx_d` in the
+docstring entry for `define_indices_d`. (These are both ≡ 0.)
+"""
 function hessian_coefficients(::Type{BSpline{Linear}}, d)
     sym, symp = symbol("c_",d), symbol("cp_",d)
     quote