API extensions, rewrite docs

docs/src/index.md

@@ -24,40 +24,50 @@ In this package,
 
 3. [`basis_at`](@ref) returns an *iterator* for evaluating basis functions at an arbitrary point inside their domain. This iterator is meant to be heavily optimized and non-allocating. [`linear_combination`](@ref) is a convenience wrapper for obtaining a linear combination of basis functions at a given point.
 
-A basis is constructed using
+## Univariate basis example
 
-1. a family on a *fixed* domain, eg [`Chebyshev`](@ref),
+We construct a basis of Chebyshev polynomials on ``[0, 4]``. This requires a transformation since their canonical domain is ``[-1,1]``. Other transformations include [`SemiInfRational`](@ref) (for ``[A, \infty]`` intervals)  and [`InfRational`](@ref).
 
-2. a grid specification like [`InteriorGrid`](@ref),
-
-3. a number of parameters that determine the number of basis functions.
+We display the domian and the dimension (number of basis functions).
+```@example univariate
+basis = Chebyshev(InteriorGrid(), 5) ∘ BoundedLinear(0, 4)
+domain(basis)
+dimension(basis)
+```
 
-A set of coordinates for a particular basis can be augmented for a wider basis with [`augment_coefficients`](@ref).
+We have chosen an interior grid, shown below. We `collect` the result for the purpose of this tutorial, since `grid` returns an iterable to avoid allocations.
+```@example univariate
+collect(grid(basis))
+```
 
-Bases have a “canonical” domain, eg ``[-1,1]`` or ``[-1,1]^n`` for Chebyshev polynomials. Use [transformations](#domains-and-transformations) for mapping to other domains.
+We can show evaluate the basis functions at a given point. Again, it is an iterable, so we `collect` to show it here.
+```@example univariate
+collect(basis_at(basis, 0.41))
+```
 
-## Examples
+We can evaluate linear combination as directly, or via partial application.
+```@example univariate
+θ = [1, 0.5, 0.2, 0.3, 0.001]               # a vector of coefficients
+x = 0.41
+linear_combination(basis, θ, x)          # combination at some value
+linear_combination(basis, θ)(x)          # also as a callable
+```
 
-### Univariate family on `[-1,1]`
+We can also evaluate *derivatives* of either the basis or the linear combination at a given point. Here we want the derivatives up to order 3.
+```@example univariate
+dx = (𝑑^2)(x)
+collect(basis_at(basis, dx))
+linear_combination(basis, θ, x)
+```
 
-```@example
-using SpectralKit
-basis = Chebyshev(EndpointGrid(), 5)        # 5 Chebyshev polynomials
-is_function_basis(basis)                    # ie we support the interface below
-dimension(basis)                            # number of basis functions
-domain(basis)                               # domain
-grid(basis)                                 # Gauss-Lobatto grid
-collect(basis_at(basis, 0.41))              # iterator for basis functions at 0.41
-collect(basis_at(basis, derivatives(0.41))) # values and 1st derivatives
-θ = [1, 0.5, 0.2, 0.3, 0.001]               # a vector of coefficients
-linear_combination(basis, θ, 0.41)          # combination at some value
-linear_combination(basis, θ)(0.41)          # also as a callable
-basis2 = Chebyshev(EndpointGrid(), 8)       # 8 Chebyshev polynomials
+Having an approximation, we can embed it in a larger basis, extending the coefficients accordingly.
+```@example univariate
+basis2 = Chebyshev(EndpointGrid(), 8) ∘ transformation(basis)       # 8 Chebyshev polynomials
 is_subset_basis(basis, basis2)              # we could augment θ …
 augment_coefficients(basis, basis2, θ)      # … so let's do it
 ```
 
-### Smolyak approximation on a transformed domain
+## Multivariate (Smolyak) approximation example
 
 ```@example
 using SpectralKit, StaticArrays

src/generic_api.jl

@@ -4,7 +4,7 @@
 
 export is_function_basis, dimension, basis_at, linear_combination, InteriorGrid,
     InteriorGrid2, EndpointGrid, grid, collocation_matrix, augment_coefficients,
-    is_subset_basis
+    is_subset_basis, transformation
 
 """
 $(TYPEDEF)
@@ -273,6 +273,10 @@ function Base.:(∘)(parent::FunctionBasis, transformation)
     TransformedBasis(parent, transformation)
 end
 
+Base.parent(basis::TransformedBasis) = basis.parent
+
+transformation(basis::TransformedBasis) = basis.transformation
+
 domain(basis::TransformedBasis) = domain(basis.transformation)
 
 dimension(basis::TransformedBasis) = dimension(basis.parent)
@@ -300,7 +304,15 @@ function Base.getindex(basis::TransformedBasis{<:MultivariateBasis}, i::Int)
     TransformedBasis(parent[i], Tuple(transformation)[i])
 end
 
-# FIXME add augmentation for transformed bases
+function is_subset_basis(basis1::TransformedBasis, basis2::TransformedBasis)
+    basis1.transformation ≡ basis2.transformation &&
+        is_subset_basis(basis1.parent, basis2.parent)
+end
+
+function augment_coefficients(basis1::TransformedBasis, basis2::TransformedBasis, θ1)
+    @argcheck is_subset_basis(basis1, basis2)
+    augment_coefficients(basis1.parent, basis2.parent, θ1)
+end
 
 function transform_to(basis::FunctionBasis, transformation, x)
     transform_to(domain(basis), transformation, x)

test/test_chebyshev.jl

@@ -73,6 +73,22 @@
     end
 end
 
+@testset "augmentation of transformed basis" begin
+    N = 5
+    M = N + 3
+    t = SemiInfRational(0.3, 0.9)
+    grid_kind = InteriorGrid()
+    basis =  Chebyshev(grid_kind, N) ∘ t
+    basis′ =  Chebyshev(grid_kind, M) ∘ t
+    @test is_subset_basis(basis, basis′)
+    for _ in 1:100
+        x = rand_in_domain(basis)
+        θ = rand(N)
+        θ′ = augment_coefficients(basis, basis′, θ)
+        @test linear_combination(basis, θ, x) ≈ linear_combination(basis′, θ′, x)
+    end
+end
+
 @testset "univariate derivatives" begin
     basis = Chebyshev(InteriorGrid(), 5)
     for (transformation, N) in ((BoundedLinear(-2, 3), 5),

test/test_generic_api.jl

@@ -35,6 +35,8 @@ end
     t = BoundedLinear(1.0, 2.0)
     @test domain(basis ∘ t) == domain(t)
     @test dimension(basis ∘ t) == dimension(basis)
+    @test parent(basis ∘ t) ≡ basis
+    @test transformation(basis ∘ t) ≡ t
     @test collect(grid(basis ∘ t)) ==
         [transform_from(domain(basis), t, x) for x in grid(basis)]
 

test/utilities.jl

@@ -42,7 +42,7 @@ function rand_in_domain(basis::SmolyakBasis{<:SmolyakIndices{N}}) where N
 end
 
 function rand_in_domain(basis::TransformedBasis)
-    (; parent, transformation)
+    (; parent, transformation) = basis
     transform_from(parent, transformation, rand_in_domain(parent))
 end