I think Derivative should behave like UniformScaling and not require axes.

I discovered the following while reviving some untested functionality in CompactBases (see also JuliaApproximation/QuasiArrays.jl#32):

julia> using IntervalSets, QuasiArrays, ContinuumArrays

julia> r = Inclusion(0.0..1.0)
Inclusion(0.0..1.0)

julia> Q = QuasiDiagonal(r)
QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}

julia> Q*Q
ERROR: MethodError: no method matching Int64(::Infinities.InfiniteCardinal{1})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:763
  (::Type{T})(::AbstractChar) where T<:Union{Int32, Int64} at char.jl:51
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] _similar_for(c::UnitRange{Int64}, #unused#::Type{Float64}, itr::Inclusion{Float64, ClosedInterval{Float64}}, #unused#::Base.HasLength)
    @ Base ./array.jl:576
  [2] _collect(cont::UnitRange{Int64}, itr::Inclusion{Float64, ClosedInterval{Float64}}, #unused#::Base.HasEltype, isz::Base.HasLength)
    @ Base ./array.jl:609
  [3] collect(itr::Inclusion{Float64, ClosedInterval{Float64}})
    @ Base ./array.jl:603
  [4] map
    @ ./tuple.jl:214 [inlined]
  [5] QuasiArray(a::ApplyQuasiMatrix{Float64, typeof(*), Tuple{QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}}})
    @ QuasiArrays ~/.julia/dev/QuasiArrays/src/quasiarray.jl:84
  [6] QuasiArray(M::ArrayLayouts.Mul{QuasiArrays.QuasiArrayLayout, QuasiArrays.QuasiArrayLayout, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}})
    @ QuasiArrays ~/.julia/dev/QuasiArrays/src/lazyquasiarrays.jl:214
  [7] _quasi_mul(M::ArrayLayouts.Mul{QuasiArrays.QuasiArrayLayout, QuasiArrays.QuasiArrayLayout, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}}, #unused#::Tuple{Inclusion{Float64, ClosedInterval{Float64}}, Inclusion{Float64, ClosedInterval{Float64}}})
    @ QuasiArrays ~/.julia/dev/QuasiArrays/src/matmul.jl:94
  [8] copy(M::ArrayLayouts.Mul{QuasiArrays.QuasiArrayLayout, QuasiArrays.QuasiArrayLayout, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}})
    @ QuasiArrays ~/.julia/dev/QuasiArrays/src/matmul.jl:96
  [9] materialize(M::ArrayLayouts.Mul{QuasiArrays.QuasiArrayLayout, QuasiArrays.QuasiArrayLayout, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}, QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/wBWqq/src/mul.jl:105
 [10] mul
    @ ~/.julia/packages/ArrayLayouts/wBWqq/src/mul.jl:106 [inlined]
 [11] *(A::QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}}, B::QuasiDiagonal{Float64, Inclusion{Float64, ClosedInterval{Float64}}})
    @ QuasiArrays ~/.julia/dev/QuasiArrays/src/matmul.jl:23
 [12] top-level scope
    @ REPL[10]:1

This is on QuasiArrays 0.5.1 and dl/QuasiArrays-v0.5 of ContinuumArrays.

@jagot any opinions on renaming Derivative as Diff? This would be to be consistent with the use of diff for differentiating a quasi-vector.

It's best to have a pull request template for the better understanding of the changes made to the project!

At the moment factorize plays the role of both a transform (a Plan) and a coordinate transformation. But sometimes its convenient to work directly with transforms. For example, at the moment I'm trying to do the analogue of fft(randn(3,3,3), 1) for a Chebyshev() basis but factorizations are inherently matrix/vec.

My plan is to add an extra function plan_transform(basis, array, dim) to support this functionality.

It would be nice if one could do:

x = Inclusion(0..1)
sum(exp.(x))

and have sum automatically expand in a basis (in this case Chebyshev or Legendre). Same with differentiation, etc.

This would need to build on the singularity detection used in ClassicalOrthogonalPolynomials.jl.

So I can define a function f(x)=x^2, x∈[0,1] by

julia> f=Inclusion(0..1).^2
(Inclusion(0..1)) .^ 2

Then I want to do composition, say, g(x)=f(f(x)). But currently I'm not allowed to do:

julia> g=f[f]
ERROR: MethodError: no method matching Int64(::Infinities.InfiniteCardinal{1})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  (::Type{T})(::BigInt) where T<:Union{Int128, Int16, Int32, Int64, Int8} at gmp.jl:356
  ...
Stacktrace:
  [1] union!(s::Set{Float64}, itr::Inclusion{Float64, IntervalSets.ClosedInterval{Int64}})
    @ Base .\abstractset.jl:93
  [2] Set{Float64}(itr::Inclusion{Float64, IntervalSets.ClosedInterval{Int64}})
    @ Base .\set.jl:10
  [3] _Set(itr::Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, #unused#::Base.HasEltype)
    @ Base .\set.jl:23
  [4] Set(itr::Inclusion{Float64, IntervalSets.ClosedInterval{Int64}})
    @ Base .\set.jl:21
  [5] issubset(l::QuasiArrays.BroadcastQuasiVector{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, Base.RefValue{Val{2}}}}, r::Inclusion{Float64, IntervalSets.ClosedInterval{Int64}})
    @ Base .\abstractset.jl:279
  [6] checkindex(#unused#::Type{Bool}, inds::Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, i::QuasiArrays.BroadcastQuasiVector{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, Base.RefValue{Val{2}}}})
    @ QuasiArrays C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\indices.jl:203
  [7] checkbounds_indices
    @ .\abstractarray.jl:642 [inlined]
  [8] checkbounds
    @ C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\abstractquasiarray.jl:275 [inlined]
  [9] checkbounds
    @ C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\abstractquasiarray.jl:287 [inlined]
 [10] view
    @ C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\subquasiarray.jl:68 [inlined]
 [11] layout_getindex
    @ C:\Users\pty\.julia\packages\ArrayLayouts\3lnt1\src\ArrayLayouts.jl:109 [inlined]
 [12] _getindex(#unused#::Type, #unused#::IndexCartesian, A::QuasiArrays.BroadcastQuasiVector{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, Base.RefValue{Val{2}}}}, I::Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, Base.RefValue{Val{2}}}}})
    @ QuasiArrays C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\abstractquasiarray.jl:384
 [13] _getindex
    @ C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\abstractquasiarray.jl:375 [inlined]
 [14] getindex(A::QuasiArrays.BroadcastQuasiVector{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, Base.RefValue{Val{2}}}}, I::QuasiArrays.BroadcastQuasiVector{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, Inclusion{Float64, IntervalSets.ClosedInterval{Int64}}, Base.RefValue{Val{2}}}})
    @ QuasiArrays C:\Users\pty\.julia\packages\QuasiArrays\bxqu9\src\abstractquasiarray.jl:370
 [15] top-level scope
    @ none:1

Can we have this feature?

It's great to have an issue template for better understanding of the issues created for bugs, feature requests for the project!

For making operators like the Airy operator it's helpful to do things like:

D^2 - diagonal(x)

Should we export diagonal to make this easy? It's a risk for conflicting exports, though I don't know of any in the moment.

We could also support

D^2 - x .* I

even though that's not supported in Base.

At the moment we do expansions as

T * (T \ f)

but it might be convenient to allow an expansion operator. This would obviously be T/T, which is short-hand for T*inv(T). So the above could be changed to

(T/T)\f

I'm surprised it doesn't work already but probably just need to make sure things are flagged as lazy:

julia> Z / Z
ERROR: MethodError: no method matching similar(::Type{Array{Float64, N} where N}, ::Tuple{QuasiArrays.Inclusion{SVector{2, Float64}, DomainSets.FixedUnitBall{SVector{2, Float64}, :closed}}, QuasiArrays.Inclusion{SVector{2, Float64}, DomainSets.FixedUnitBall{SVector{2, Float64}, :closed}}})
Closest candidates are:
  similar(::QuasiArrays.AbstractQuasiArray{T, N} where N, ::Tuple) where T at /Users/sheehanolver/.julia/packages/QuasiArrays/Eszna/src/abstractquasiarray.jl:295
  similar(::AbstractBlockArray{T, N} where N, ::Tuple) where T at /Users/sheehanolver/.julia/packages/BlockArrays/g01fJ/src/abstractblockarray.jl:38
  similar(::AbstractArray{T, N} where N, ::Tuple) where T at abstractarray.jl:734
  ...
Stacktrace:
 [1] similar(A::ArrayLayouts.Rdiv{ContinuumArrays.BasisLayout, ContinuumArrays.BasisLayout, Zernike{Float64}, Zernike{Float64}}, #unused#::Type{Float64}, axes::Tuple{QuasiArrays.Inclusion{SVector{2, Float64}, DomainSets.FixedUnitBall{SVector{2, Float64}, :closed}}, QuasiArrays.Inclusion{SVector{2, Float64}, DomainSets.FixedUnitBall{SVector{2, Float64}, :closed}}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/75WvP/src/ldiv.jl:17
 [2] similar(A::ArrayLayouts.Rdiv{ContinuumArrays.BasisLayout, ContinuumArrays.BasisLayout, Zernike{Float64}, Zernike{Float64}}, #unused#::Type{Float64})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/75WvP/src/ldiv.jl:18
 [3] similar(A::ArrayLayouts.Rdiv{ContinuumArrays.BasisLayout, ContinuumArrays.BasisLayout, Zernike{Float64}, Zernike{Float64}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/75WvP/src/ldiv.jl:19
 [4] copy
   @ ~/.julia/packages/ArrayLayouts/75WvP/src/ldiv.jl:21 [inlined]
 [5] materialize
   @ ~/.julia/packages/ArrayLayouts/75WvP/src/ldiv.jl:22 [inlined]
 [6] rdiv
   @ ~/.julia/packages/ArrayLayouts/75WvP/src/ldiv.jl:87 [inlined]
 [7] /(A::Zernike{Float64}, B::Zernike{Float64})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/Eszna/src/matmul.jl:38
 [8] top-level scope
   @ REPL[32]:1

I recently added jacobimatrix(P,n) to ClassicalOrthogonalPolynomials.jl as it was sort of needed for n x n truncations that preserved the Tridiagonal structure. But perhaps this interface makes sense for other operators like weaklaplacian and grammatrix?

I discovered the following, not-critical bug:

using ContinuumArrays

S = LinearSpline{Float64}(0.0:3)
x = range(0,stop=3,length=11)
# Evaluate basis at finer grid
χ = S[x,:]*S'

u = S*rand(size(S,2))

uv = χ*u
display(uv) # Displays vector of #undef's
uv[1] # Retrieves correct value

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

What do people think of

diff(f, order=3) # 3rd derivative
diff(f, order=(1,1)) # partial derivative f_{xy}

Behind the scenes it will wrap in a Val and because of constant propagation this should maintain type-stability..

Background:
I want to find all zero crossings of a function, expanded over a basis. I decided I'd give IntervalRootFinding.jl a shot. In essence, IntervalRootFinding.roots needs to evaluate the function (and its derivative) on an interval x and expects an interval y back which holds the minimum and maximum values attained by the function (derivative) on x.

Now, since my function is expressed as a linear combination of basis functions, my best idea so far is to find the corresponding interval y_i for the ith basis function, multiply this interval by the expansion coefficient (which will flip the interval if the coefficient is negative) and sum all the minima and maxima. We cannot simply use the union of the intervals since two functions on top each other would necessarily increase the interval. At the same time, for two disjoint functions, adding the intervals would make the bound looser than it necessarily needs to be. Is there a way to find a tighter bound?

How would the interface for finding the extrema of each basis function over an interval look like? If B isa Basis, with continuous first dimension, extrema(B, interval) could e.g. return a vector of intervals, but I am not sure if this is the best way. It would also be dependent on which particular implementation of Interval one chooses (IntervalSets.jl, IntervalArithmetic.jl, etc.).

Finally, is there a more clever way to find roots than this?

I want to define the p-norm for MulQuasiArrays for my finite differences like so:

const FDVector{T,B<:AbstractFiniteDifferences} = MulQuasiArray{T,2,<:Mul{<:Tuple,<:Tuple{<:B,<:AbstractVector}}}
const FDMatrix{T,B<:AbstractFiniteDifferences} = MulQuasiArray{T,2,<:Mul{<:Tuple,<:Tuple{<:B,<:AbstractMatrix}}}
const FDVecOrMat{T,B} = Union{FDVector{T,B},FDMatrix{T,B}}

function LinearAlgebra.norm(v::FDVecOrMat, p::Real=2)
    B,c = v.mul.factors
    norm(c, p)*(step(B)^(inv(p)))
end

function LinearAlgebra.normalize!(v::FDVecOrMat, p::Real=2)
    v.mul.factors[2][:] /= norm(v,p)
    v
end

norm works as intended, but normalize! leaves v unchanged (but if I debug-print the norm in normalize! it works). What am I doing wrong? Is it a scoping/ownership issue?

A common pattern to expand a function is:

P = Legendre()
x = axes(P,1)
f = P / P \ exp.(x)

But perhaps a shortcut would be useful so this is one line:

expand(Legendre(), exp)

That is, add the function:

function expand(P, f)
    x = axes(P,1)
    P / P \ f.(x)
end

At the moment we have the following setup for linear operators / functionals that are not built out of elementary operations, where f isa AbstractQuasiVector and T isa AbstractQuasiMatrix:

1. Derivative(x) * f and Derivative(x) * T
2. sum(f) and sum(T; dims=1)
   As we consider other linear operators such as cumsum it might be a good idea to see how to make this consistent

Version (2) is in some sense "canonical" as it exists in Base. So the question is what to do with Version (1). We could have something like

struct Linear{T,Axes,F} <: AbstractQuasiMatrix{T}
    f::F
    axes::Axes
end

*(L::Linear, v::AbstractQuasiVector) = L.f(v)
*(L::Linear, v::AbstractQuasiMatrix) = L.f(v; dims=1)
axes(L::Linear) = L.axes

const Diff{T,Axis} = Linear{T, Axis,typeof(diff)} # Replaces Derivative
const Cumsum{T,Axis} = Linear{T, Axis,typeof(cumsum)}
const Sum{T,Axis} = Linear{T,Axis,typeof(sum)}

Diff(x) = Linear{Float64}(diff, (x,x))
Cumsum(x) = Linear{Float64}(cumsum, (x,x))
Sum(x) = Linear{Float64}(diff, (Base.OneTo(1),x))

How this relates to axes-free alternatives is another question: #22

Isn’t a FunctionSpace (almost by definition) just a Space with evaluation defined (that respects the linearity of the space)?

It almost feels like Fun is misnomer, as the space could be something other than a FunctionSpace. Perhaps its more general to do the following:

# In FunctionSpaces.jl
struct Vec{S<:Space, T, VT}
   space::S
   coefficients::VT
end

# In ApproxFun.jl
const Fun{S, T, VT} = Vec{S<:FunctionSpace, T, VT}

@jagot do you use the splines here, or should they be moved out and this made a "minimal" package?

EDIT: This was my original motivation, see next post for a more generalized take on the problem.

The Schrödinger equation in spherical coordinates

is singular at the origin. The solution has to obey the following boundary conditions:

In finite-differences, it is common to adjust the upper-left corner of the derivative matrices for \ell = 0 to implement this and preserve the convergence order (see refs below). I wonder how I could accomplish this with ContinuumArrays.jl, namely I wish to dispatch an expression like

R = ...
D = Derivative(axes(R,1))
Δ = R'D'D*R

correctly such that I get the correct boundary conditions? One idea could be to introduce

abstract type AbstractDerivative{T} <: LazyQuasiMatrix{T} end

struct Derivative{T,D} <: AbstractDerivative{T}
    axis::Inclusion{T,D}
end
...

# Then I could have my own type
struct SingularDerivative{T} <: AbstractDerivative{T}
    axis::Inclusion{T,D}
end

• If R is aware of SingularDerivative
  • If \ell == 0, apply correction
  • Else, return uncorrected matrix
• If R is not aware of SingularDerivative, proceed as normal

The problem with this approach is that suddenly all basis libraries will need to dispatch on AbstractDerivative.

A better alternative would of course be if I could somehow attach boundary conditions to the derivative operator:

D = Derivative(axes(R,1), :singular, :regular)

Are there better approaches?

References

• Schafer, K. J., Gaarde, M. B., Kulander, K. C., Sheehy, B., &
  DiMauro, L. F. (2000). Calculations of strong field multiphoton
  processes in alkali metal atoms. AIP Conference Proceedings, 525(1),
  45–58. http://dx.doi.org/10.1063/1.1291925

• Schafer, K. J. (2009). Numerical methods in strong field physics. In
  T. Brabec (Eds.), (pp. 111–145). : Springer.

• Muller, H. G. (1999). An Efficient Propagation Scheme for the
  Time-Dependent Schrödinger equation in the Velocity Gauge. Laser
  Physics, 9(1), 138–148.

• Patchkovskii, S., & Muller, H. (2016). Simple, accurate, and
  efficient implementation of 1-electron atomic time-dependent
  Schrödinger equation in spherical coordinates. Computer Physics
  Communications, 199(nil),
  153–169. http://dx.doi.org/10.1016/j.cpc.2015.10.014

We should add support for plotting using grids/transforms. Anyone have an opinion on whether to support:

👍 Makie
🎉 Plots
🚀 Both

Please vote below. @jagot @TSGut

Personally I think Makie will be really beneficial for complicated geometries like disks/triangles/etc. but perhaps Plots.jl can do this too...

I was introduced as Continous Linear Algebra by Nick Trefethen which demonstrated the issue of using monomials a basis using chebfun, so i expected similar functionality from ApproxFun which is apperently better addressed here.

I was failing to get a result for rank([ 1 sin(x)^2 cos(x)^2]) which would be 2 but there is a lot more that can be done with continous linear algebra.

What would be a sensible interface for computing integrals? If 'B' is a Basis (or suitable QuasiMatrix), c the expansion coefficients of a function, and f = B*c the function expanded over B, I'm thinking

• integral(f) would give me the integral of the function,
• integral(B) would give me an adjoint vector with the integral of each individual basis function, which I can dot with c to get the same as integral(f),
• integral(B[Inclusion(a..b),:]) would integrate the basis functions on a subinterval a..b,
• integral(B*QuasiDiagonal(f.(axes(B,1))) would integrate each basis function, weighted by the function f(x).

Instead of integral, I guess sum would be a good name, but that's already taken for something I am not completely sure what it does (maybe exactly this?).

In any case, it tells me I need to Override for ContinuumArrays.BasisLayout(). If this is what I want, what is the canonical way to do this?

Could we claim that two WeightedBasis are equal iff both the weights and the basis are equal?

==(A::WeightedBasis, B::WeightedBasis) = A.args == B.args

If so, we can omit this line in OrthogonalPolynomialsQuasi.jl/jacobi.jl:

==(A::WeightedJacobi, B::WeightedJacobi) = A.args == B.args

I had been thinking of operators in the kernel sense, that is,
$$
(K u)(x) = \int_a^b K(x,y) u(y) dy
$$
So that K[x,y] is equivalent to K(x,y).

But then we are using diagonal operators. One might expect QuasiDiagonal(one.(Inclusion(-1...1))) to be the identity operator, and thats essentially how we've been using it. The issue is that the identity operator corresponds to K(x,y) = δ(y-x) so is actual infinite on the diagonal.

Is this the convention we want? E.g., QuasiDiagonal(sqrt.(1-x.^2))[0.1,0.1] == ∞?

At the moment the following makes a lazy *, but I think it should be left as a view:

julia> S = LinearSpline{Float64}(0.0:3);

julia> S[:,2:end-1] # Lazy S * Project
QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Spline{1,Float64},BandedMatrices.BandedMatrix{Int64,FillArrays.Ones{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}}}(*, (Spline{1,Float64}([0.0, 1.0, 2.0, 3.0]), [0 0; 1 0; 0 1; 0 0]))

julia> view(S,:,2:size(S,2)-1) # proposal: returns this
QuasiArrays.SubQuasiArray{Float64,2,Spline{1,Float64},Tuple{Inclusion{Float64,IntervalSets.Interval{:closed,:closed,Float64}},UnitRange{Int64}},false}(Spline{1,Float64}([0.0, 1.0, 2.0, 3.0]), (Inclusion(0.0..3.0), 2:3), 0, 0)

I believe this is what you were asking for @jagot.

The steps to make it work are

• Support lazy multiplication with views
• Fix getindex for lazy multiplication of views
• Lower D*view(S,:,kr) to view(D*S,:,kr)

I will add to the below list, as I think of more things.

Derivatives

D = Derivative(axes(B,1))

Since ∂' = -∂ (anti-symmetric), I think it's worth to point out that a weak Laplacian is -B'*D'*D*B. Furthermore, I assume that B'*(D*D*B) is a strong (?) Laplacian, i.e. the right basis is differentiated twice and is then projected on the left basis.

Higher derivatives

This is not very common in my applications, but I guess a D^n operator would be nice, instead of D*D*D*D...*D*D. For weak vs strong derivatives, I guess the user would write B'*(D^a)'*D^b*B.

Diagonal operators

This expresses the kind of operator that only locally depends on the coordinate, e.g. a potential, that when acting on a function, yields another function that's multiplied by the potential: V*f -> g(x) = V(x)*f(x). This is different from function interpolation below, although superficially similar in orthogonal bases (e.g. finite-differences); in e.g. B-splines, such operators become banded matrices. At the moment, my syntax for this is V = Matrix(x -> -inv(x), B) for e.g. the Coulomb potential, or V = Matrix(x -> -inv(x), B_1, B_2) for differing left and right bases (see below).

Function interpolation

Right now, to expand a function on a basis, I have defined B \ f::Function, which evaluates the function f on the quadrature points (or equivalently) and then does the least-square approximation with respect to the basis functions on the same points.

Properties of bases ("Internals")

Some of the below properties are not applicable to all bases, particularly not global bases, I guess, but are nevertheless handy to have for investigation and logging purposes (Which order did I use for this calculation?)

• order polynomial order of basis functions, defined as polynomial degree + 1; usually decides convergence rate of spatial derivatives/integrals. E.g. finite-differences with the [-1 2 -1] stencil for the 2nd derivative, can be viewed as 2nd order piecewise polynomials, i.e. tent functions.
• locs & weights locations and weights of quadrature points used in e.g. B-splines and finite-elements.

Variable transformations

This is something we discussed some time ago; in e.g. atomic physics, it's common to work with log(r) instead of r to increase accuracy of wavefunctions. However, one must then work with another set of equations. It would be cool to stick a variable transformation "between" the equations and the basis, such that neither the user writing the equations, nor the basis implementing basis functions and derivatives etc, need to care about the distribution of grid points. You mentioned that auto-differentiation could possibly be used for this.

Multiple bases

Sometimes it's useful to express two different functions in separate bases. If these functions are coupled via an equation system, the off-diagonal terms will need to go from one basis to another. Right now, I'm working on this for B-splines, where I use two bases of order k and k-1 to avoid spurious eigenstates of the Dirac Hamiltonian. I can do this easily, since I use the same knot set and the same set of quadrature points/weights, but for two completely different bases, e.g. B-splines and finite-differences, how would one accomplish e.g. B_1'*D*B_2, without forcing each basis implementation to know about all other possible bases?

How do we represent change of variables? E.g., P_n(2x-1) for Legendre on 0..1.

I think the "right" way to do this is to treat change of variables also as quasi-matrices. That is the map u -> (x -> u(2x-1)) is a linear operator from functions on -1..1 to functions on 0..1. Formally this can be thought of as convolution with a delta function:
$$
\int_{-1}^1 \delta(y - (2x-1)) u(y) dy
$$
This hints that we can represent such an operator as:

x, y = Inclusion(0..1), Inclusion(-1..1)
R = δ.(y' .- (2.*x .- 1))
R * Legendre() 

How does this interact with other operators like derivative? I think quite nicely: we would just need a simplifying relationship so that:

D = Derivative(x)
D*R == 2R*D

Then full materialize would automatically reduce:

D*R*Legendre() == 2R*D*Legendre() == 2R*Jacobi(1,1)*B

where B is some banded matrix.

PS: this assumes the current DiracDelta <: AbstractQuasiVector being replaced with a δ that behaves like a function, and then δ.(Inclusion(d) .- x) replacing DiracDelta(d,x).

PPS: We would need to make sure constant multiplication is "passed through" to the matrix.

PPPS: The simpler approach is to do what ApproxFun did and encode the domain inside Legendre. This has the effect that every space needs to know about coordinate transformations. This worked but I believe the R approach above to be more versatile (it completely decouples coordinate transformations from bases), and so crazy that I really want to to do it 🤣

Assume a linear operator Ĥ with a discrete and continuous
spectrum. Then any function may be expanded over its eigenfunctions:

When calculating inner products and matrix elements between such
expansions, one encounters inner products of the kinds

Usually, the eigenfunctions of the discrete spectrum have compact
support, or are at least very localized (they may have an
exponentially decaying tail, which is negligible). It would be nice to
be able to calculate such inner products, without having to represent
the discrete eigenfunctions on the entire interval. Maybe something
like this?:

B <: AbstractQuasiMatrix
B̃ = B[1:100]
i = B̃*rand(100)
k = B*rand(size(B,2))

# This is the syntax I'd like to automagically work:
i'k

As a follow-up to this, I'll need to solve the Poisson problem on the
entire interval, with source terms of the form

where the source term will be non-zero only on the interval where
is.

Say A isa Basis and B isa Basis, where axes(A,1) != axes(B,1), then we get the following error:

julia> A'B
ERROR: DimensionMismatch: Second axis of A, Inclusion(0.0..2.0), and first axis of B, Inclusion(0.0..2.5) must match

However, if I have a function $f$ that is expanded over the basis A, $|f\rangle = \sum_i c_i|a_i\rangle$, and I want to project it onto the basis B, $|f\rangle = \sum_j d_i|b_j\rangle$, then the above basis overlap is precisely the one I need, since $d_j = \sum_i c_i \langle b_j|a_i\rangle$, or equivalently $\mathbf{d}=A^HB\mathbf{c}$.

To me it is unimportant that the axes are different, i.e. I think it is up to the user to make sure that the end result makes sense (e.g. projecting onto a "larger" basis, or maybe I am interested in the subspace, etc).

I can of course override this on a case-by-case basis, but it would be convenient if this just worked, and the only thing you'd need to provide is the calculation of basis function overlaps $\langle b_j|a_i\rangle$. Then it would be very easy to transform between e.g. B-splines and finite-elements or -differences.

What do you think?

EDIT: It seems it will be difficult to work around this on a case-by-case basis, since I cannot circumvent the dimension check easily:

julia> ApplyQuasiArray(*, A', B)
ERROR: DimensionMismatch: Second axis of A, Inclusion(0.0..2.0), and first axis of B, Inclusion(0.0..2.5) must match
Stacktrace:
 [1] _check_mul_axes(A::QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, B::QuasiArrays.SubQuasiArray{Float64, 2, FEDVR{Float64, Float64, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}, Tuple{Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, UnitRange{Int64}}, false})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:91
 [2] check_mul_axes
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:107 [inlined]
 [3] check_applied_axes
   @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:22 [inlined]
 [4] instantiate
   @ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:63 [inlined]
 [5] ApplyQuasiArray
   @ ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:53 [inlined]
 [6] ApplyQuasiArray
   @ ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:54 [inlined]
 [7] ApplyQuasiArray(M::LazyArrays.Applied{LazyArrays.MulStyle, typeof(*), Tuple{QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, QuasiArrays.SubQuasiArray{Float64, 2, FEDVR{Float64, Float64, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}, Tuple{Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, UnitRange{Int64}}, false}}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:55
 [8] ApplyQuasiArray(::Function, ::QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, ::Vararg{Any})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:59

I assume that's what's wrong with this

julia> using QuasiArrays, ContinuumArrays

julia> import ContinuumArrays: @simplify

julia> struct MyOp{T,Ax} <: QuasiArrays.LazyQuasiMatrix{T}
           axis::Ax
       end

julia> @simplify function *(A::MyOp, P::Basis{T}) where T
           P * Diagonal(1:∞)
       end
ERROR: LoadError: AssertionError: (qt.args[1]).head == :call
Stacktrace:
 [1] var"@simplify"(__source__::LineNumberNode, __module__::Module, qt::Any)
   @ ContinuumArrays ~/.julia/packages/ContinuumArrays/tg4mj/src/operators.jl:21
in expression starting at REPL[5]:1

Maybe it's not an issue because macros work with types?

I needed to add copy and copyto! methods for MulQuasiArrays, and nothing worked until I realized that they were not imported from Base in QuasiArrays.jl, effectively making the copy methods in abstractquasiarray.jl hidden. It's easy enough to add more imports to QuasiArrays.jl, but recently I've started explicitly prepending the name of the module to the method I'm overloading, to decrease confusion, e.g. Base.copyto!(...) =. What do you think about this?

Related to #8

I started thinking again how to implement overlapping domains; when solving PDEs like the Schrödinger equation, different but equivalent formulations are more/less stiff in different parts of space. The formulations are related via a unitary transformation. One can then solve two problems in slightly overlapping domains and communicating the overlap between them every time step.

Setting up the equations is easy, as soon as the correct AbstractQuasiMatrices are known, the question is rather how to split the whole domain into subdomains. Imagine

R = FiniteDifferences(0.0..100, δx=1.0)

I now want to split R into subdomains R1 and R2, covering 0.0..30 and 20.0..100, respectively. I also want all materialized matrices to be of appropriate sizes, i.e. not full matrices with zeros outside their respective subdomains. This can be accomplished by

R1 = R[:,1:31]
R2 = R[:,21:end]

but then the first axes of R1 and R2 will be "wrong" and if R is not a uniform grid, as in this case, it is not so trivial to determine which basis functions need to be included. If on the other hand I do

R1 = R[Inclusion(0.0..30),:]
R2 = R[Inclusion(20.0..100),:]

then the first axes are correct but all the basis functions are still there, so the materialized matrices will still be too large.

• What would be a good interface for this?
• This obviously only makes sense for local bases, e.g. something like orthogonal polynomials would not fit into this. It's also not clear how this would translate to B-splines that have compact support, but are not orthogonal, i.e. the overlap matrix is banded.

Want to do:

T = ChebyshevT()
x = axes(T,1)
K = (x,y) -> exp(x*cos(y)) # some kernel
f = expand(T, exp) # some function
K.(x, x') * f # this is integration wrt a kernel

We have some experiments for 2D quasivectors. Here K.(x, x') is a quasi-matrix. Do we want to support bases as 4-tensors: i.e. T²[x,y,k,j] so that if X isa AbstractMatrix we can write K = T² * X ?

Functions in this basis are represented by a lazy multiplication by a basis and a vector of coefficients:

julia> f = L*[1,2,3,4,5,6]
QuasiArrays.ApplyQuasiArray{Float64,1,typeof(*),Tuple{Spline{1,Float64},Array{Int64,1}}}(*, (Spline{1,Float64}([0.0, 0.2, 0.4, 0.6, 0.8, 1.0]), [1, 2, 3, 4, 5, 6]))

julia> axes(f)
(Inclusion(0.0..1.0),)

julia> f[0.1]
1.5

However, if it was f(0.1) at the end, not only would it look more like a function, but f.([0.1,0.2]) would then work.

• h-FEM. TetGen.jl might provide a good way of creating meshes
• p-FEM. In 1D, these can be built from (1-x^2)*P_n^(1,1)(x) combined with splines. For triangles, we need to port operators from MultivariateOrthogonalPolynomials.jl
• Spherical harmonics.
• Ultraspherical. This links in with InfiniteArrays.jl

This is a follow-up to #14, but more speculatory.

At the moment, I'm working with atoms with spherical symmetry, which allows me to reduce my differential equations to one-dimensional radial problems, with the angular coordinates treated exactly using angular momentum algebra. In the future, I will probably want/need to go to molecules which possess less/other symmetries. There are many different ways of treating molecules, but the dominating by far is linear combinations of (multi-centred) 3d Gaussians (actually polynomials multiplied by Gaussians). It would be interesting to see if we could fit these type of basis functions within the ContinuumArrays framework. If we get this to work, a Hartree–Fock code could conceivably be used for both atoms and molecules without modifications (some tensor magic would need to happen).

I am by no means well-versed on the topic, and would need to read up quite a bit (although I have a lot of good references). @jarvist could possibly be interested in this as well.

At the moment new operators are defined via @simplify. This was probably me being too clever for my own good and is very hard on the compiler.

With a bit of thought we can redesign it so new operators can be added and used without having to mess with @simplify. For example, rather than overloading @simplify *(::Derivative, ::MyBasis) we can instead overload _diff(::MyBasis, dims) and then a user can access this via diff(MyBasis(); dims=1) without any appeal to the simplify machinery.

Note this probably needs more thought: one normally does things like

T,U = ChebyshevT(), ChebyshevU()
U \ diff(T; dims=1)

But the second command still goes through the simplify machinery. But maybe this is fine.

Sometimes, I need to orthogonalize a MulQuasiVector against a set of
other MulQuasiVectors:

Currently, I do it like this:

const RadialOrbital{T,B<:AbstractQuasiMatrix} = MulQuasiArray{T,1,<:Mul{<:Tuple,<:Tuple{<:B,<:AbstractVector}}}

struct Projector{ΦT,B<:AbstractQuasiMatrix,RO<:RadialOrbital{ΦT,B}}
    ϕs::Vector{RO}
end

projectout!(y::RO, ::Nothing) where RO = y

function projectout!(y::RO, projector::Proj) where {RO,Proj<:Projector}
    yc = y.mul.factors[2]

    for ϕ in projector.ϕs
        c = ϕ'y
        yc .-= c*ϕ.mul.factors[2]
        # y -= c*ϕ # It would be nice if this worked
    end
    y
end

Could we get y -= c*ϕ syntax working? It would be nice if that did not have to be implemented for every kind of basis, but just in ContinuumArrays.

This (and more complicated examples) have worked previously, but no longer:

julia> using IntervalSets, QuasiArrays, ContinuumArrays

julia> r = Inclusion(0.0..50.0)
Inclusion(0.0..50.0)

julia> v = 40.0 .< r .≤ 50.0
BroadcastQuasiArray{Bool,1,typeof(&),Tuple{BroadcastQuasiArray{Bool,1,typeof(<),Tuple{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}},Base.Broadcast.Broadcasted{QuasiArrays.LazyQuasiArrayStyle{1},Nothing,typeof(<=),Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},Float64}}}}(&, (BroadcastQuasiArray{Bool,1,typeof(<),Tuple{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}(<, (40.0, Inclusion(0.0..50.0))), Base.Broadcast.Broadcasted(<=, (Inclusion(0.0..50.0), 50.0))))

julia> v[1.0]
ERROR: MethodError: no method matching error_if_canonical_getindex(::IndexLinear, ::BroadcastQuasiArray{Bool,1,typeof(<),Tuple{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}, ::Tuple{QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndex{1,Tuple{Float64}}})
Closest candidates are:
  error_if_canonical_getindex(::IndexStyle, ::AbstractArray, ::Any...) at abstractarray.jl:992
  error_if_canonical_getindex(::IndexCartesian, ::AbstractQuasiArray{T,N}, ::Tuple) where {T, N} at /Users/jagot/.julia/packages/QuasiArrays/3FYr1/src/abstractquasiarray.jl:400
  error_if_canonical_getindex(::IndexLinear, ::AbstractArray, ::Int64) at abstractarray.jl:988
Stacktrace:
 [1] _getindex at /Users/jagot/.julia/packages/QuasiArrays/3FYr1/src/abstractquasiarray.jl:378 [inlined]
 [2] getindex(::BroadcastQuasiArray{Bool,1,typeof(<),Tuple{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndex{1,Tuple{Float64}}) at /Users/jagot/.julia/packages/QuasiArrays/3FYr1/src/abstractquasiarray.jl:374
 [3] _broadcast_getindex at ./broadcast.jl:584 [inlined]
 [4] _getindex at ./broadcast.jl:627 [inlined]
 [5] _broadcast_getindex at ./broadcast.jl:603 [inlined]
 [6] getindex at /Users/jagot/.julia/packages/QuasiArrays/3FYr1/src/quasibroadcast.jl:116 [inlined]
 [7] getindex at /Users/jagot/.julia/packages/QuasiArrays/3FYr1/src/quasibroadcast.jl:119 [inlined]
 [8] getindex(::BroadcastQuasiArray{Bool,1,typeof(&),Tuple{BroadcastQuasiArray{Bool,1,typeof(<),Tuple{Float64,Inclusion{Float64,Interval{:closed,:closed,Float64}}}},Base.Broadcast.Broadcasted{QuasiArrays.LazyQuasiArrayStyle{1},Nothing,typeof(<=),Tuple{Inclusion{Float64,Interval{:closed,:closed,Float64}},Float64}}}}, ::Float64) at /Users/jagot/.julia/packages/QuasiArrays/3FYr1/src/lazyquasiarrays.jl:131
 [9] top-level scope at REPL[13]:1

  [7ae1f121] ContinuumArrays v0.2.4
  [8197267c] IntervalSets v0.5.1
  [c4ea9172] QuasiArrays v0.2.1

I will try to produce the more complicated example, that I need to have working.