JuliaApproximation / QuasiArrays.jl

A package for representing quasi-arrays
MIT License
12 stars 2 forks source link

Update for New ArrayLayouts #38

Closed dlfivefifty closed 4 years ago

dlfivefifty commented 4 years ago

@jagot this is in queue for being merged and completely revamps the simplification setup (which has been moved to https://github.com/JuliaArrays/LazyArrays.jl/pull/121)

This should make @simplify type stable, which should improve compile time and make it more robust.

jagot commented 4 years ago

Very nice! I did not follow very closely, but is the interface the same? I guess I will notice when all CompactBases.jl tests break :)

dlfivefifty commented 4 years ago

Hmm, for CompactBases.jl it should be , but let me a fork and do some testing

dlfivefifty commented 4 years ago

Tried it but realised I don't understand @materialize. I think we'll just have to work on updating together

codecov[bot] commented 4 years ago

Codecov Report

Merging #38 into master will increase coverage by 0.29%. The diff coverage is 69.01%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #38      +/-   ##
==========================================
+ Coverage   39.98%   40.28%   +0.29%     
==========================================
  Files          17       17              
  Lines        1258     1204      -54     
==========================================
- Hits          503      485      -18     
+ Misses        755      719      -36     
Impacted Files Coverage Δ
src/QuasiArrays.jl 70.00% <ø> (+26.25%) :arrow_up:
src/abstractquasiarray.jl 32.23% <ø> (-0.25%) :arrow_down:
src/quasiadjtrans.jl 43.67% <0.00%> (ø)
src/inv.jl 33.33% <33.33%> (-3.34%) :arrow_down:
src/quasiarray.jl 84.61% <60.00%> (-9.33%) :arrow_down:
src/matmul.jl 48.57% <71.73%> (+5.84%) :arrow_up:
src/lazyquasiarrays.jl 69.73% <90.00%> (+1.44%) :arrow_up:
src/indices.jl 41.60% <0.00%> (-0.92%) :arrow_down:
... and 2 more

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 d7e95cc...502b019. Read the comment docs.

jagot commented 4 years ago

Does this PR mean that we get B*A*inv(B) where B<:Basis to just work when applying to B*c, i.e. B*A*c? If so, 🎉 , but how do we handle non-orthogonal bases, where A can take you out of B (as I try to describe here: https://juliaapproximation.github.io/CompactBases.jl/stable/bsplines/theory/#Solving-equations-1)?

In CompactBases.jl I added the stop-gap measure operator_metric and LinearOperator to automatically apply the inverse metric in those cases it's necessary (mainly B-splines), but it maybe feels a bit too magical to do so when computing B*A*inv(B)*B*c.

Re update of CompactBases.jl, your branch passes CI, so I guess that's good?

dlfivefifty commented 4 years ago

Should do but we will likely need to iterate a bit. I don't see the problem with non-orthogonal bases, these are used all the time in OrthogonalPolynomialQuasi.jl.

jagot commented 4 years ago

It's actually not so much about non-orthogonal spaces as it is about which space your operator takes you to, i.e. differentiating a cubic spline formally gives you a quadratic spline, but since in computations you usually want to stay within one space, you project onto the basis again, which then means you have to solve by the mass matrix/metric.

For finite-differences and finite-elements, this is not an issue, since the operators are defined to take you back to the same space.

One the one hand, I expect B*A*inv(B)*B*c to be in the correct basis (which would imply the solution step), on the other hand you may want to avoid the solution step for computational efficiency, if you're going to compute an inner product directly afterwards, i.e. <x | A | y> ~ x*S*(S\A*y) == x*A*y (whereas an finite-differences it would be x*S*A*y). At the moment, I've hacked around this with matrix_element_metric, which I don't think is so beautiful.

Re @materialize, I have to remember what I did, but it's mostly a souped-up version of @simplify which supports multiple arguments.

dlfivefifty commented 4 years ago

yes that's exactly the point of ultraspherical spectral method:

D T = U (...)