juliaapproximation / approxfunorthogonalpolynomials.jl Goto Github PK
View Code? Open in Web Editor NEWSupport for orthogonal polynomial-based spaces in ApproxFun
License: MIT License
approxfunorthogonalpolynomials.jl's People
Contributors
![dependabot[bot] avatar](https://avatars.githubusercontent.com/in/29110?v=4 "dependabot[bot]")
![github-actions[bot] avatar](https://avatars.githubusercontent.com/in/15368?v=4 "github-actions[bot]")
Stargazers
Watchers
Forkers
putianyi889 tomtrogdon gyoshi standardgalactic luapulu stevengj jishnub benjione lucasaschenbach jondeuceapproxfunorthogonalpolynomials.jl's Issues
Refactor tests
We should have as goal a 10-minute test suite per job on CI, per JuliaApproximation package. The tests in the ApproxFun ecosystem evolved from
- a collection of tricksy gotcha's to make function approximation fun and intuitive
- systematic test functions from https://github.com/JuliaApproximation/ApproxFunBase.jl/blob/master/src/testing.jl Probably there are ways to pare down the tests while ensuring the same level of support. Huge test suites needlessly clog up CI infrastructure and slow down development.
Sizes of Multiplication operators in Ultraspherical spaces vs others are different infinities
julia> size(Multiplication(Fun(), Chebyshev()))
(ℵ₀, ℵ₀)
julia> size(Multiplication(Fun(), Legendre()))
(ℵ₀, ℵ₀)
julia> size(Multiplication(Fun(), Ultraspherical(1)))
(∞, ∞)
julia> size(Multiplication(Fun(), Ultraspherical(2)))
(ℵ₀, ℵ₀)This seems to happen because the Toeplitz-Hankel operators are defined on Sequence spaces, and
julia> ApproxFunBase.dimension(SequenceSpace())
∞However, this seems to be an implementation detail. Should size use the dimension of the domain/range spaces instead?
julia> ApproxFunBase.dimension(Ultraspherical(1))
ℵ₀Ultraspherical conversion to lower order domain error consistency
Hi there, I was going through some examples and trying to learn how to work with ApproxFun when I ran into this.
I am not sure if this is a superficial issue of invalid inputs not being caught, or if there is some deeper problems, but here are my observations:
Note: The issue is the Conversion() method, but I use an Operator with → here because it results in an immediately evaluated representation. The following issues occurred with all the Operators I tried.
I presume the following should not be allowed to evaluate
julia> Operator(Fun()) : Ultraspherical(2) → Ultraspherical(1)
TimesOperator : Ultraspherical(2) → Ultraspherical(1)
NaN NaN 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱ ⋱and a stack overflow error is not ideal for the following
julia> Operator(Fun()) : Ultraspherical(3) → Ultraspherical(1.9)
ERROR: StackOverflowError:
Stacktrace:
[1] Conversion(::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunOrthogonalPolynomials\jojxS\src\Spaces\Ultraspherical\UltrasphericalOperators.jl:141
[2] Conversion(::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunOrthogonalPolynomials\jojxS\src\Spaces\Ultraspherical\UltrasphericalOperators.jl:148 (repeats 14488 times)
[3] promoterangespace(::ApproxFunBase.ConcreteMultiplication{Chebyshev{ChebyshevInterval{Float64},Float64},Ultraspherical{Int64,ChebyshevInterval{Float64},Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}, ::Ultraspherical{Int64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunBase\FYFgS\src\Operators\spacepromotion.jl:83
[4] promoterangespace at C:\Users\Johannes\.julia\packages\ApproxFunBase\FYFgS\src\Operators\spacepromotion.jl:79 [inlined]
[5] →(::ApproxFunBase.ConcreteMultiplication{Chebyshev{ChebyshevInterval{Float64},Float64},Ultraspherical{Int64,ChebyshevInterval{Float64},Float64},Float64}, ::Ultraspherical{Float64,ChebyshevInterval{Float64},Float64}) at C:\Users\Johannes\.julia\packages\ApproxFunBase\FYFgS\src\Operators\spacepromotion.jl:74
[6] top-level scope at none:0If we want to catch these two cases, I believe changing the elseifs as in the following should be sufficient?
function Conversion(A::Ultraspherical,B::Ultraspherical)
a=order(A); b=order(B)
if b==a
ConversionWrapper(Operator(I,A))
# elseif a<b≤a+1 || b<a≤b+1
elseif -1≤b-a≤1 && b≠1
ApproxFun.ConcreteConversion(A,B)
# elseif b ≠ 1
elseif b-a>1
d=domain(A)
ConversionWrapper(TimesOperator(Conversion(Ultraspherical(b-1,d),B),
Conversion(A,Ultraspherical(b-1,d))))
else
error("Cannot convert from $A to $B")
end
endI also noted some inconsistent behaviour of getindex() when passing Ultrasphericals of Float and/or Rational order (inconsistency of when it throws an error, not inconsistencies in evaluation), but I suppose that is a separate issue.
TagBot trigger issue
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!
`points(Ultraspherical(1, -1..1), n)` returns points on Chebyshev grid of the first kind and not the second kind
using ApproxFunOrthogonalPolynomials
# For reference
chebyshevpoints(Float64, 3, Val(1)) # [0.866, 0.0, -0.866]
chebyshevpoints(Float64, 3, Val(2)) # [1.0, 0.0, -1.0]
# These should be points corresponding to the Chebyshev polynomials of the second kind, I thought?
points(Ultraspherical(1, -1..1), 3) # [0.866, 0.0, -0.866]ApproxFun+IntervalArithmetic: Add support for constructor
The following would ideally work:
julia> using ApproxFun
julia> import IntervalArithmetic: @interval
julia> Base.sinpi(x::IntervalArithmetic.Interval) = sin(π*x)
julia> Fun(exp, @interval(0)..@interval(1))
ERROR: MethodError: no method matching plan_chebyshevtransform(::Array{IntervalArithmetic.Interval{Float64},1})
Closest candidates are:
plan_chebyshevtransform(::Array{D<:DualNumbers.Dual,1}; kind) where D<:DualNumbers.Dual at /Users/solver/Projects/ApproxFun.jl/src/Extras/dualnumbers.jl:35
plan_chebyshevtransform(::AbstractArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},1}; kind) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/solver/.julia/packages/FastTransforms/vEjxF/src/chebyshevtransform.jl:29
plan_chebyshevtransform(::AbstractArray{T<:Union{Complex{BigFloat}, BigFloat},1}; kind) where T<:Union{Complex{BigFloat}, BigFloat} at /Users/solver/Projects/ApproxFun.jl/src/Extras/fftGeneric.jl:49
Stacktrace:
[1] plan_transform(::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Array{IntervalArithmetic.Interval{Float64},1}) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:65
[2] transform(::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Array{IntervalArithmetic.Interval{Float64},1}) at /Users/solver/Projects/ApproxFunBase.jl/src/Space.jl:427
[3] default_Fun(::Type{IntervalArithmetic.Interval{Float64}}, ::Function, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Array{IntervalArithmetic.Interval{Float64},1}, ::Type{Val{false}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:44
[4] default_Fun(::ApproxFunBase.DFunction, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Int64, ::Type{Val{false}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:57
[5] default_Fun(::Function, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}, ::Int64) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:72
[6] default_Fun(::ApproxFunBase.DFunction, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:121
[7] Fun(::Function, ::Chebyshev{Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}},IntervalArithmetic.Interval{Float64}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:177
[8] Fun(::Function, ::Interval{:closed,:closed,IntervalArithmetic.Interval{Float64}}) at /Users/solver/Projects/ApproxFunBase.jl/src/constructors.jl:173
[9] top-level scope at none:0Ideally, the transform would return a data structure representing an infinite vector with exponential decay. It's actually already possible to construct such a thing using InfiniteArrays.jl:
julia> BroadcastArray(IntervalArithmetic.Interval, -2.0 .^(-(1:∞)), 2.0 .^(-(1:∞)))
BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}} with indices OneToInf():
[-0.5, 0.5]
[-0.25, 0.25]
[-0.125, 0.125]
[-0.0625, 0.0625]
[-0.03125, 0.03125]
[-0.015625, 0.015625]
[-0.0078125, 0.0078125]
[-0.00390625, 0.00390625]
[-0.00195313, 0.00195313]
[-0.000976563, 0.000976563]
[-0.000488282, 0.000488282]
[-0.000244141, 0.000244141]
[-0.000122071, 0.000122071]
[-6.10352e-05, 6.10352e-05]
[-3.05176e-05, 3.05176e-05]
[-1.52588e-05, 1.52588e-05]
[-7.6294e-06, 7.6294e-06]
[-3.8147e-06, 3.8147e-06]
⋮ This would be combined with a finite vector of coefficients using Vcat. Some mathematical thought is needed on how to calculate this automatically, which would require bounding in a Bernstein ellipse.
@dpsanders FYI.
Spurious test failure
This test failure showed up in a log:
Check resizing: Test Failed at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:226
Expression: norm(((Laplacian() * u + 100u) - (A * u)[2]).coefficients) < 1.0e-9
Evaluated: 1.3708085290831766e-9 < 1.0e-9
Stacktrace:
[1] top-level scope at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:226
[2] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
[3] top-level scope at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:176
[4] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
[5] top-level scope at /root/.julia/packages/ApproxFunOrthogonalPolynomials/Cxm84/test/PDETest.jl:5
I don't think it is a real error so it might make sense to loosen the tolerance a bit.
ApproxFun+IntervalArithmetic: Support evaluation with infinite vectors
It's possible to represent functions with infinite coefficients using InfiniteArrays.jl, however, evaluation is broken:
julia> f = Fun(Chebyshev(), -2.0 .^(-(1:∞)))
Fun(Chebyshev(),[-0.5, -0.25, -0.125, -0.0625, -0.03125, -0.015625, -0.0078125, -0.00390625, -0.00195313, -0.000976563 … ])
julia> f(0.1)
ERROR: ArgumentError: Cannot create range starting at infinity
Stacktrace:
[1] (::Colon)(::InfiniteArrays.Infinity, ::Int64, ::Int64) at /Users/solver/Projects/InfiniteArrays.jl/src/infrange.jl:10
[2] chebyshev_clenshaw(::BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/LinearAlgebra/clenshaw.jl:243
[3] clenshaw(::Chebyshev{ChebyshevInterval{Float64},Float64}, ::BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:73
[4] evaluate at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/PolynomialSpace.jl:23 [inlined]
[5] evaluate at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:216 [inlined]
[6] (::Fun{Chebyshev{ChebyshevInterval{Float64},Float64},Float64,BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}})(::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:220
[7] top-level scope at none:0It's actually a bit simpler to think about in the interval arithmetic setting where we only need a bound on evaluation (otherwise we would need to talk about tolerances). That is, in the format that #2 would return:
julia> f = Fun(Chebyshev(), BroadcastArray(IntervalArithmetic.Interval, -2.0 .^(-(1:∞)), 2.0 .^(-(1:∞))))
Fun(Chebyshev(),IntervalArithmetic.Interval{Float64}[[-0.5, 0.5], [-0.25, 0.25], [-0.125, 0.125], [-0.0625, 0.0625], [-0.03125, 0.03125], [-0.015625, 0.015625], [-0.0078125, 0.0078125], [-0.00390625, 0.00390625], [-0.00195313, 0.00195313], [-0.000976563, 0.000976563] … ])
julia> f(0.1)
ERROR: ArgumentError: Cannot create range starting at infinity
Stacktrace:
[1] (::Colon)(::InfiniteArrays.Infinity, ::Int64, ::Int64) at /Users/solver/Projects/InfiniteArrays.jl/src/infrange.jl:10
[2] chebyshev_clenshaw(::BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/LinearAlgebra/clenshaw.jl:243
[3] clenshaw(::Chebyshev{ChebyshevInterval{Float64},Float64}, ::BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:73
[4] evaluate at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/PolynomialSpace.jl:23 [inlined]
[5] evaluate at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:216 [inlined]
[6] (::Fun{Chebyshev{ChebyshevInterval{Float64},Float64},IntervalArithmetic.Interval{Float64},BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}})(::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:220
[7] top-level scope at none:0For Chebyshev it should be easy to work out bounds and override chebyshev_clenshaw(::BroadcastArray{...}, x) to return these bounds. We could also then override chebyshev_clenshaw(a::Vcat, x) to do Clenshaw recursively, using the result of chebyshev_clenshaw(a.arrays[end], x) to determine the initial condition for the rest of Clenshaw. This would give rigorous evaluation of functions.
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.