feat(model): Taylor models with rigorous remainder bounds (Makino ch. 4) - #31
Open
andreapasquale94 wants to merge 5 commits into
Open
feat(model): Taylor models with rigorous remainder bounds (Makino ch. 4)#31andreapasquale94 wants to merge 5 commits into
andreapasquale94 wants to merge 5 commits into
Conversation
Add a new header-only module tax::model implementing remainder-enhanced differential algebra (Taylor models) after K. Makino's PhD thesis (MSUCL-1093, 1998), chapter 4: - Interval<T>: outward-rounded interval arithmetic (constexpr bit_cast nextUp/nextDown, 1 ulp per arithmetic op, 2 ulps on libm endpoints), sharp even-power rule (5.4), enclosures for exp/log/sqrt/sin/cos/ sinh/cosh with extremum detection for the trigonometric ones. - TaylorModel<T, N, M> (alias TM<N, M>): dense TaylorExpansion polynomial part + remainder interval + expansion point + domain box, guaranteeing f(x) in P(x - x0) + I on the domain. Range bounds via per-order interval sums; addition/subtraction per (4.5); multiplication folds the degree->N excess of the product into the remainder (5.3.2); antiderivation per (4.12); interval-scalar operands supported. - Intrinsics per the 4.3.2 recipe (constant-part split, Horner series in TM arithmetic, Lagrange remainder over c + hull(0, B(Pbar) + I)): exp, log, sqrt, isqrt, reciprocal, square, pow, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, and division. Domain violations throw std::domain_error. asin/acos/atan use the direct Taylor form with the thesis derivative recursions instead of the addition formulas (which need branch conditions the thesis leaves implicit). - Tests: interval unit tests, TM arithmetic/intrinsic containment sweeps, antiderivation, and thesis worked examples kept as regression oracles (reciprocal model (4.15)/(4.16), Table 4.2 remainder values, Table 4.3 interval-vs-TM comparison, the 4.4.2 enclosure functions, and the 5.5.2 verified double integral). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EjvS2WURjgthURqEtJimmS
- guide/models.md: task-oriented walkthrough (why not plain intervals, factories, arithmetic/intrinsics, bounds, verified integration). - reference/models.md: complete Interval + TaylorModel API surface, rounding contract, domain conditions, exception summary. - internals/taylor-models.md: full implementation details — outward rounding mechanics, sin/cos extremum detection, range bounding, multiplication excess algorithm, per-intrinsic coefficient recursions and Lagrange enclosures, inverse-trig deviation from the thesis, antiderivation, rigor contract, complexity, and an assessment of what thesis ch. 5 adds (per-order bound caching, sharper bounders, coefficient-error sweeping) with a suggested implementation order. - Wire the pages into mkdocs nav and section index tables. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EjvS2WURjgthURqEtJimmS
…olicy
Add tax::model::Bounder { Naive, Quadratic } selecting how TaylorModel's
polynomial part is bounded, defaulting to the sharper Quadratic strategy
for bound() / polynomialBound(). Implements the diagonal case of the
thesis's exact quadratic bounder (§5.4.3):
- boundDiagonalQuadratic: exact 1-D enclosure of q_ii*h_i^2 + g_i*h_i over
the domain interval via monotone endpoint hull plus a folded-in vertex
value when the parabola vertex may lie inside; recovers interior minima
the naive order-sum misses (e.g. (h-0.3)^2 on [-1,1]: naive [-0.51,1.69],
quadratic exact [0,1.69]).
- quadraticRangeBound: sums the exact per-variable diagonal bounds, keeps
cross terms (exact via independent interval products) and orders >= 3 on
the naive sum, and intersects with the naive bound so the result is a
valid enclosure that is never wider than naive. Drops the order-3 1/x+x
total bound from 0.15138 to 0.15013 (toward the thesis 0.14988).
Range-bound machinery (DomainPowers, polyRangeBound, orderRangeBound,
excessProductBound) moves into the new model/bounders.hpp; remainder
propagation deliberately stays on the naive sum (as the thesis does), so
the ch. 4 remainder oracles remain pinned. New tests in test_bounders.cpp
cover interior/exterior vertices, tightening vs naive, enclosure validity,
multivariate cross terms, and low-order agreement. Docs (guide, reference,
internals, CLAUDE.md) updated to describe the strategy and the remaining
full-multivariate upgrade.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EjvS2WURjgthURqEtJimmS
…e helpers Add the Taylor-model surface needed to build a verified ODE integrator on top (Makino-Berz Picard iteration r = r0 + integral F(r) dt over TMs in the initial-condition and time variables). - fix<I>(v) / fix(i,v): partial evaluation — evaluate one variable at a scalar coordinate exactly (a point adds no remainder), collapsing that axis while the others stay symbolic. The step-continuation primitive: fixing the time axis at the step endpoint yields the end-of-step state as a function of the initial conditions. - retarget<I>(x0,dom): reset a collapsed variable's expansion point/domain (guarded to models independent of that variable) so the time slot can be recycled for the next step's [0,h] window. - compose(G,H) (model/compose.hpp): substitute a TM-vector into a TM (and a vector of outer models), via multivariate Horner in TM arithmetic + I_G, with a rigorous inner-range subset outer-domain check. Flow/map composition for long-time integration and Poincare maps. - model/io.hpp: operator<< / to_string for TaylorModel; value/bound/jacobian over a std::array<TM,D> state (jacobian = state-transition matrix from the polynomial parts). No Eigen NumTraits<TM> — a domain-carrying TM(0) literal is incompatible with real-domain TMs, so state vectors use std::array. Only the antiderivation (Picard) direction is needed; the derivation operator remains intentionally absent (thesis §4.3.3). Tests: fix/retarget/ compose/io units (test_ode_primitives.cpp) and an end-to-end harmonic- oscillator integrator proving the pipeline — flow matches the analytic Taylor series exactly, encloses the true solution across the (IC, time) box, and multi-step fix/retarget continuation recovers cos/sin and the STM (test_ode_integration.cpp). Docs (guide 'Toward an ODE integrator', reference, internals, CLAUDE.md, README) updated. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EjvS2WURjgthURqEtJimmS
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #31 +/- ##
==========================================
- Coverage 96.57% 96.51% -0.06%
==========================================
Files 80 98 +18
Lines 3444 4739 +1295
Branches 549 803 +254
==========================================
+ Hits 3326 4574 +1248
- Misses 47 61 +14
- Partials 71 104 +33 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Make TaylorModel a first-class Eigen scalar so ODE state vectors and state-transition matrices can be Eigen::Matrix<TM, ...> instead of hand-rolled std::array. The obstacle is that Eigen synthesises Scalar(0)/Scalar(1) literals (in setZero, Identity, reductions, products) with no domain, which a strict domain-carrying TM would reject when combined with a real-domain model. The fix is mathematically principled rather than a special case: - A TM built from a bare scalar (new implicit TaylorModel(T) ctor) is a *domain-agnostic constant* (abstract_ flag): a constant function is valid over every domain, so it adopts its partner's expansion point/domain in any binary operation. isAbstractConstant()/overDomain()/asAbstractConstant() expose and re-home it. - detail::reconcile promotes an abstract operand onto the concrete partner's domain; detail::keepAbstract threads the flag through the scalar/interval operators so a literal-derived constant never decays into a concrete degenerate-domain landmine. The concrete-concrete fast path is unchanged (guarded by two bool reads) and still throws on mismatched domains. - model/eigen.hpp: Eigen::NumTraits<TaylorModel> specialization, TMVec/TMMat aliases, an Eigen-returning variables() factory, and Eigen overloads of value/bound/jacobian (state-transition matrix). Supports Matrix<TM> storage, +/-/*, matrix and matrix-vector products (the variational STM * state pattern), transpose, Identity/Zero. Unordered, so Eigen's fuzzy isApprox and pivoted decompositions are unsupported (documented). Tests (test_eigen.cpp): abstract-constant reconciliation + rigor, matrix- vector products with literals, Identity/Zero, matrix-matrix products, the vector helpers, and a variational harmonic-oscillator step. Docs (guide, reference, internals, CLAUDE.md, README) updated; the earlier 'no NumTraits' notes are corrected. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EjvS2WURjgthURqEtJimmS
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Add a new header-only module tax::model implementing remainder-enhanced
differential algebra (Taylor models) after K. Makino's PhD thesis
(MSUCL-1093, 1998), chapter 4:
nextUp/nextDown, 1 ulp per arithmetic op, 2 ulps on libm endpoints),
sharp even-power rule (5.4), enclosures for exp/log/sqrt/sin/cos/
sinh/cosh with extremum detection for the trigonometric ones.
part + remainder interval + expansion point + domain box, guaranteeing
f(x) in P(x - x0) + I on the domain. Range bounds via per-order interval
sums; addition/subtraction per (4.5); multiplication folds the
degree->N excess of the product into the remainder (5.3.2);
antiderivation per (4.12); interval-scalar operands supported.
TM arithmetic, Lagrange remainder over c + hull(0, B(Pbar) + I)):
exp, log, sqrt, isqrt, reciprocal, square, pow, sin, cos, tan, asin,
acos, atan, sinh, cosh, tanh, and division. Domain violations throw
std::domain_error. asin/acos/atan use the direct Taylor form with the
thesis derivative recursions instead of the addition formulas (which
need branch conditions the thesis leaves implicit).
antiderivation, and thesis worked examples kept as regression oracles
(reciprocal model (4.15)/(4.16), Table 4.2 remainder values, Table 4.3
interval-vs-TM comparison, the 4.4.2 enclosure functions, and the
5.5.2 verified double integral).
Co-Authored-By: Claude Fable 5 noreply@anthropic.com
Claude-Session: https://claude.ai/code/session_01EjvS2WURjgthURqEtJimmS