Support algebra operations on expansions

[dlfivefifty]

@jagot FYI, this will be a big change so let me know your thoughts.

At the moment it's just a quick prototype: I want to make it trait based, so that, for example, view(f, 0.5:0.5:1) is also treated as an expansion. I think this is doable via a BroadcastStyle, instead of overloading broadcasted directly.

[jagot]

I think it looks a lot like the ad hoc thing I have in Atoms.jl: https://github.com/JuliaAtoms/Atoms.jl/blob/compact-bases/src/radial_orbitals.jl#L2 (and also in CompactBases.jl: https://github.com/JuliaApproximation/CompactBases.jl/blob/master/src/types.jl), so I think it is a good idea to upstream the type definitions to ContinuumArrays.jl. Maybe also add the higher-order tensor case, i.e. const Expansion{T,N,Space<:Basis,Coeffs<:AbstractArray{T,N}} = ApplyQuasiVector{T,typeof(*),<:Tuple{Space,Coeffs}}? This would represent a 1-d curve in N-d space, if I'm thinking correctly.

The algebra is also very nice to have! Is D^2 different from D'D? I would assume so, except for Dirichlet-0 boundary conditions.

Will the comparison of Expansions work with strict equality?

ContinuumArrays.jl/src/bases/bases.jl

Lines 192 to 197 in 0c94353

	@eval function ==(f::Expansion, g::Expansion)
	S,c = arguments(f)
	T,d = arguments(g)
	ST = broadcastbasis(+, S, T)
	(ST \ S) * c == (ST \ T) * d
	end

I.e. what type should broadcastbasis return?

And what is the difference between + and - in this case?:

ContinuumArrays.jl/src/bases/bases.jl

Line 181 in 0c94353

broadcastbasis(::typeof(-), a, b) = broadcastbasis(+, a, b)

[dlfivefifty]

The algebra is also very nice to have! Is D^2 different from D'D?

Yes. D^2 doesn't really make sense for low order FEM bases: it's more appropriate for Fourier and Chebyshev.

Will the comparison of Expansions work with strict equality?

Yes: if ST = T then ST \T is Eye(∞).

And what is the difference between + and - in this case?

The idea was that one would only need to overload one of these.

[dlfivefifty]

Maybe also add the higher-order tensor case

I don't follow: higher dimensions is a property of the axes. Unless you are proposing we support tensor algebra, that is, QuasiArray{T,3} can act on a QuasiMatrix{T}? That would require a lot of thought.

[jagot]

Yeah, I realized I was wrong; rather I meant supporting

const Expansion{T,N,Space<:Basis,Coeffs<:AbstractVecOrMat{T}} =
    ApplyQuasiVector{T,typeof(*),<:Tuple{Space,Coeffs}}

Then, if Coeffs<:AbstractVector, it would correspond to a single function expansion, whereas if Coeffs<:AbstractMatrix, each column of the matrix could be taken to be a separate function expansion (which is how I typically use it), or indeed a 1d curve in N-d space. Basically like the last figure here: https://juliaapproximation.github.io/CompactBases.jl/stable/bsplines/splines/

I think for functions in higher dimensions (this is really #15), i.e. beyond curves, you need to increase the number of quasi-axes as you increase the number of discrete axes (if that makes sense). I.e. a function mapping from 2d -> 1d would have 2 quasi-axes and 2 "basis function" axes (which we could denote e.g. (2,2)->1). These could either be constructed from tensor products of two (1,1)->1 quasi-matrices, or a basis that is inherently (2,2)->1, such as RBFs.

Also mentioned in #15 is the case (1,1)->2, i.e a vector-valued function.

[dlfivefifty]

If the coefficients are a matrix, then the basis needs to be a 3-tensor, with 3-tensor times matrix multiplication defined. otherwise * doesn't make sense.

[jagot]

Well, then I would simply have Expansions as a convenience alias where the coefficients are represented by a matrix, simply meaning a collection of function expansions.

EDIT: The reason being if C is the matrix of expansion coeffs, and χ = R[x,:], χ*C will give me a matrix where each column corresponds to a reconstructed function.

[dlfivefifty]

Algebra has to make sense: if A is a Basis <: QuasiMatrix and C is a Matrix, then A*C is a QuasiMatrix with axes axes(A*C) == (axes(A,1), axes(C,2)). So yes it sounds fine, provided that each column is A*C[:,j].

Where things are more complicated were if you wanted spherical harmonics with the coefficients represented as a matrix. This would not makes sense as the basis is different for each column. So in that case it would be best to wrap the matrix as a block-vector and avoid the need to define tensors.

[jagot]

Fair enough. For spherical harmonics, the number of coefficients would be different though, for each irreducible representation? I.e. for ell=0,1,2,... the dimensions would be 1,3,5,.... I guess in addition to storing the coefficients as a (column) block-vector, you would need to store the bases for each irrep in a (row) block-vector?

[dlfivefifty]

Just as a block-quasi-matrix, see https://github.com/JuliaApproximation/SphericalHarmonics.jl

[codecov]

Codecov Report

Merging #54 into master will increase coverage by 9.79%.
The diff coverage is 74.60%.

@@            Coverage Diff             @@
##           master      #54      +/-   ##
==========================================
+ Coverage   69.64%   79.43%   +9.79%     
==========================================
  Files           4        4              
  Lines         336      355      +19     
==========================================
+ Hits          234      282      +48     
+ Misses        102       73      -29     

Impacted Files	Coverage Δ
src/operators.jl	69.64% <42.85%> (+27.11%)	⬆️
src/bases/bases.jl	71.73% <64.28%> (-0.59%)	⬇️
src/bases/splines.jl	90.27% <85.71%> (-5.04%)	⬇️
src/ContinuumArrays.jl	88.76% <100.00%> (+13.42%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c0d13f4...c591fd0. Read the comment docs.