Use ContinuumArrays.jl

[dlfivefifty]

No description provided.

[wsshin]

I wanted to test this PR, but the unit tests fail. For example, the following code taken from the unit tests generates an error:

julia> using ApproxFunBase, ContinuumArrays, OrthogonalPolynomialsQuasi

julia> f = Fun(Chebyshev(), [1.,2.,3.])
ERROR: MethodError: no method matching InfiniteArrays.Infinity(::Int64)
...

Could you make the unit tests work?

Also, most of the required packages described in Project.toml are outdated. It will be nice to bump up the versions of those packages and update the code if necessary. Once the unit tests start working again with the current Project.toml, they will reveal issues with the updated Project.toml, and interested contributors (including me) will be able to attempt to fix the compatibility issues.

[dlfivefifty]

Ah yes it's really old... might work if OrthogonalPolynomialsQuasi -> ClassicalOrthogonalPolynomials

[dlfivefifty]

@jishnub Could we have a Zoom call at some point to discuss this PR? I want to see how it fits into your plans.

I still need to decide what to do about operators. One plan is to make operators "linear functions of functions", which support algebraic manipulations and domains and ranges, but now there is no longer a notion of getindex as all of that is left to ContinuumArrays.jl. That is, we would do the following:

f = Fun(exp, ChebyshevT())
D = Derivative(ChebyshevT())
@test domain(D) == ChebyshevT() # Should this be more pedantic and be ColumnSpan(ChebyshevT())?
@test range(D) == ChebyshevU()

@test D(f) == f'
N = D*f; @test N(g) == g' * f
N = f*D; @test N(g) == f*g'
N = D ∘ D; @test N(f) == f'' # Note this actually makes sense since ColumnSpan(ChebyshevU()) ⊂ ColumnSpan(ChebyshevT())
@test D^2 == D ∘ D # debatable overloading

B = Dirichlet(Chebyshev())
@test B(f) == [f(-1),f(1)]

# Airy example
L = [B; D^2 - x*I]
@test L(f) == [f(-1); f(1); f'' - x*f]

u = L \ [airyai(-1); airyai(1); 0] # debatable overloading

Behind the scenes the last line would be lowered to

T = domain(L)
C = range(L)
D̃ = Diff(axes(T,1)) # Currently ContinuumArrays.Derivative but can change name
L̃ = [T[[begin,end],:]; C \ D̃^2 * T - C \ (x .* T)] # This an ∞ Matrix, which takes the place of `getindex`
c = L̃ \  [airyai(-1); airyai(1); zeros(∞)] # this is an ∞ vector
Fun(T, c)

But we may want to abandon operators completely to focus on auto-differentiation (which would also support nonlinear).

[jishnub]

Sorry for the delayed reply, I've been a bit preoccupied. I think I would have time to look into this in the second half of October and November. Perhaps we can schedule something for a weekend?

[dlfivefifty]

Sounds good, just propose a time. I'll try to get the conflicts fixed (though we might want to think of a different way of utilising ContinuumArrays.jl than what I did here)