allow chebinterp to omit the order and determine it adaptively (#32)

* allow chebinterp to omit the order and determine it adaptively, similar to ApproxFun

* use powers of 2 for repeatedly doubled orders to speed up DCTs

Author: stevengj

README.md

@@ -21,7 +21,17 @@ via `c(y)`.  Alternatively, you can combine `chebpoints` and `chebinterp` into a
 ```jl
 c = chebinterp(f, order, lb, ub)
 ```
-which evaluates the given function at Chebyshev points for you.
+which evaluates the given function at Chebyshev points for you.  In fact, you can omit the `order`
+and simply call
+```jl
+c = chebinterp(f, lb, ub)
+```
+and it will adaptively double the order until a given relative tolerance is attained
+(specified by the optional `tol` keyword argument, default to `5*eps` in the
+given floating-point precision).  The starting order for the doubling can be specified
+by a `min_order` keyword (which defaults to `8` in each dimension), and a maximum
+can be specified by the `max_order` keyword.  This feature is best for *smooth* functions.
+(More sophisticated versions of this idea are implemented by [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl).)
 
 We also provide a function `chebgradient(c,y)` that returns a tuple `(c(y), ∇c(y))` of
 the interpolant and its gradient at a point `y`.  (You can also use automatic differentiation, e.g. via the [ForwardDiff.jl package](https://github.com/JuliaDiff/ForwardDiff.jl),

src/interp.jl

@@ -37,6 +37,8 @@ chebpoints(order::Integer, lb::Real, ub::Real) =
 
 # O(n log n) method to compute Chebyshev coefficients
 function chebcoefs(vals::AbstractArray{<:Number,N}) where {N}
+    all(isfinite, vals) || throw(DomainError("non-finite interpolant value"))
+
      # type-I DCT, except for size-1 dimensions where we want identity
     kind = map(n -> n > 1 ? FFTW.REDFT00 : FFTW.DHT, size(vals))
     coefs = FFTW.r2r(vals, kind)
@@ -73,7 +75,7 @@ infnorm(x::Number) = abs(x)
 infnorm(x::AbstractArray) = maximum(abs, x)
 
 function droptol(coefs::Array{<:Any,N}, tol::Real) where {N}
-    abstol = sum(infnorm, coefs) * tol # absolute tolerance = L1 norm * tol
+    abstol = maximum(infnorm, coefs) * tol # absolute tolerance = Linf norm * tol
     all(c -> infnorm(c) ≥ abstol, coefs) && return coefs # nothing to drop
 
     # compute the new size along each dimension by checking
@@ -103,6 +105,7 @@ epsbounds(lb, ub) = eps(float(real(promote_type(eltype(lb), eltype(ub)))))
 """
     chebinterp(vals, lb, ub; tol=eps)
     chebinterp(f::Function, order, lb, ub; tol=eps)
+    chebinterp(f::Function, lb, ub; tol=eps, min_order=(8,..), max_order=(typemax(Int),...))
 
 Given a multidimensional array `vals` of function values (at
 points corresponding to the coordinates returned by `chebpoints`),
@@ -115,6 +118,12 @@ or a tuple of integers), and it will call [`chebpoints`](@ref) for you
 to obtain the Chebyshev points and then compute `vals` by evaluating
 `f` at these points.
 
+If a function `f` is supplied and the `order` argument is omitted, it will adaptively
+determine the order by repeatedly doubling it until `tol` is achieved
+or `max_order` is reached, starting at `min_order` (which defaults to `8` or
+a tuple of `8`s in each dimension; this might need to be increased for
+highly oscillatory functions).   This feature is best used for smooth functions.
+
 This object `c = chebinterp(...)` can be used to
 evaluate the interpolating polynomial at a point `x` via
 `c(x)`.
@@ -146,6 +155,11 @@ of vectors are painful to construct.)
 chebinterp_v1(vals::AbstractArray{T}, lb, ub; tol::Real=epsvals(vals)) where {T<:Number} =
     chebinterp(dropdims(reinterpret(SVector{size(vals,1),T}, Array(vals)), dims=1), lb, ub; tol=tol)
 
+#####################################################################################################
+# supplying the function rather than the values
+
+# supplying both function and order:
+
 chebinterp(f::Function, order::NTuple{N,Int}, lb::SVector{N,<:Real}, ub::SVector{N,<:Real}; tol::Real=epsbounds(lb,ub)) where {N} =
     chebinterp(map(f, chebpoints(order, lb, ub)), lb, ub; tol=tol)
 
@@ -154,3 +168,36 @@ chebinterp(f::Function, order::NTuple{N,Int}, lb::AbstractArray{<:Real}, ub::Abs
 
 chebinterp(f::Function, order::Integer, lb::Real, ub::Real; tol::Real=epsbounds(lb,ub)) =
     chebinterp(x -> f(@inbounds x[1]), (Int(order),), SVector(lb), SVector(ub); tol=tol)
+
+## adaptively determine the order by repeated doublying, ala chebfun or approxfun:
+
+function chebinterp(f::Function, lb::SVector{N,<:Real}, ub::SVector{N,<:Real};
+                    tol::Real=5epsbounds(lb,ub),
+                    min_order::NTuple{N,Int}=ntuple(i->8,Val{N}()),
+                    max_order::NTuple{N,Int}=ntuple(i->typemax(Int),Val{N}())) where {N}
+    tol > 0 || throw(ArgumentError("tolerance $tol must be > 0"))
+    all(min_order .> 0) || throw(ArgumentError("minimum order $min_order must be > 0"))
+    all(max_order .>= 0) || throw(ArgumentError("maximum order $max_order must be ≥ 0"))
+    order = min.(min_order, max_order)
+    while true
+        # in principle we could re-use function evaluations when doubling the order,
+        # but that would greatly complicate the code and only saves a factor of 2
+        c = chebinterp(map(f, chebpoints(order, lb, ub)), lb, ub; tol=tol)
+        order_done = (size(c.coefs) .- 1 .< order) .| (order .== max_order)
+        all(order_done) && return c
+        order = ifelse.(order_done, order, min.(max_order, nextpow.(2, order .* 2)))
+    end
+end
+
+function chebinterp(f::Function, lb::AbstractArray{<:Real}, ub::AbstractArray{<:Real};
+                    tol::Real=5epsbounds(lb,ub),
+                    min_order=fill(8, length(lb)), max_order=fill(typemax(Int), length(lb)))
+    N = length(lb)
+    N == length(ub) == length(min_order) == length(max_order) || throw(DimensionMismatch("dimensions must all == $N"))
+    Base.require_one_based_indexing(min_order, max_order)
+    chebinterp(f, SVector{N}(lb), SVector{N}(ub);
+               tol=tol, min_order=ntuple(i -> Int(min_order[i]), N), max_order=ntuple(i -> Int(max_order[i]), N))
+end
+
+chebinterp(f::Function, lb::Real, ub::Real; tol::Real=5epsbounds(lb,ub), min_order::Integer=8, max_order::Integer=typemax(Int)) =
+    chebinterp(x -> f(@inbounds x[1]), SVector(lb), SVector(ub); tol=tol, min_order=(Int(min_order),), max_order=(Int(max_order),))

src/roots.jl

@@ -88,7 +88,7 @@ to `5eps` where `eps` is the precision of `c`.
 function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Integer=50)
     tol > 0 || throw(ArgumentError("tolerance $tol for truncating coefficients must be > 0"))
     maxsize > 0 || throw(ArgumentError("maxsize $maxsize must be > 0"))
-    abstol = sum(abs, c.coefs) * tol # absolute tolerance = L1 norm * tol
+    abstol = maximum(abs, c.coefs) * tol # absolute tolerance = Linf norm * tol
     n = something(findlast(c -> abs(c) ≥ abstol, c.coefs), 1)
     if n <= maxsize
         λ = hesseneigvals!(UpperHessenberg(_colleague_matrix(@view c.coefs[1:n])))

test/runtests.jl

@@ -25,11 +25,13 @@ Random.seed!(314159) # make chainrules tests deterministic
         test_frule(interp, x1, rtol=sqrt(eps(T)), atol=sqrt(eps(T)))
         test_rrule(interp, x1, rtol=sqrt(eps(T)), atol=sqrt(eps(T)))
 
-        @test roots(chebinterp(x -> 1.1, 10, 0,1)) ≈ Float64[]
-        @test roots(chebinterp(x -> x - 0.1, 10, 0,1)) ≈ [0.1]
-        @test roots(chebinterp(x -> (x - 0.1)*(x-0.2), 10, 0,1)) ≈ [0.1, 0.2]
-        @test roots(chebinterp(cos, 10000, 0, 1000))/pi ≈ (0:317) .+ T(0.5)
-        @test colleague_matrix(chebinterp(x -> (x - 0.1)*(x-0.2), 10, 0,1)) ≈ T[0.5 -0.24; 0.5 -0.2]
+        @test roots(chebinterp(x -> 1.1, 0,1)) ≈ Float64[]
+        @test roots(chebinterp(x -> x - 0.1, 0,1)) ≈ [0.1]
+        p = chebinterp(x -> (x - 0.1)*(x-0.2), 0,1)
+        @test length(p.coefs) == 3 # degree-2 polynomial
+        @test roots(p) ≈ [0.1, 0.2]
+        @test roots(chebinterp(cos, 10^4, 0, 1000))/pi ≈ (0:317) .+ T(0.5)
+        @test colleague_matrix(chebinterp(x -> (x - 0.1)*(x-0.2), 0,1)) ≈ T[0.5 -0.24; 0.5 -0.2]
     end
 end
 
@@ -42,13 +44,16 @@ end
         @test eltype(x) == SVector{2,T}
         interp = chebinterp(f.(x), lb, ub)
         interp2 = chebinterp(f, (48,39), lb, ub)
+        interp3 = chebinterp(f, lb, ub)
+        @test all(size(interp3.coefs) .< (40,30))
+        @test all(size(chebinterp(f, lb, ub, max_order=(18, 4)).coefs) .≤ (19,5))
         @test interp isa FastChebInterp.ChebPoly{2,T,T}
         @test interp2.coefs ≈ interp.coefs
         interp0 = chebinterp(f.(x), lb, ub, tol=0)
         @test repr("text/plain", interp0) == "ChebPoly{2,$T,$T} order (48, 39) polynomial on [-0.3,0.9] × [0.1,1.2]"
         @test ndims(interp) == 2
         x1 = T[0.2, 0.8] # not too close to boundary or test_frule can fail by taking a big FD step
-        @test interp(x1) ≈ f(x1)
+        @test interp(x1) ≈ f(x1) ≈ interp3(x1)
         @test interp(x1) ≈ interp0(x1) rtol=10eps(T)
         @test all(n -> n[1] < n[2], zip(size(interp.coefs), size(interp0.coefs)))
         @test chebgradient(interp, x1) ≈′ (f(x1), ∇f(x1))