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

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 (#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. #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 #8 about moving ContinuousDPs' interpolation machinery upstream.
• Complementarity with QuantEcon/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)

Start by reviewing src/point_eval.jl and test/test_point_eval.jl, then inspect how src/cdp.jl imports and uses the kernel. Propose and iterate on the API with BasisMatrices.jl, including the Interpoland integration questions. Once an upstream release provides the functionality, remove the local implementation and tests, import the upstream equivalents, and bump the BasisMatrices compat bound.