Design review report

[kdw503]

---

Phase 3 — Coherence checks

3a. Scope coherence — The package is narrow and well-focused: one factory function, one type, quadratic B-spline caching. No unrelated utilities. One mild
question: the "quadratic B-spline" algorithm is hardwired — the caching wrapper doesn't allow substituting a different interpolation kernel. Whether
that's a scope constraint or an oversight depends on intent.

3b. Type hierarchy — CachedInterpolation{T, N, M, O, K} has five type parameters with a silent constraint M = N + K. O (origin offset) is packed into the
type signature despite being semantically a runtime value; its only apparent purpose is to propagate into axes without storing a field. The
AbstractInterpolation inheritance is load-bearing (it forces getindex to exist) but may be more obligation than benefit for a caching wrapper.

3c. Overlapping operations — itp(xs...) and itp[xs...] do the same thing for integer inputs; the latter just forwards to the former. The distinction is
implicit and underdocumented.

3d. Implied but absent operations — gradient requires a prior call to itp(ys...) at the same coordinates; calling it in isolation silently produces wrong
results. There is no combined evaluate_and_gradient function that would make this safe by construction. gradient! delegates to gradient and then copies,
so the in-place variant provides no allocation saving over the out-of-place variant.

3c. Overlapping operations — itp(xs...) and itp[xs...] do the same thing for integer inputs; the latter just forwards to the former. The distinction is
implicit and underdocumented.

3d. Implied but absent operations — gradient requires a prior call to itp(ys...) at the same coordinates; calling it in isolation silently produces wrong
results. There is no combined evaluate_and_gradient function that would make this safe by construction. gradient! delegates to gradient and then copies,
so the in-place variant provides no allocation saving over the out-of-place variant.

3e. Abstraction level — Consistent. Everything exported is directly user-facing.

3f. Composability — The shared coefs buffer across all tiles creates a non-obvious constraint: concurrent access to different tiles of the same
cachedinterpolators array is data-race unsafe. This is not documented. The gradient precondition further breaks composability — you cannot safely call
gradient(itp, ys...) without knowing what itp was last evaluated at.

3g. Relationship to Base/stdlib — Extends Interpolations.gradient/gradient! and subtypes AbstractInterpolation. No collision with Base names. The
inherited AbstractArray obligations (size, axes, getindex) are satisfied, though getindex for non-integer inputs is not clearly defined.

3h. Error and failure model — The gradient precondition is silently violated (wrong results, not an error). Boundary behavior (evaluation near the edge of
the array where a 3-element neighborhood would go out of bounds) is not mentioned.

3i. Missing public annotations — The inner cachedinterpolators function-barrier method is not exported or annotated public, appears in no user-facing
documentation, and is only visible in tests via qualified access (CachedInterpolations.cachedinterpolators).

---

Phase 4 — Report

---

Likely design issues

1. Gradient precondition is an invisible footgun (src/CachedInterpolations.jl, gradient method).
   Interpolations.gradient(itp, ys...) is documented to require a prior itp(ys...) call at the same coordinates. Nothing in the type system or the call
   signature enforces this. A user who reaches for gradient in isolation — a completely natural thing to do with any AbstractInterpolation — gets silently
   wrong results. The standard fix in similar packages is a combined value_and_gradient function (or making the evaluation-then-gradient sequence atomic), so
   that the precondition cannot be violated.

2. gradient! is not actually in-place (gradient! method).
   The implementation delegates to gradient and then copies into g. This means it allocates just as much as the out-of-place version, yet it presents the !
   interface that conventionally signals "no allocation." A user who calls gradient! expecting to avoid allocations in a hot loop will be silently
   disappointed.

3. Silent constraint M = N + K in type parameters.
   CachedInterpolation{T, N, M, O, K} carries both M and N/K even though M is determined by N + K. Julia cannot enforce this relationship at the type level,
   so a user constructing one directly (rather than via cachedinterpolators) can create an incoherent instance. If M is needed at all, it could be computed
   inside the struct via a helper rather than appearing as a free parameter.

4. The shared coefs buffer is undocumented.
   All CachedInterpolation objects produced by one cachedinterpolators call share a single coefs buffer. Accessing two tiles alternately will thrash the
   cache. Accessing them from two threads concurrently is a data race. Neither constraint appears in the documentation or in the type's fields/docstring.

---

Design questions

Q1. Why subtype AbstractInterpolation rather than wrapping one?
AbstractInterpolation <: AbstractArray imports array-interface obligations (getindex, size, axes). getindex with integer arguments is implemented, but it
just calls the functor — so integer indexing and floating-point evaluation are the same operation. Was the goal to satisfy the Interpolations.jl gradient
extension machinery, or to make CachedInterpolation usable anywhere an AbstractArray is expected? If the latter, the implementation seems incomplete (no
setindex!, no iteration fallbacks). If the former, it might be cleaner to hold an AbstractInterpolation internally and dispatch gradient on the wrapper
rather than subtyping.

Q2. Why is the origin offset O a type parameter rather than a field?
Encoding a runtime value as a type parameter is occasionally done in Julia for performance (constant-folding in generated code), but it also means that
changing the origin requires constructing a new type, and that two otherwise-identical CachedInterpolation objects with different origins have
incompatible types. Was this a deliberate performance choice, and if so, is the gain measurable?

Q3. Should cachedinterpolators accept something other than a bare Array?
The factory signature is cachedinterpolators(parent::Array, N). This excludes memory-mapped arrays wrapped in types that are <:AbstractArray but not
<:Array (e.g., Mmap.MmapArray if it uses a different type). The stated motivation for the package is memory-mapped files — so this annotation may
inadvertently exclude the primary use case on some Julia versions or with some mmap wrappers.

Q4. Hardwired quadratic B-spline — is that a constraint or an oversight?
The caching mechanism is general (any kernel could be cached), but the interpolation kernel is fixed. Was the intent to eventually support other kernels
(cubic, linear), or is "quadratic B-spline" a defining characteristic of the package?

---

Observations

• The inner @noinline cachedinterpolators(buffer, parent, origin, center, sztiles) function-barrier exists purely to work around a type-instability in
  sztiles. This is a known Julia idiom but is entirely invisible to users; it could be documented with a brief comment for future maintainers.
• The getindex method is retained because AbstractArray requires it, and the code comments acknowledge it may be deprecated. If AbstractInterpolation
  inheritance were dropped, getindex could be removed cleanly.
• CachedInterpolation is a mutable struct. Given that center is the only field that changes (cache invalidation), a const annotation on the other fields
  (coefs, parent, tileindex) would communicate intent clearly and is now idiomatic in Julia 1.8+.

---

Overall assessment

CachedInterpolations is a small, focused package with a clear and useful purpose. The core idea — a mutable buffer per tile that short-circuits repeated
disk reads — is sound and well-implemented. The design works for its primary use case (sequential iterative optimization over one tile at a time).

The two highest-leverage changes would be: (1) eliminating the gradient precondition footgun by providing an atomic evaluate_and_gradient function (or by
making gradient trigger its own evaluation when the cache is stale), and (2) documenting and/or restructuring the shared-buffer semantics so that
threading constraints are explicit. Everything else — the type-parameter redundancy, the gradient! allocation issue, the Array vs. AbstractArray
annotation — is either a minor inconsistency or a design question that depends on intent.

---

Phase 5 — Values clarification

Before I write the plan, I'd like to understand your intent on a few tensions the review surfaced. Brief answers are fine.

1. Scope and audience: Is the package specifically for the "one tile at a time, sequential optimizer" use case, or do you want it to be general enough for
   concurrent tile access or arbitrary AbstractArray parents (not just Array)? Is supporting kernels other than quadratic B-spline in scope?
2. Central abstraction: Should CachedInterpolation remain a subtype of AbstractInterpolation (accepting the full AbstractArray obligation that comes with
   it), or would you prefer it to be a standalone type that holds an AbstractInterpolation internally and just extends gradient?
3. Gradient API: The precondition footgun (Update to the new Interpolations #1 in "Likely design issues") could be fixed three ways: (a) add a combined evaluate_and_gradient function; (b)
   make gradient internally call itp(ys...) when the cache is stale (removing the precondition entirely); or (c) leave the precondition but document it more
   prominently. Which direction fits the intended usage pattern?
4. Error and failure model: For out-of-bounds evaluation and the gradient precondition violation — should these throw errors, or is silent wrong-result
   behavior acceptable given the performance-critical context?
5. O as type parameter: Was encoding the origin offset in the type signature a deliberate performance choice you'd like to keep, or something you'd
   revisit?