Move the point-evaluation kernel (src/point_eval.jl) upstream to BasisMatrices.jl

[oyamad]

Background

#93 added src/point_eval.jl: a non-allocating single-point interpolant evaluation kernel used in the inner loop of _s_wise_max!. Previously, each objective-function call inside Optim.maximize invoked funeval, which has no fast path for a single point: the scalar method wraps the point in a matrix (fill(x, 1, 1)) and constructs a full BasisMatrix. Replacing those calls with the kernel reduced the memory footprint of a 50-iteration VFI from 700 MiB / 5.8M allocations to 434 KiB / 5.4K allocations and roughly halved wall-clock time (see the benchmark table in #93; addresses #73 and most of the cost behind #74).

The file is deliberately self-contained — it contains no DP-specific types and speaks only BasisMatrices' language — because this functionality conceptually belongs upstream: it is generic interpolation machinery, useful to anyone evaluating a fitted function point-wise in a loop (optimizers, solvers, simulations). This issue tracks moving it upstream to BasisMatrices.jl.

What the kernel provides (order 0, implemented)

• PointEvalCache(p::BasisParams): per-dimension cache with a preallocated buffer for basis-function values at a point; specialized kernels for ChebParams (Chebyshev recurrence on the unscaled point), SplineParams (de Boor recurrence on the augmented knot sequence, replicating lookup(augbreaks, x, 3) semantics), and LinParams (two-weight linear interpolation, both the evennum fast path and the lookup path); other basis families fall back to a generic (allocating) evalbase-based implementation.
• point_evalbase!(cache, x::Real) -> (first, nvals): writes the values of the basis functions with nonzero support at x into cache.vals, exploiting sparsity for splines (k+1 values) and piecewise linear (2 values).
• FunEvalCache(basis::Basis{N}) and funeval_point!(cache, c, x) -> Float64: tensor-product contraction over the per-dimension kernels for an N-dimensional point; agrees with funeval(c, basis, x) to machine precision, including for points outside the interpolation domain (each family reproduces evalbase's extrapolation behavior exactly).
• Agreement tests in test/test_point_eval.jl: 260 cases across basis families (Cheb; splines of degree 1–3 with even and uneven breaks; Lin with evennum zero and nonzero), mixed 2-D/3-D tensor bases, points at/inside/outside domain boundaries, plus zero-allocation assertions. These tests would transplant to BasisMatrices.jl essentially as-is, and until the move happens they guard against divergence from upstream funeval conventions.

Design decisions on record

• The caches are Float64-only and funeval_point! restricts coefficients to AbstractVector{Float64}, so that use with other numeric types (BigFloat, AD dual numbers) fails at dispatch instead of silently downcasting. Generic-eltype support would mean parameterizing the caches; whether that is worth the complexity is an API question best settled upstream.
• Only derivative order 0 is implemented. The (first, nvals)-into-buffer contract extends directly to derivative orders via the standard Chebyshev/B-spline derivative recurrences (e.g. point_evalbase!(cache, x, order)); this is what derivative-based inner optimization (Optimization #7) needs. Since the interpolant's derivative is known in closed form, the intended pattern there is chain-rule composition — AD applied to the user's f and g only, with the interpolant gradient supplied by an order-1 kernel — rather than running AD through the basis recurrences.
• FunEvalCache is mutable scratch and not thread-safe; parallel state loops need one cache per thread (relevant to a future threading workspace, and to any user-facing result type that embeds a cache, cf. set_eval_nodes! #10).

Suggested upstream shape

Beyond a direct port, two natural integrations on the BasisMatrices.jl side:

• An Interpoland fast path, so itp(x) at a single point uses the kernel instead of building a BasisMatrix — this also connects to the old discussion in Interp #8 about moving ContinuousDPs' interpolation machinery upstream.
• Complementarity with In-Place Basis Matrix and funeval BasisMatrices.jl#60, which requests in-place batch updates of basis matrices at a fixed set of nodes: that covers the fixed-grid batch workflow, this kernel covers the scattered point-wise workflow where the evaluation points change on every call; a design accommodating both may be worth considering.

Plan

1. Propose the move at BasisMatrices.jl (issue link to follow) and iterate on API naming/shape there.
2. Once the functionality is available in a released BasisMatrices.jl version: delete src/point_eval.jl and test/test_point_eval.jl here, import the upstream equivalents in src/cdp.jl, and bump the BasisMatrices compat bound.

🤖 Drafted with Claude Code (Claude Fable 5)

[oyamad]

Update since opening:

• Derivative point evaluation is now implemented downstream as well (ENH: Add derivative point evaluation via coefficient differentiation #98: DerivFunEvalCache, via coefficient differentiation using derivative_op), so the upstream proposal covers derivative orders with a working implementation, not just a plan.
• A finding relevant to the upstream design (Proposal: non-allocating single-point evaluation (funeval fast path) BasisMatrices.jl#82): funeval's Direct/SplineSparse path errors for LinParams derivatives ("derivatives un-supported right now"), while the dense evalbase supports them.
• The kernels now also back the inner solvers (PERF: Solve the inner maximization via the first-order condition #99, ENH: Add multi-dimensional continuous actions #105) and the policy-evaluation assembly (PERF: Assemble the policy-evaluation system with point kernels, sparse when possible #100), which raises the value of the move; the set to migrate is unchanged: src/point_eval.jl plus its agreement tests.