Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Controlling tolerance during function approximation #53

Open
wsshin opened this issue Jan 22, 2022 · 1 comment
Open

Controlling tolerance during function approximation #53

wsshin opened this issue Jan 22, 2022 · 1 comment

Comments

@wsshin
Copy link

wsshin commented Jan 22, 2022

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.

@wsshin wsshin changed the title Controlling tolerance in during function approximation Controlling tolerance during function approximation Jan 22, 2022
@dlfivefifty
Copy link
Member

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}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants