[ctc] Fix CtcDeriv propagation without explicit gates

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

src/core/contractors/codac2_CtcDeriv.h

@@ -79,27 +79,80 @@ 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);
+
+        // Keep boundary values locally to propagate contractions across slices.
+        // Only gates already present in the TDomain are updated afterwards.
+        std::vector<std::shared_ptr<Slice<T>>> sx;
+        std::vector<std::shared_ptr<const Slice<T>>> sv;
+        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));
+          }
 
-        if((_time_propag & TimePropag::FWD) == TimePropag::FWD)
+        if(sx.empty())
+          return;
+
+        std::vector<T> 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_beg ; it != std::next(it_end) ; it++)
+          T envelope = sx[i]->codomain();
+
+          if constexpr(std::is_same_v<T,Interval>)
+            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<T,IntervalVector>)
           {
-            auto sx = x.slice(it);
-            if(!sx->is_gate())
-              this->contract(*sx, *v.slice(it), ctc_indices);
+            std::vector<Index> 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 gates explicitly present in the TDomain.
+        // Temporary boundary values are discarded after contraction.
+        for(std::size_t i = 0 ; i < sx.size() ; i++)
         {
-          for(auto it = it_end ; it != std::prev(it_beg) ; it--)
-          {
-            auto sx = x.slice(it);
-            if(!sx->is_gate())
-              this->contract(*sx, *v.slice(it), ctc_indices);
-          }
+          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);
         }
       }
 

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<typename OtherDerived>
   inline std::pair<IntervalMatrix,Interval>

tests/core/contractors/codac2_tests_CtcDeriv.cpp

@@ -12,6 +12,7 @@
 #include <codac2_SlicedTube.h>
 #include <codac2_Approx.h>
 #include <iomanip>
+#include <iterator>
 
 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<Interval> x(tdomain, Interval());
+      SlicedTube<Interval> 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<Interval> x(tdomain, Interval());
+    SlicedTube<Interval> 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<Interval> x(tdomain, Interval());
+    SlicedTube<Interval> 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.);

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.)