new chebinterp(f::Function, order, lb, ub) method (#30)

Author: stevengj

README.md

@@ -17,7 +17,11 @@ x = chebpoints(order, lb, ub) # an array of `SVector` from [StaticArrays.jl](htt
 c = chebinterp(f.(x), lb, ub)
 ```
 and then evaluate the interpolant for a point `y` (a vector)
-via `c(y)`.
+via `c(y)`.  Alternatively, you can combine `chebpoints` and `chebinterp` into a single call:
+```jl
+c = chebinterp(f, order, lb, ub)
+```
+which evaluates the given function at Chebyshev points for you.
 
 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),
@@ -49,6 +53,7 @@ Here is an example interpolating the (highly oscillatory) 1d function `f(x) = si
 f(x) = sin(2x + 3cos(4x))
 x = chebpoints(200, 0, 10)
 c = chebinterp(f.(x), 0, 10)
+# alternatively, c = chebinterp(f, 200, 0, 10)
 ```
 We can then compare the exact function and its interpolant at a set of points:
 ```jl
@@ -90,6 +95,7 @@ g(x) = sin(x[1] + cos(x[2]))
 lb, ub = [1,3], [2, 4] # lower and upper bounds of the domain, respectively
 x = chebpoints((10,20), lb, ub)
 c = chebinterp(g.(x), lb, ub)
+# alternatively, c = chebinterp(g, (10,20), lb, ub)
 ```
 Let's evaluate the interpolant at an arbitrary point `(1.2, 3.4)` and compare it to the exact value:
 ```jl

src/interp.jl

@@ -6,8 +6,10 @@ function chebpoint(i::CartesianIndex{N}, order::NTuple{N,Int}, lb::SVector{N}, u
     @. lb + (1 + cos(T($SVector($Tuple(i))) * π / $SVector(ifelse.(iszero.(order),2,order)))) * (ub - lb) * $(T(0.5))
 end
 
-chebpoints(order::NTuple{N,Int}, lb::SVector{N}, ub::SVector{N}) where {N} =
-    [chebpoint(i,order,lb,ub) for i in CartesianIndices(map(n -> n==0 ? (1:1) : (0:n), order))]
+function chebpoints(order::NTuple{N,Int}, lb::SVector{N}, ub::SVector{N}) where {N}
+    all(≥(0), order) || throw(ArgumentError("invalid negative order $order"))
+    return [chebpoint(i,order,lb,ub) for i in CartesianIndices(map(n -> n==0 ? (1:1) : (0:n), order))]
+end
 
 """
     chebpoints(order, lb, ub)
@@ -26,13 +28,12 @@ the order in that direction.)
 function chebpoints(order, lb, ub)
     N = length(order)
     N == length(lb) == length(ub) || throw(DimensionMismatch())
-    all(≥(0), order) || throw(ArgumentError("invalid negative order $order"))
     return chebpoints(NTuple{N,Int}(order), SVector{N}(lb), SVector{N}(ub))
 end
 
 # return array of scalars in 1d
 chebpoints(order::Integer, lb::Real, ub::Real) =
-    first.(chebpoints(order, SVector(lb), SVector(ub)))
+    first.(chebpoints((Int(order),), SVector(lb), SVector(ub)))
 
 # O(n log n) method to compute Chebyshev coefficients
 function chebcoefs(vals::AbstractArray{<:Number,N}) where {N}
@@ -97,17 +98,24 @@ end
 
 # precision for float(vals), handling arrays of numbers and arrays of arrays of numbers
 epsvals(vals) = eps(float(real(eltype(eltype(vals)))))
+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)
 
 Given a multidimensional array `vals` of function values (at
 points corresponding to the coordinates returned by `chebpoints`),
 and arrays `lb` and `ub` of the lower and upper coordinate bounds
 of the domain in each direction, returns a Chebyshev interpolation
 object (a `ChebPoly`).
 
-This object `c = chebinterp(vals, lb, ub)` can be used to
+Alternatively, one can supply a function `f` and an `order` (an integer
+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.
+
+This object `c = chebinterp(...)` can be used to
 evaluate the interpolating polynomial at a point `x` via
 `c(x)`.
 
@@ -137,3 +145,12 @@ 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)
+
+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)
+
+chebinterp(f::Function, order::NTuple{N,Int}, lb::AbstractArray{<:Real}, ub::AbstractArray{<:Real}; tol::Real=epsbounds(lb,ub)) where {N} =
+    chebinterp(f, order, SVector{N}(lb), SVector{N}(ub); tol=tol)
+
+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)

test/runtests.jl

@@ -14,11 +14,13 @@ Random.seed!(314159) # make chainrules tests deterministic
         x = chebpoints(48, lb, ub)
         @test eltype(x) == T
         interp = chebinterp(f.(x), lb, ub, tol=0)
+        interp2 = chebinterp(f, 48, lb, ub, tol=0)
         @test interp isa FastChebInterp.ChebPoly{1,T,T}
+        @test interp2.coefs ≈ interp.coefs
         @test repr("text/plain", interp) == "ChebPoly{1,$T,$T} order (48,) polynomial on [-0.3,0.9]"
         @test ndims(interp) == 1
         x1 = T(0.2)
-        @test interp(x1) ≈ f(x1)
+        @test interp(x1) ≈ f(x1) ≈ interp2(x1)
         @test chebgradient(interp, x1) ≈′ (f(x1), f′(x1))
         test_frule(interp, x1, rtol=sqrt(eps(T)), atol=sqrt(eps(T)))
         test_rrule(interp, x1, rtol=sqrt(eps(T)), atol=sqrt(eps(T)))
@@ -33,7 +35,9 @@ end
         x = chebpoints((48,39), lb, ub)
         @test eltype(x) == SVector{2,T}
         interp = chebinterp(f.(x), lb, ub)
+        interp2 = chebinterp(f, (48,39), lb, ub)
         @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