From e0b616bcd9a4f995e80a98e5c65e2dfd3bd8b601 Mon Sep 17 00:00:00 2001 From: Simon Rohou Date: Wed, 7 Oct 2026 08:52:11 +0200 Subject: [PATCH 1/3] [doc] updated inversion.rst --- doc/manual/manual/linear/inversion.rst | 3 +++ src/core/matrices/codac2_inversion.h | 6 ++---- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/manual/manual/linear/inversion.rst b/doc/manual/manual/linear/inversion.rst index 5d9f024bb..be4afb802 100644 --- a/doc/manual/manual/linear/inversion.rst +++ b/doc/manual/manual/linear/inversion.rst @@ -39,4 +39,7 @@ Reliable inversions of matrices :project: codac .. doxygenfunction:: codac2::infinite_sum_enclosure + :project: codac + +.. doxygenfunction:: codac2::cofactor_matrix_enclosure :project: codac \ No newline at end of file diff --git a/src/core/matrices/codac2_inversion.h b/src/core/matrices/codac2_inversion.h index 94f3df0a9..216871ec7 100644 --- a/src/core/matrices/codac2_inversion.h +++ b/src/core/matrices/codac2_inversion.h @@ -119,11 +119,9 @@ namespace codac2 /** - * \brief compute the interval cofactor matrix, as well as the - * determinant, of a square matrix + * \brief Compute the interval cofactor matrix, as well as the determinant, of a square matrix * \param A A matrix expression, possibly interval. - * \return The enclosure of the cofactor matrix (first) - * and determinant (second) + * \return The enclosure of the cofactor matrix (first) and determinant (second) */ template inline std::pair From 66632002f55a4eb398c27e5f153b25840c40fd9b Mon Sep 17 00:00:00 2001 From: Simon Rohou Date: Thu, 8 Oct 2026 21:44:57 +0200 Subject: [PATCH 2/3] [ctc] Fix CtcDeriv propagation without explicit gates --- src/core/contractors/codac2_CtcDeriv.h | 72 +++++++++++++++++++++----- 1 file changed, 60 insertions(+), 12 deletions(-) diff --git a/src/core/contractors/codac2_CtcDeriv.h b/src/core/contractors/codac2_CtcDeriv.h index a53b557eb..01393bb56 100644 --- a/src/core/contractors/codac2_CtcDeriv.h +++ b/src/core/contractors/codac2_CtcDeriv.h @@ -82,24 +82,72 @@ namespace codac2 auto it_beg = x.tdomain()->sample(t.lb(),true); auto it_end = x.tdomain()->sample(t.ub(),true); - if((_time_propag & TimePropag::FWD) == TimePropag::FWD) - { - for(auto it = it_beg ; it != std::next(it_end) ; it++) + // Store gates locally during contraction to preserve boundary information + // and ensure efficient propagation between slices. + // Only explicitly defined gates are retained in the TDomain. + // Temporary gates are discarded after contraction. + std::vector>> sx; + std::vector>> sv; + for(auto it = it_beg ; it != std::next(it_end) ; it++) + if(!it->is_gate()) { - auto sx = x.slice(it); - if(!sx->is_gate()) - this->contract(*sx, *v.slice(it), ctc_indices); + sx.push_back(x.slice(it)); + sv.push_back(v.slice(it)); } - } - if((_time_propag & TimePropag::BWD) == TimePropag::BWD) + if(sx.empty()) + return; + + std::vector gates; + gates.reserve(sx.size()+1); + gates.push_back(sx.front()->input_gate()); + for(const auto& s : sx) + gates.push_back(s->output_gate()); + + auto contract_slice = [&](std::size_t i) { - for(auto it = it_end ; it != std::prev(it_beg) ; it--) + T envelope = sx[i]->codomain(); + + if constexpr(std::is_same_v) + contract(sx[i]->t0_tf(), envelope, gates[i], gates[i+1], + sv[i]->codomain(), _time_propag, _fast_mode); + + else if constexpr(std::is_same_v) { - auto sx = x.slice(it); - if(!sx->is_gate()) - this->contract(*sx, *v.slice(it), ctc_indices); + std::vector ctc_indices_(ctc_indices); + if(ctc_indices_.empty()) + { + ctc_indices_.resize(sx[i]->size()); + std::iota(ctc_indices_.begin(), ctc_indices_.end(), 0); + } + + for(auto j : ctc_indices_) + contract(sx[i]->t0_tf(), envelope[j], gates[i][j], gates[i+1][j], + sv[i]->codomain()[j], _time_propag, _fast_mode); } + + sx[i]->set(envelope, false); + }; + + if((_time_propag & TimePropag::FWD) == TimePropag::FWD) + for(std::size_t i = 0 ; i < sx.size() ; i++) + contract_slice(i); + + if((_time_propag & TimePropag::BWD) == TimePropag::BWD) + for(std::size_t i = sx.size() ; i-- > 0 ; ) + contract_slice(i); + + // Update only the gates explicitly defined in the TDomain. + // Temporary gates are discarded after contraction. + for(std::size_t i = 0 ; i < sx.size() ; i++) + { + auto prev = sx[i]->prev_slice(); + if(prev && prev->is_gate()) + prev->set(gates[i], false); + + auto next = sx[i]->next_slice(); + if(next && next->is_gate()) + next->set(gates[i+1], false); } } From df2dd36bf449df8089834bd14a983b3957f3cbe4 Mon Sep 17 00:00:00 2001 From: Simon Rohou Date: Thu, 8 Oct 2026 22:30:19 +0200 Subject: [PATCH 3/3] [ctc] Add CtcDeriv regression tests for propagation without explicit gates --- src/core/contractors/codac2_CtcDeriv.h | 27 +++--- .../contractors/codac2_tests_CtcDeriv.cpp | 82 +++++++++++++++++++ .../core/contractors/codac2_tests_CtcDeriv.py | 60 ++++++++++++++ 3 files changed, 158 insertions(+), 11 deletions(-) diff --git a/src/core/contractors/codac2_CtcDeriv.h b/src/core/contractors/codac2_CtcDeriv.h index 01393bb56..3202e650e 100644 --- a/src/core/contractors/codac2_CtcDeriv.h +++ b/src/core/contractors/codac2_CtcDeriv.h @@ -79,17 +79,22 @@ namespace codac2 { assert_release(TDomain::are_same(x.tdomain(), v.tdomain())); Interval t = x.tdomain()->t0_tf() & _tdomain; - auto it_beg = x.tdomain()->sample(t.lb(),true); - auto it_end = x.tdomain()->sample(t.ub(),true); + if(t.is_empty() || t.is_degenerated()) + return; + + // Refine only when the restricted domain cuts through a slice. + // No explicit gates are created. + if(t.lb() > x.tdomain()->t0_tf().lb()) + x.tdomain()->sample(t.lb(), false); + if(t.ub() < x.tdomain()->t0_tf().ub()) + x.tdomain()->sample(t.ub(), false); - // Store gates locally during contraction to preserve boundary information - // and ensure efficient propagation between slices. - // Only explicitly defined gates are retained in the TDomain. - // Temporary gates are discarded after contraction. + // Keep boundary values locally to propagate contractions across slices. + // Only gates already present in the TDomain are updated afterwards. std::vector>> sx; std::vector>> sv; - for(auto it = it_beg ; it != std::next(it_end) ; it++) - if(!it->is_gate()) + for(auto it = x.tdomain()->begin(); it != x.tdomain()->end(); ++it) + if(!it->is_gate() && it->is_subset(t)) { sx.push_back(x.slice(it)); sv.push_back(v.slice(it)); @@ -136,9 +141,9 @@ namespace codac2 if((_time_propag & TimePropag::BWD) == TimePropag::BWD) for(std::size_t i = sx.size() ; i-- > 0 ; ) contract_slice(i); - - // Update only the gates explicitly defined in the TDomain. - // Temporary gates are discarded after contraction. + + // Update only gates explicitly present in the TDomain. + // Temporary boundary values are discarded after contraction. for(std::size_t i = 0 ; i < sx.size() ; i++) { auto prev = sx[i]->prev_slice(); diff --git a/tests/core/contractors/codac2_tests_CtcDeriv.cpp b/tests/core/contractors/codac2_tests_CtcDeriv.cpp index 11683a5b0..07eb5a27a 100644 --- a/tests/core/contractors/codac2_tests_CtcDeriv.cpp +++ b/tests/core/contractors/codac2_tests_CtcDeriv.cpp @@ -12,6 +12,7 @@ #include #include #include +#include using namespace std; using namespace codac2; @@ -226,6 +227,87 @@ TEST_CASE("CtcDeriv") CHECK(sv->codomain() == Interval(-1,oo)); } + SECTION("Propagation without internal gates") + { + for(auto direction : {TimePropag::FWD, TimePropag::BWD, TimePropag::FWD_BWD}) + { + auto tdomain = create_tdomain({0,3}, 1., false); + SlicedTube x(tdomain, Interval()); + SlicedTube v(tdomain, Interval(1)); + + auto it_first = tdomain->begin(); + auto it_middle = std::next(it_first); + auto it_last = std::prev(tdomain->end()); + + // Constrain whole slices; setting a value at an instant would add a gate. + if(direction != TimePropag::BWD) + x.slice(it_first)->set(Interval(0,1)); + if(direction != TimePropag::FWD) + x.slice(it_last)->set(Interval(2,3)); + + REQUIRE(tdomain->nb_tslices() == 3); + CHECK_FALSE(tdomain->all_gates_defined()); + + CtcDeriv(direction, false).contract(x, v); + + if(direction == TimePropag::FWD) + { + CHECK(x.slice(it_first)->codomain() == Interval(0,1)); + CHECK(x.slice(it_middle)->codomain().is_subset(Interval(0,2))); + CHECK(x.slice(it_last)->codomain().is_subset(Interval(1,3))); + } + else if(direction == TimePropag::BWD) + { + CHECK(x.slice(it_first)->codomain().is_subset(Interval(-1,2))); + CHECK(x.slice(it_middle)->codomain().is_subset(Interval(0,3))); + CHECK(x.slice(it_last)->codomain() == Interval(2,3)); + } + else + { + CHECK(x.slice(it_first)->codomain() == Interval(0,1)); + CHECK(x.slice(it_middle)->codomain().is_subset(Interval(0,3))); + CHECK(x.slice(it_last)->codomain() == Interval(2,3)); + } + + // Contracting a tube must not create any explicit gate. + CHECK(tdomain->nb_tslices() == 3); + CHECK_FALSE(tdomain->all_gates_defined()); + for(const auto& ts : *tdomain) + CHECK_FALSE(ts.is_gate()); + } + } + + SECTION("No new gates in an existing partially gated TDomain") + { + auto tdomain = create_tdomain({0,3}, 1., false); + SlicedTube x(tdomain, Interval()); + SlicedTube v(tdomain, Interval(1)); + x.set(Interval(0), 0.); // Existing gate must be preserved. + + REQUIRE(tdomain->nb_tslices() == 4); + CtcDeriv(TimePropag::FWD_BWD, false).contract(x, v); + CHECK(tdomain->nb_tslices() == 4); + CHECK(x(0.) == Interval(0)); + for(const auto& ts : *tdomain) + if(ts.is_gate()) + CHECK(ts.lb() == 0.); + } + + SECTION("No new gates with a restricted temporal domain") + { + auto tdomain = create_tdomain({0,4}, 1., false); + SlicedTube x(tdomain, Interval()); + SlicedTube v(tdomain, Interval(1)); + CtcDeriv ctc(TimePropag::FWD_BWD, false); + ctc.restrict_tdomain(Interval(1,3)); + + REQUIRE(tdomain->nb_tslices() == 4); + ctc.contract(x, v); + CHECK(tdomain->nb_tslices() == 4); + for(const auto& ts : *tdomain) + CHECK_FALSE(ts.is_gate()); + } + SECTION("Test fwd") { auto tdomain = create_tdomain({0,6},1.); diff --git a/tests/core/contractors/codac2_tests_CtcDeriv.py b/tests/core/contractors/codac2_tests_CtcDeriv.py index aa9cf9337..57c55afb3 100644 --- a/tests/core/contractors/codac2_tests_CtcDeriv.py +++ b/tests/core/contractors/codac2_tests_CtcDeriv.py @@ -213,6 +213,66 @@ def test_CtcDeriv(self): self.assertTrue(sx.codomain() == Interval(-5,4)) self.assertTrue(sv.codomain() == Interval(-1,oo)) + # Propagation without internal gates + + for direction in (TimePropag.FWD, TimePropag.BWD, TimePropag.FWD_BWD): + tdomain = create_tdomain([0,3], 1., False) + x = SlicedTube(tdomain, Interval()) + v = SlicedTube(tdomain, Interval(1)) + sx = [s for s in x if not s.is_gate()] + + # Constrain whole slices, without creating pointwise gates. + if direction != TimePropag.BWD: + sx[0].set(Interval(0,1)) + if direction != TimePropag.FWD: + sx[-1].set(Interval(2,3)) + + self.assertEqual(len(sx), 3) + self.assertFalse(tdomain.all_gates_defined()) + + CtcDeriv(direction, False).contract(x, v) + + if direction == TimePropag.FWD: + self.assertEqual(sx[0].codomain(), Interval(0,1)) + self.assertTrue(sx[1].codomain().is_subset(Interval(0,2))) + self.assertTrue(sx[2].codomain().is_subset(Interval(1,3))) + elif direction == TimePropag.BWD: + self.assertTrue(sx[0].codomain().is_subset(Interval(-1,2))) + self.assertTrue(sx[1].codomain().is_subset(Interval(0,3))) + self.assertEqual(sx[2].codomain(), Interval(2,3)) + else: + self.assertEqual(sx[0].codomain(), Interval(0,1)) + self.assertTrue(sx[1].codomain().is_subset(Interval(0,3))) + self.assertEqual(sx[2].codomain(), Interval(2,3)) + + # Contracting a tube must not create any explicit gate. + slices = list(x) + self.assertEqual(len(slices), 3) + self.assertFalse(tdomain.all_gates_defined()) + self.assertTrue(all(not s.is_gate() for s in slices)) + + # Existing gates must be preserved without creating new ones. + tdomain = create_tdomain([0,3], 1., False) + x = SlicedTube(tdomain, Interval()) + v = SlicedTube(tdomain, Interval(1)) + x.set(Interval(0), 0.) + self.assertEqual(len(list(x)), 4) + CtcDeriv(TimePropag.FWD_BWD, False).contract(x, v) + self.assertEqual(len(list(x)), 4) + self.assertEqual(x(0.), Interval(0)) + self.assertEqual(len([s for s in x if s.is_gate()]), 1) + + # Restricted contraction at existing slice boundaries must not create gates. + tdomain = create_tdomain([0,4], 1., False) + x = SlicedTube(tdomain, Interval()) + v = SlicedTube(tdomain, Interval(1)) + ctc = CtcDeriv(TimePropag.FWD_BWD, False) + ctc.restrict_tdomain(Interval(1,3)) + self.assertEqual(len(list(x)), 4) + ctc.contract(x, v) + self.assertEqual(len(list(x)), 4) + self.assertTrue(all(not s.is_gate() for s in x)) + # Test fwd tdomain = create_tdomain([0,6],1.)