JuliaApproximation / HarmonicOrthogonalPolynomials.jl

A Julia package for working with spherical harmonic expansions
MIT License
24 stars 2 forks source link

Controlling tolerance during function approximation #53

Open wsshin opened 2 years ago

wsshin commented 2 years ago

It will be nice to have a means to control the tolerance during function approximation. In other words, using the example code in README, I would like to have something like

julia> S = SphericalHarmonic();

julia> 𝐱 = axes(S,1);

julia> f = 𝐱 -> ((x,y,z) = 𝐱; exp(x)*cos(y*sin(z)));

julia> f̃ = S * ( \(S, f.(𝐱); tolerance=1e-6) ); # expansion of f in spherical harmonics with tolerance

where the last line specifies the tolerance during \, which is motivated by this section of the ApproxFun.jl documentation; see the last line of the section.

Not sure in which package this feature needs to be implemented: this package, ContinuumArrays.jl, or ClassicalOrthogonalPolynomials.jl? If basic directions are provided, I don't mind generating a PR myself, so please advise.

dlfivefifty commented 2 years ago

At the moment support for infinite axes is in ClassicalOrthogonalPolynomials.jl, so the transform is here:

https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/blob/main/src/adaptivetransform.jl

The design of how to do this is not clear as it might be better to allow users to specify their own convergence tests. There is also a question of relative and absolute tolerance.

Though I guess for now passing a keyword as you suggest is fine. Unfortunately because the code is designed on "traits" there will be a lot of functions across different packages. Below is an abridged debugger walkthrough of what needs to be modified, e.g. by changing the first one to:

 \(A::AbstractQuasiArray, B::AbstractQuasiArray; kwds...) = ldiv(A, B; kwds...)

Here is the walkthrough:

In \(A, B) at /Users/solver/.julia/packages/QuasiArrays/1I3a9/src/matmul.jl:34

>34  @inline \(A::AbstractQuasiArray, B::AbstractQuasiArray) = ldiv(A,B)

In ldiv(A, B) at /Users/solver/.julia/packages/ArrayLayouts/tg77x/src/ldiv.jl:86
>86  @inline ldiv(A, B) = materialize(Ldiv(A,B))

In materialize(M) at /Users/solver/.julia/packages/ArrayLayouts/tg77x/src/ldiv.jl:22
>22  @inline materialize(M::$Typ) = copy(instantiate(M))

In copy(L) at /Users/solver/.julia/packages/ContinuumArrays/5vkIF/src/bases/bases.jl:271
>271  copy(L::Ldiv{<:AbstractBasisLayout,<:BroadcastLayout}) = transform_ldiv_if_columns(L)

In transform_ldiv_if_columns(L) at /Users/solver/.julia/packages/ContinuumArrays/5vkIF/src/bases/bases.jl:270
>270  transform_ldiv_if_columns(L) = transform_ldiv_if_columns(L, axes(L.B,2))

In transform_ldiv_if_columns(L, #unused#) at /Users/solver/.julia/packages/ContinuumArrays/5vkIF/src/bases/bases.jl:269
>269  transform_ldiv_if_columns(L, ::OneTo) = transform_ldiv(L.A,L.B)

In transform_ldiv(A, B) at /Users/solver/.julia/packages/ContinuumArrays/5vkIF/src/bases/bases.jl:261
>261  transform_ldiv(A, B) = transform_ldiv(A, B, size(A))

In transform_ldiv(A, f, #unused#) at /Users/solver/.julia/packages/ClassicalOrthogonalPolynomials/Ak9Ta/src/adaptivetransform.jl:1
>1  transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  =
 2      adaptivetransform_ldiv(convert(AbstractQuasiArray{promote_type(T,V)}, A), f)

In adaptivetransform_ldiv(A, f) at /Users/solver/.julia/packages/ClassicalOrthogonalPolynomials/Ak9Ta/src/adaptivetransform.jl:16
 16  function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiVector{V}) where {U,V}