Skip to content

ApproxFun+IntervalArithmetic: Support evaluation with infinite vectors #3

@dlfivefifty

Description

@dlfivefifty

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:0

It'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:0

For 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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions