diff --git a/Project.toml b/Project.toml index 3a68b98..4763971 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CompactBases" uuid = "2c0377a8-7469-4ebd-be0f-82e501f20078" authors = ["Stefanos Carlström "] -version = "0.3.12" +version = "0.3.13" [deps] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" diff --git a/docs/make.jl b/docs/make.jl index f60cfa1..dffecec 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -35,6 +35,7 @@ makedocs(; ] ], ], + "Basis transforms" => "basis_transforms.md" ], repo="https://github.com/JuliaApproximation/CompactBases.jl/blob/{commit}{path}#L{line}", sitename="CompactBases.jl", diff --git a/docs/plots.jl b/docs/plots.jl index f71fe6d..fcf2de6 100644 --- a/docs/plots.jl +++ b/docs/plots.jl @@ -269,6 +269,75 @@ function diagonal_operators() savedocfig("diagonal_operators") end +function basis_transforms() + function compare_bases(x, bases, functions) + m = length(bases) + n = length(functions) + + cs = [B \ f.(axes(B,1)) + for B in bases, f in functions] + es = extrema.(cs) + ymins = [minimum(first.(es[i,:])) for i = 1:m] + ymaxs = [maximum(last.(es[i,:])) for i = 1:m] + + cfigure("basis comparison", figsize=(10,16)) do + for i = 1:m + B = bases[i] + a = axes(B,1) + χ = B[x,:] + xl = [mean_position(x, view(χ, :, j)) for j in 1:size(B,2)] + csubplot(m,n+1,(i,1), nox=i exp(-(x-1.0).^2/(2*0.2^2)) + u = x -> (0.5 < x < 1.5)*one(x) + + x = range(-1.5, stop=3.5, length=1001) + + compare_bases(x, (A,B,C,D,E), (gauss,u)) + + savedocfig("basis_transforms") +end + macro echo(expr) println(expr) :(@time $expr) @@ -285,3 +354,4 @@ include("bspline_plots.jl") include("fd_plots.jl") @echo densities() @echo diagonal_operators() +@echo basis_transforms() diff --git a/docs/src/basis_transforms.md b/docs/src/basis_transforms.md new file mode 100644 index 0000000..eb91ec2 --- /dev/null +++ b/docs/src/basis_transforms.md @@ -0,0 +1,221 @@ +# Transforms between different bases + +Sometimes, we have a function ``f`` expanded over the basis functions +of one basis ``A``, +```math +f_A = A\vec{c} +``` +and wish to re-express as an expansion over the +basis functions of another basis ``B``, +```math +f_B = B\vec{d}. +``` +We thus want to solve ``f_A = f_B`` for ``\vec{d}``; this is very +easy, and follows the same reasoning as in [Solving +equations](@ref). We begin by projecting into the space ``B`` by using +the projector ``BB^H``: +```math +BB^HA\vec{c} = +BB^HB\vec{d}. +``` +This has to hold for any ``B``, as long as the basis functions of +``B`` are linearly independent, and is then equivalent to +```math +B^HA\vec{c} = +B^HB\vec{d}, +``` +the solution of which is +```math +\vec{d} = +(B^HB)^{-1} +B^HA +\vec{c} +\defd +S_{BB}^{-1} +S_{BA} +\vec{c}. +``` + +[`basis_overlap`](@ref) computes ``S_{BA}=B^HA`` and +[`basis_transform`](@ref) the full transformation matrix +``(B^HB)^{-1}B^HA``; sometimes it may be desirable to perform the +transformation manually in two steps, since the combined matrix may be +dense whereas the constituents may be relatively sparse. For +orthogonal bases, ``B^HB`` is naturally diagonal. + +There is no requirement that the two bases span the same intervals, +i.e. `axes(A,1)` and `axes(B,1)` need not be fully overlapping. Of +course, if they are fully disjoint, the transform matrix will be +zero. Additionally, some functions may not be well represented in a +particular basis, e.g. the step function is not well approximated by +B-splines (of order ``>1``), but poses no problem for +finite-differences, except at the discontinuity. To judge if a +transform is faithful thus amount to comparing the transformed +expansion coefficients with those obtained by expanding the function +in the target basis directly, instead of comparing reconstruction +errors. + +## Examples + +We will now illustrate transformation between B-splines, FEDVR, and +various finite-differences: +```julia +julia> A = BSpline(LinearKnotSet(7, 0, 2, 11)) +BSpline{Float64} basis with typename(LinearKnotSet)(Float64) of order k = 7 on 0.0..2.0 (11 intervals) + +julia> B = FEDVR(range(0.0, 2.5, 5), 6)[:,2:end-1] +FEDVR{Float64} basis with 4 elements on 0.0..2.5, restricted to elements 1:4, basis functions 2..20 ⊂ 1..21 + +julia> C = StaggeredFiniteDifferences(0.03, 0.3, 0.01, 3.0) +Staggered finite differences basis {Float64} on 0.0..3.066973760326056 with 90 points with spacing varying from 0.030040496962651875 to 0.037956283408020486 + +julia> D = FiniteDifferences(range(0.25,1.75,length=31)) +Finite differences basis {Float64} on 0.2..1.8 with 31 points spaced by Δx = 0.05 + +julia> E = BSpline(LinearKnotSet(9, -1, 3, 15))[:,2:end] +BSpline{Float64} basis with typename(LinearKnotSet)(Float64) of order k = 9 on -1.0..3.0 (15 intervals), restricted to basis functions 2..23 ⊂ 1..23 +``` + +We use these bases to expand two test functions: +```math +f_1(x) = \exp\left[ +-\frac{(x-1)^2}{2\times0.2^2} +\right] +``` +and +```math +f_2(x) = \theta(x-0.5) - \theta(x-1.5), +``` +where ``\theta(x)`` is the [Heaviside step +function](https://en.wikipedia.org/wiki/Heaviside_step_function). + +```julia +julia> gauss(x) = exp(-(x-1.0).^2/(2*0.2^2)) +gauss (generic function with 1 method) + +julia> u(x) = (0.5 < x < 1.5)*one(x) +u (generic function with 1 method) +``` + +If for example, we first expand ``f_1`` in B-splines (``A``), and then wish to +project onto the staggered finite-differences (``C``), we do the +following: + +```julia +julia> c = A \ gauss.(axes(A, 1)) +17-element Vector{Float64}: + -0.0003757237742430358 + 0.001655433792947186 + -0.003950561899347059 + 0.00714774775995981 + -0.011051108016208545 + 0.0173792655408227 + 0.04090160980771093 + 0.6336420285181446 + 1.3884690430556483 + 0.6336420285181449 + 0.04090160980771065 + 0.01737926554082302 + -0.011051108016208311 + 0.007147747759958865 + -0.00395056189934631 + 0.0016554337929467365 + -0.00037572377424269145 + +julia> S_CA = basis_transform(C, A) +90×17 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1530 stored entries: +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⣿⣿⣿⣿⣿⣿⣿⣿⡇ +⠛⠛⠛⠛⠛⠛⠛⠛⠃ + +julia> S_CA*c +90-element Vector{Float64}: + 3.805531055898412e-5 + 1.6890846360940058e-5 + -4.5403814886370545e-5 + -4.121425310424971e-5 + 1.2245471085512564e-5 + 6.482348697494439e-5 + 8.896755347867799e-5 + 9.703404969825733e-5 + 0.0001290654548866452 + 0.00023543152355538763 + 0.0004663267687150117 + ⋮ + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + +julia> basis_transform(C, A, c) # Or transform without forming matrix +90-element Vector{Float64}: + 3.80553105589841e-5 + 1.6890846360940065e-5 + -4.540381488637052e-5 + -4.121425310424972e-5 + 1.2245471085512601e-5 + 6.482348697494439e-5 + 8.896755347867799e-5 + 9.703404969825729e-5 + 0.00012906545488664526 + 0.0002354315235553877 + 0.00046632676871501176 + ⋮ + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 +``` + +In the plots below, we show each basis (labelled ``A-E``), and their +respective reconstructions of ``f_1`` and ``f_2`` as solid red lines, +with the true functions shown as dashed black lines. For each +basis–function combination, shown are also the expansion coefficients +as red points joined by solid lines, as well as the expansion +coefficients transformed from the other bases, as labelled in the +legends. As an example, we see that first expanding in +finite-differences ``D`` and then transforming to B-splines ``A`` is +numerically unstable, leading to very large coefficients at the edges +of ``D``. Going _to_ finite-differences is always stable, even though +the results may not be accurate. + +![Transforms between bases](../figures/basis_transforms.svg) + +## Reference + +```@docs +basis_overlap +basis_transform +``` diff --git a/docs/src/overview.md b/docs/src/overview.md index d20713d..8558169 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -74,6 +74,9 @@ julia> S = B'B ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1 + +julia> orthogonality(B) +CompactBases.Orthogonal() ``` In contrast, a non-orthogonal basis such as B-splines has a _banded_ @@ -111,6 +114,9 @@ julia> S = B'B ⋅ ⋅ 9.79816e-7 0.00036737 0.0100953 0.0349662 0.0464947 0.0246054 0.00347002 ⋅ ⋅ ⋅ 1.65344e-6 0.000811839 0.00747795 0.0246054 0.0331746 0.0139286 ⋅ ⋅ ⋅ ⋅ 1.32275e-5 0.000365961 0.00347002 0.0139286 0.0222222 + +julia> orthogonality(B) +CompactBases.NonOrthogonal() ``` where the bandwidth depends on the B-spline order. diff --git a/src/CompactBases.jl b/src/CompactBases.jl index 01c6787..efc4d17 100644 --- a/src/CompactBases.jl +++ b/src/CompactBases.jl @@ -50,6 +50,7 @@ using RecipesBase vandermonde(B::Basis) = B[locs(B),:] include("distributions.jl") +include("orthogonality.jl") include("restricted_bases.jl") include("types.jl") include("unit_vectors.jl") @@ -78,6 +79,8 @@ include("inner_products.jl") include("densities.jl") include("linear_operators.jl") +include("basis_transforms.jl") + """ centers(B) diff --git a/src/basis_transforms.jl b/src/basis_transforms.jl new file mode 100644 index 0000000..bf80df8 --- /dev/null +++ b/src/basis_transforms.jl @@ -0,0 +1,181 @@ +# # * General case + +@doc raw""" + basis_transform(A, B) + +Compute the (possibly dense) basis transform matrix to go from `B` to +`A`, via [`basis_overlap`](@ref) for orthogonal basis `A`: + +```math +B^HA +``` +""" +basis_transform(::Ortho, A, B) = basis_overlap(A, B) + +@doc raw""" + basis_transform(A, B) + +Compute the (possibly dense) basis transform matrix to go from `B` to +`A`, via [`basis_overlap`](@ref) for non-orthogonal basis `A`: + +```math +(B^HB)^{-1}B^HA +``` +""" +basis_transform(::NonOrthogonal, A, B) = (A'A) \ basis_overlap(A, B) + +""" + basis_transform(A, B) + +Compute the (possibly dense) basis transform matrix to go from `B` to +`A`, via [`basis_overlap`](@ref), dispatching on the +[`orthogonality`](@ref) trait of `A`. +""" +basis_transform(A, B) = basis_transform(orthogonality(A), A, B) + +# * To finite-differences + +""" + basis_overlap(A::FiniteDifferencesOrRestricted, B) + +Transforming to finite-differences is simple: just evaluate the basis +functions of the source basis on the grid points, possibly weighting +them. +""" +function basis_overlap(A::FiniteDifferencesOrRestricted, B) + xa = locs(A) + wa = weights(A) + + B[xa,:] ./ wa +end + +function basis_transform(A::FiniteDifferencesOrRestricted, B, c::AbstractVector) + xa = locs(A) + wa = weights(A) + + (B[xa,:]*c) ./ wa +end + +# * To FEDVR + +""" + basis_overlap(A::FEDVROrRestricted, B) + +Transforming to FEDVR is much the same as for finite-differences; +evaluate at the quadrature nodes and scale by the weights. However, +this has to be done per-element, so we feed it through the +interpolation routines already setup for this. +""" +function basis_overlap(A::FEDVROrRestricted, B) + T = promote_type(eltype(A), eltype(B)) + S = Matrix{T}(undef, size(A,2), size(B,2)) + + xa = axes(A, 1) + for j ∈ axes(B,2) + # This is ugly; it is necessary since B[x,j] checks that x ⊆ + # axes(B, 1), whereas B[[x],j:j] does not, and for the purposes of the + # transform, we wish to evaluate the basis functions as vanishing + # outside the axis, which the latter does. + f = x -> first(B[[x],j:j]) + S[:,j] .= A \ f.(xa) + end + + S +end + +function basis_transform(A::FEDVROrRestricted, B, c::AbstractVector) + xa = axes(A, 1) + # This is ugly; it is necessary since B[x,:] checks that x ⊆ + # axes(B, 1), whereas B[[x],:] does not, and for the purposes of the + # transform, we wish to evaluate the basis functions as vanishing + # outside the axis, which the latter does. + f = x -> dot(B[[x],:],c) + A \ f.(xa) +end + +# * To B-splines +# ** From polynomial bases + +function _basis_transform_2_bspline_common(A::BSplineOrRestricted, + B::Union{BSplineOrRestricted,FEDVROrRestricted}) + k = order(A) + k′ = max(k, maximum(order(B))) + N = num_quadrature_points(k, k′-k+2) + + x,w = lgwt(knotset(A), N) + W = Diagonal(w) + + x, A[x,:]'W +end + +""" + basis_overlap(A::BSplineOrRestricted, B::Union{BSplineOrRestricted,FEDVROrRestricted}) + +From a piecewise polynomial basis, i.e. B-splines or FEDVR, we can +use Gaussian quadrature to compute the overlap matrix elements +exactly. It is possible that we could derive better results, if we +used the fact that FEDVR basis functions are orthogonal in the sense +of the Gauss–Lobatto quadrature, but we leave that for later. +""" +function basis_overlap(A::BSplineOrRestricted, B::Union{BSplineOrRestricted,FEDVROrRestricted}) + T = promote_type(eltype(A), eltype(B)) + S = Matrix{T}(undef, size(A,2), size(B,2)) + + x,χA = _basis_transform_2_bspline_common(A, B) + + for j ∈ axes(B,2) + S[:,j] .= χA*B[x,j] + end + + S +end + +function basis_transform(A::BSplineOrRestricted, B::Union{BSplineOrRestricted,FEDVROrRestricted}, + c::AbstractVector) + T = promote_type(eltype(A), eltype(B), eltype(c)) + d = zeros(T, size(A, 2)) + + x,χA = _basis_transform_2_bspline_common(A, B) + + for j ∈ axes(B,2) + d .+= c[j] .* (χA*B[x,j]) + end + + (A'A) \ d +end + +# ** From finite-differences (or any other basis) + +""" + basis_transform(A::BSplineOrRestricted, B::AbstractFiniteDifferences) + +This works via Vandermonde interpolation, which is potentially +unstable. In this case, we do not provide [`basis_overlap`](@ref). +""" +function basis_transform(A::BSplineOrRestricted, B::AbstractFiniteDifferences) + T = promote_type(eltype(A), eltype(B)) + S = Matrix{T}(undef, size(A,2), size(B,2)) + + x = locs(B) + χA = A[x,:] + for j ∈ axes(B,2) + S[:,j] .= χA \ Vector(B[x,j]) + end + + S +end + +function basis_transform(A::BSplineOrRestricted, B::AbstractFiniteDifferences, c::AbstractVector) + T = promote_type(eltype(A), eltype(B)) + N = size(B,2) + + x = locs(B) + χA = A[x,:] + χB = B[x,:] + + χA \ (χB*c) +end + +# * Exports + +export basis_transform, basis_overlap diff --git a/src/bsplines.jl b/src/bsplines.jl index 615075b..bfa34cf 100644 --- a/src/bsplines.jl +++ b/src/bsplines.jl @@ -127,7 +127,9 @@ Base.hash(B::BSpline, h::UInt) = hash(B.t, h) distribution(B::BSpline) = distribution(B.t) -order(B::BSplineOrRestricted) = order(unrestricted_basis(B).t) +knotset(B::BSpline) = B.t +knotset(B::RestrictedBSpline) = knotset(unrestricted_basis(B)) +order(B::BSplineOrRestricted) = order(knotset(B)) function show(io::IO, B::BSpline{T}) where T write(io, "BSpline{$(T)} basis with $(B.t)") diff --git a/src/fedvr.jl b/src/fedvr.jl index e31ebad..069efb4 100644 --- a/src/fedvr.jl +++ b/src/fedvr.jl @@ -72,6 +72,7 @@ assert_compatible_bases(A::FEDVROrRestricted, B::FEDVROrRestricted) = parent(A) == parent(B) || throw(ArgumentError("Can only multiply FEDVRs of the same order and sharing the same nodes")) +orthogonality( ::FEDVROrRestricted) = Orthonormal() distribution(B::FEDVR) = distribution(B.t) order(B::FEDVR) = B.order diff --git a/src/finite_differences.jl b/src/finite_differences.jl index 25132b4..4738bbd 100644 --- a/src/finite_differences.jl +++ b/src/finite_differences.jl @@ -11,6 +11,11 @@ size(B::AbstractFiniteDifferences) = (ℵ₁, length(locs(B))) distribution(B::AbstractFiniteDifferences) = distribution(locs(B)) +orthogonality(::Uniform, ::AbstractFiniteDifferences) = Orthogonal() +orthogonality(::NonUniform, ::AbstractFiniteDifferences) = Orthonormal() +orthogonality(B::AbstractFiniteDifferences) = orthogonality(distribution(B), B) +orthogonality(B::RestrictedFiniteDifferences) = orthogonality(parent(B)) + # step is only well-defined for uniform grids step(::Uniform, B::AbstractFiniteDifferences) = step(locs(B)) # All-the-same, we define step for non-uniform grids to make the mass diff --git a/src/orthogonality.jl b/src/orthogonality.jl new file mode 100644 index 0000000..38efa2d --- /dev/null +++ b/src/orthogonality.jl @@ -0,0 +1,59 @@ +""" + BasisOrthogonality + +Base type for the orthogonality trait of a basis. +""" +abstract type BasisOrthogonality end + +""" + Ortho <: BasisOrthogonality + +Intermediate type for the trait for all orthogonal bases. +""" +abstract type Ortho <: BasisOrthogonality end + +""" + Orthogonal <: Ortho + +Concrete type for the trait of orthogonal bases, i.e. with diagonal metric. +""" +struct Orthogonal <: Ortho end + +""" + Orthonormal <: Ortho + +Concrete type for the trait of orthonormal bases, i.e. with metric `I`. +""" +struct Orthonormal <: Ortho end + +""" + NonOrthogonal <: BasisOrthogonality + +Concrete type for the trait of non-orthogonal bases, i.e. with non-diagonal metric. +""" +struct NonOrthogonal <: BasisOrthogonality end + +""" + orthogonality(B) + +Returns the [`BasisOrthogonality`](@ref) trait of the basis `B`. + +# Examples + +```julia +julia> orthogonality(BSpline(LinearKnotSet(7, 0, 1, 11))) +CompactBases.NonOrthogonal() + +julia> orthogonality(FEDVR(range(0, stop=1, length=11), 7)) +CompactBases.Orthonormal() + +julia> orthogonality(StaggeredFiniteDifferences(0.1, 0.3, 0.1, 10.0)) +CompactBases.Orthonormal() + +julia> orthogonality(FiniteDifferences(range(0, stop=1, length=11))) +CompactBases.Orthogonal() +``` +""" +function orthogonality end + +export orthogonality diff --git a/src/restricted_bases.jl b/src/restricted_bases.jl index e3adc61..a80b0eb 100644 --- a/src/restricted_bases.jl +++ b/src/restricted_bases.jl @@ -65,6 +65,8 @@ IntervalSets.leftendpoint(B::RestrictedQuasiArray) = IntervalSets.rightendpoint(B::RestrictedQuasiArray) = rightendpoint(parent(B)) +orthogonality(::BasisOrRestricted) = NonOrthogonal() + # TODO: This is invalid for non-orthogonal bases such as B-splines, # since it selects as many quadrature nodes as there are basis # functions in the restriction. Each B-spline is actually associated diff --git a/test/basis_transforms.jl b/test/basis_transforms.jl new file mode 100644 index 0000000..d6cb4df --- /dev/null +++ b/test/basis_transforms.jl @@ -0,0 +1,57 @@ +@testset "Basis transforms" begin + A = BSpline(LinearKnotSet(7, 0, 2, 11)) + A′ = A[:,2:end-1] + B = FEDVR(range(0.0, stop=2.5, length=5), 6)[:,2:end-1] + C = StaggeredFiniteDifferences(0.03, 0.3, 0.01, 3.0) + D = FiniteDifferences(range(0.25, stop=1.75, length=31)) + E = BSpline(LinearKnotSet(9, -1, 3, 15))[:,2:end] + F = FiniteDifferences(range(1, stop=4, step=0.1)) + + ex(B,f) = B \ f.(axes(B,1)) + bt(A,B,f) = basis_transform(A,B,ex(B,f)) + bt2(A,B,f) = basis_transform(A,B)*ex(B,f) + + gauss = x -> exp(-(x-1.0).^2/(2*0.2^2)) + # u = x -> (0.5 < x < 1.5)*one(x) + + for (A,B,f,rtol) in [(A,A,gauss,1e-13), + (A′,A,gauss,1e-3), + (B,A,gauss,5e-3), + (C,A,gauss,5e-3), + (D,A,gauss,5e-3), + (E,A,gauss,1e-1), + (F,A,gauss,5e-3), + + (A,A′,gauss,1e-3), + (A′,A′,gauss,1e-13), + (B,A′,gauss,5e-3), + (C,A′,gauss,5e-3), + (D,A′,gauss,5e-3), + (E,A′,gauss,1e-1), + (F,A′,gauss,5e-3), + + (A,B,gauss,1e-1), + (A′,B,gauss,1e-1), + (B,B,gauss,1e-13), + (C,B,gauss,1e-2), + (D,B,gauss,1e-2), + (E,B,gauss,1e-1), + (F,B,gauss,1e-2), + + (A,C,gauss,1e-2), + (A′,C,gauss,5e-3), + (B,C,gauss,1e-2), + (C,C,gauss,1e-13), + (D,C,gauss,1e-2), + (F,C,gauss,1e-2)] + c = ex(A, f) + d = bt(A, B, f) + d2 = bt2(A, B, f) + + (norm(c-d)/norm(c) ≥ rtol || norm(c-d2)/norm(c) ≥ rtol) && + @error "Basis transform accuracy test failed" A B rtol norm(c-d)/norm(c) norm(c-d)/norm(c) < rtol norm(c-d2)/norm(c) norm(c-d2)/norm(c) < rtol + + @test d ≈ c rtol=rtol + @test d2 ≈ c rtol=rtol + end +end diff --git a/test/orthogonality.jl b/test/orthogonality.jl new file mode 100644 index 0000000..ccd3043 --- /dev/null +++ b/test/orthogonality.jl @@ -0,0 +1,17 @@ +@testset "Orthogonality" begin + @test orthogonality(BSpline(LinearKnotSet(7, 0, 1, 11))) == + CompactBases.NonOrthogonal() + + @test orthogonality(FEDVR(range(0, stop=1, length=11), 7)) == + CompactBases.Orthonormal() + + SFD = StaggeredFiniteDifferences(0.1, 0.3, 0.1, 10.0) + @test orthogonality(SFD) == CompactBases.Orthonormal() + @test orthogonality(SFD[:,2:end-1]) == CompactBases.Orthonormal() + + @test orthogonality(StaggeredFiniteDifferences(11, 0.1)) == + CompactBases.Orthogonal() + + @test orthogonality(FiniteDifferences(range(0, stop=1, length=11))) == + CompactBases.Orthogonal() +end diff --git a/test/runtests.jl b/test/runtests.jl index 5d702e0..28c350d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,4 +39,6 @@ include("derivative_accuracy_utils.jl") include("inner_products.jl") include("densities.jl") include("linear_operators.jl") + include("orthogonality.jl") + include("basis_transforms.jl") end