From 58573fd01dc9f86b41f535cdcfcbd31f36a551e7 Mon Sep 17 00:00:00 2001
From: SimonRohou <simon.rohou@gmail.com>
Date: Mon, 4 Dec 2023 18:07:48 +0100
Subject: [PATCH] [codac2] templated paving with weak_ptrs

---
 src/core/2/cn/codac2_ContractorNetwork.cpp |   8 +-
 src/core/2/cn/codac2_ContractorNetwork.h   |   2 +-
 src/core/2/cn/codac2_ContractorNode.h      |   6 +-
 src/core/2/cn/codac2_DomainNode.h          |   6 +
 src/core/2/contractors/codac2_CtcInter.h   |  11 +-
 src/core/2/contractors/codac2_CtcPaver.h   |  13 +-
 src/core/2/domains/paving/codac2_Paving.h  | 219 ++++++++++++---------
 7 files changed, 151 insertions(+), 114 deletions(-)

diff --git a/src/core/2/cn/codac2_ContractorNetwork.cpp b/src/core/2/cn/codac2_ContractorNetwork.cpp
index 052b17081..b2d67a0c8 100644
--- a/src/core/2/cn/codac2_ContractorNetwork.cpp
+++ b/src/core/2/cn/codac2_ContractorNetwork.cpp
@@ -29,7 +29,7 @@ namespace codac2
     _auto_fixpoint = !disable;
   }
 
-  void ContractorNetwork::contract(bool verbose)
+  double ContractorNetwork::contract(bool verbose)
   {
     if(verbose)
     {
@@ -45,7 +45,7 @@ namespace codac2
       shared_ptr<ContractorNodeBase> current_ctc = _stack.front();
       _stack.pop_front();
 
-      auto contracted_doms = current_ctc->call_contract();
+      auto contracted_doms = current_ctc->call_contract(false);
 
       for(auto& d : contracted_doms)
         for(auto& ci : d->contractors())
@@ -56,8 +56,10 @@ namespace codac2
         }
     }
     
+    double elapsed_time = (double)(clock()-t_start)/CLOCKS_PER_SEC;
     if(verbose)
-      std::cout << "  Constraint propagation time: " << (double)(clock() - t_start)/CLOCKS_PER_SEC << "s" << std::endl;
+      std::cout << "  Constraint propagation time: " << elapsed_time << "s" << std::endl;
+    return elapsed_time;
   }
   
   void ContractorNetwork::reset_all_vars()
diff --git a/src/core/2/cn/codac2_ContractorNetwork.h b/src/core/2/cn/codac2_ContractorNetwork.h
index 2c4256213..d7a0ab2d6 100644
--- a/src/core/2/cn/codac2_ContractorNetwork.h
+++ b/src/core/2/cn/codac2_ContractorNetwork.h
@@ -18,7 +18,7 @@ namespace codac2
       void add(const std::shared_ptr<ContractorNodeBase>& ctc);
       void add_ctc_to_stack(const std::shared_ptr<ContractorNodeBase>& ctc);
       void disable_auto_fixpoint(bool disable = true);
-      void contract(bool verbose = true);
+      double contract(bool verbose = true);
 
       void reset_all_vars();
 
diff --git a/src/core/2/cn/codac2_ContractorNode.h b/src/core/2/cn/codac2_ContractorNode.h
index 9e0ede9cd..2d2d379ff 100644
--- a/src/core/2/cn/codac2_ContractorNode.h
+++ b/src/core/2/cn/codac2_ContractorNode.h
@@ -18,7 +18,7 @@ namespace codac2
 
       virtual ~ContractorNodeBase() = default;
       virtual size_t nb_args() const = 0;
-      virtual std::list<std::shared_ptr<DomainNodeBase>> call_contract() = 0;
+      virtual std::list<std::shared_ptr<DomainNodeBase>> call_contract(bool verbose = false) = 0;
       virtual Contractor* raw_ptr() const = 0;
       virtual void associate_domains(std::vector<std::shared_ptr<DomainNodeBase>>& cn_domains) = 0;
       virtual std::string contractor_class_name() const = 0;
@@ -60,8 +60,10 @@ namespace codac2
           l.push_back(d);
       }
 
-      std::list<std::shared_ptr<DomainNodeBase>> call_contract()
+      std::list<std::shared_ptr<DomainNodeBase>> call_contract(bool verbose = false)
       {
+        if(verbose)
+          std::cout << "Calling " << contractor_class_name() << std::endl;
         std::list<std::shared_ptr<DomainNodeBase>> contracted_doms;
 
         std::apply(
diff --git a/src/core/2/cn/codac2_DomainNode.h b/src/core/2/cn/codac2_DomainNode.h
index 2673166ba..a9c38f25c 100644
--- a/src/core/2/cn/codac2_DomainNode.h
+++ b/src/core/2/cn/codac2_DomainNode.h
@@ -103,6 +103,12 @@ namespace codac2
           _x.get().reset();
       }
 
+      void reset(const T& x)
+      {
+        if constexpr(_is_var)
+          _x.get().reset(x);
+      }
+
       constexpr bool is_var() const
       {
         return _is_var;
diff --git a/src/core/2/contractors/codac2_CtcInter.h b/src/core/2/contractors/codac2_CtcInter.h
index 667cd7a6f..0fec6b43b 100644
--- a/src/core/2/contractors/codac2_CtcInter.h
+++ b/src/core/2/contractors/codac2_CtcInter.h
@@ -28,13 +28,10 @@ namespace codac2
 
       void contract(Paving<N>& x, IntervalVector_<N>& y)
       {
-        IntervalVector_<N> y_ = IntervalVector_<N>::empty_set();
-        for(auto& li : x.leaves_list())
-        {
-          y_ |= li->_x & y;
-          //li->_x &= y;
-        }
-        y = y_;
+        auto y_ = y;
+        y.set_empty();
+        for(const auto& l : x.leaves())
+          y |= y_ & l.lock()->box();
       }
 
       make_available_in_cn__templated(CtcInter<N>)
diff --git a/src/core/2/contractors/codac2_CtcPaver.h b/src/core/2/contractors/codac2_CtcPaver.h
index 59d0beb96..9d8318e8f 100644
--- a/src/core/2/contractors/codac2_CtcPaver.h
+++ b/src/core/2/contractors/codac2_CtcPaver.h
@@ -21,9 +21,9 @@ namespace codac2
   {
     public:
 
-      static bool compare_paving(Paving<N>* first, Paving<N>* second)
+      static bool compare_paving(const std::weak_ptr<PavingNode<N>>& first, const std::weak_ptr<PavingNode<N>>& second)
       {
-        return first->box().volume() > second->box().volume();
+        return first.lock()->_x.volume() > second.lock()->_x.volume();
       }
 
       CtcPaver(ContractorOnBox<N>& ctc, size_t max_leaves)
@@ -34,15 +34,16 @@ namespace codac2
 
       void contract(Paving<N>& x)
       {
-        list<Paving<N>*> l_subpav = x.leaves_list();
+        auto l_subpav = x.leaves();
         size_t nb_leaves = l_subpav.size();
 
         while(!l_subpav.empty())
         {
           l_subpav.sort(compare_paving);
 
-          Paving<N>* q = l_subpav.front();
+          std::shared_ptr<PavingNode<N>> q = l_subpav.front().lock();
           l_subpav.pop_front();
+          assert(q.use_count() > 1);
 
           _ctc.contract(q->_x);
 
@@ -53,8 +54,8 @@ namespace codac2
           {
             q->bisect();
             nb_leaves++;
-            l_subpav.push_back(q->_left.get());
-            l_subpav.push_back(q->_right.get());
+            l_subpav.push_back(q->left());
+            l_subpav.push_back(q->right());
           }
         }
       }
diff --git a/src/core/2/domains/paving/codac2_Paving.h b/src/core/2/domains/paving/codac2_Paving.h
index ef2d873bf..35ab27b89 100644
--- a/src/core/2/domains/paving/codac2_Paving.h
+++ b/src/core/2/domains/paving/codac2_Paving.h
@@ -21,30 +21,31 @@
 
 namespace codac2
 {
-  template<class P,int N=Dynamic>
-  class PavingBase : virtual public Domain
+  template<int N=Dynamic>
+  class Paving;
+
+  template<int N=Dynamic>
+  class PavingNode : public std::enable_shared_from_this<PavingNode<N>>
   {
     public:
 
-      explicit PavingBase(const IntervalVector_<N>& x)
-       : _x(x)
-      {
-
-      }
+      explicit PavingNode(Paving<N>& paving, const IntervalVector_<N>& x)
+       : _paving(paving), _x(x)
+      { }
 
-      virtual ~PavingBase() = default;
+      virtual ~PavingNode() = default;
 
       const IntervalVector_<N>& box() const
       {
         return _x;
       }
 
-      std::shared_ptr<P> left()
+      std::shared_ptr<PavingNode<N>> left()
       {
         return _left;
       }
 
-      std::shared_ptr<P> right()
+      std::shared_ptr<PavingNode<N>> right()
       {
         return _right;
       }
@@ -63,131 +64,159 @@ namespace codac2
         return not _left && not _right;
       }
 
-      double volume() const
-      {
-        if(is_leaf())
-          return _x.volume();
-        double v = 0.;
-        if(_left) v += _left->volume();
-        if(_right) v += _right->volume();
-        return v;
-      }
-
-      virtual void bisect(float ratio = 0.49)
+      void bisect(float ratio = 0.49)
       {
         assert(Interval(0.,1.).interior_contains(ratio));
         assert(is_leaf() && "only leaves can be bisected");
         assert(_x.is_bisectable());
         auto p = _x.bisect(ratio);
-        _left = std::make_shared<P>(p.first);
-        _right = std::make_shared<P>(p.second);
+
+        _left = std::make_shared<PavingNode<N>>(_paving, p.first);
+        _right = std::make_shared<PavingNode<N>>(_paving, p.second);
+
+        auto wp = this->weak_from_this();
+        _paving._leaves.remove_if([wp](std::weak_ptr<PavingNode<N>> p)
+        {
+          auto swp = wp.lock(), sp = p.lock();
+          if(swp && sp) return swp == sp;
+          return false;
+        });
+
+        _paving._leaves.push_back(_left->weak_from_this());
+        _paving._leaves.push_back(_right->weak_from_this());
       }
 
-      IntervalVector_<N> hull_box() const
+    protected:
+
+      template<int N_>
+      friend class Paving;
+
+      void recurse__boxes_list(std::list<std::reference_wrapper<const IntervalVector_<N>>>& l, const IntervalVector_<N>& intersect = IntervalVector_<N>()) const
       {
-        if(is_leaf())
-          return _x;
-        auto hull = IntervalVector_<N>::empty_set();
-        if(_left) hull |= _left->hull_box();
-        if(_right) hull |= _right->hull_box();
-        return hull;
+        if(is_leaf() && !_x.is_empty() && _x.intersects(intersect))
+          l.push_back(std::cref(_x));
+        else
+        {
+          if(_left) _left->recurse__boxes_list(l);
+          if(_right) _right->recurse__boxes_list(l);
+        }
       }
 
-      std::list<std::reference_wrapper<const IntervalVector_<N>>> boxes_list(const IntervalVector_<N>& intersect = IntervalVector_<N>()) const
+    public: // todo
+
+      IntervalVector_<N> _x;
+      //std::weak_ptr<PavingNode<N>> _parent = nullptr;
+      std::shared_ptr<PavingNode<N>> _left = nullptr, _right = nullptr;
+      Paving<N>& _paving;
+  };
+
+  template<int N>
+  class Paving : public Domain
+  {
+    public:
+
+      explicit Paving(size_t n)
+       : Paving<N>(IntervalVector_<N>(n))
+      { }
+
+      explicit Paving(const IntervalVector_<N>& x)
+       : _tree(std::make_shared<PavingNode<N>>(*this, IntervalVector_<N>(x)))
       {
-        std::list<std::reference_wrapper<const IntervalVector_<N>>> l;
-        boxes_list_push(l, intersect);
-        return l;
+        _leaves.push_back(_tree->weak_from_this());
       }
 
-      std::list<P*> leaves_list()
+      Paving(const Paving<N>&) = delete;
+      Paving& operator=(const Paving<N>&) = delete;
+
+      size_t nb_leaves() const
       {
-        std::list<P*> l;
-        leaves_list_push(l);
-        return l;
+        return _leaves.size();
       }
 
-      size_t nb_leaves() const
+      IntervalVector_<N> hull_box() const
       {
-        if(is_leaf())
-          return 1;
-        else
-          return (_left ? _left->nb_leaves() : 0) +
-            (_right ? _right->nb_leaves() : 0);
+        IntervalVector_<N> hull = IntervalVector_<N>::empty_set();
+        for(const auto& l : _leaves)
+          hull |= l.lock()->box();
+        return hull;
+      }
+
+      bool is_empty() const
+      {
+        for(const auto& l : _leaves)
+          if(!l.lock()->box()->is_empty())
+            return false;
+        return true;
+      }
+
+      double volume() const
+      {
+        double vol = 0.;
+        for(const auto& l : _leaves)
+          vol += l.lock()->box().volume();
+        return vol;
       }
       
       DomainVolume dom_volume() const
       {
         DomainVolume vol;
-        dom_volume_push(vol);
+        for(const auto& l : _leaves)
+          vol.push_back(l.lock()->box().volume());
         return vol;
       }
-      
-      template<class P_,int N_>
-      friend std::ostream& operator<<(std::ostream& os, const PavingBase<P_,N_>& p);
 
-    protected:
-
-      void boxes_list_push(std::list<std::reference_wrapper<const IntervalVector_<N>>>& l, const IntervalVector_<N>& intersect = IntervalVector_<N>()) const
+      std::list<std::weak_ptr<PavingNode<N>>>& leaves()
       {
-        if(is_leaf() && !_x.is_empty() && _x.intersects(intersect))
-          l.push_back(std::cref(_x));
-        else
-        {
-          if(_left) _left->boxes_list_push(l);
-          if(_right) _right->boxes_list_push(l);
-        }
+        return _leaves;
       }
 
-      void leaves_list_push(std::list<P*>& l)
+      std::list<std::reference_wrapper<const IntervalVector_<N>>> boxes_list(const IntervalVector_<N>& intersect = IntervalVector_<N>()) const
       {
-        if(is_leaf() && !_x.is_empty())
-          l.push_back(dynamic_cast<P*>(this));
-        else
-        {
-          if(_left) _left->leaves_list_push(l);
-          if(_right) _right->leaves_list_push(l);
-        }
+        std::list<std::reference_wrapper<const IntervalVector_<N>>> l;
+        _tree->recurse__boxes_list(l, intersect);
+        return l;
       }
 
-      void dom_volume_push(DomainVolume& vol) const
-      {
-        if(is_leaf())
-          vol.push_back(_x.volume());
-        else
-        {
-          if(_left) _left->dom_volume_push(vol);
-          if(_right) _right->dom_volume_push(vol);
-        }
-      }
+      template<int N_>
+      friend class PavingNode;
+      
+      template<int N_>
+      friend std::ostream& operator<<(std::ostream& os, const Paving<N_>& p);
 
-    public: // todo
 
-      IntervalVector_<N> _x;
-      std::shared_ptr<P> _left = nullptr, _right = nullptr;
-  };
+    public: // todo
 
-  template<class P_,int N_>
-  inline std::ostream& operator<<(std::ostream& os, const PavingBase<P_,N_>& p)
-  {
-    os << "Paving";
-    return os;
-  }
+      using leaves_container = std::list<std::weak_ptr<PavingNode<N>>>;
+      leaves_container _leaves;
+      std::shared_ptr<PavingNode<N>> _tree; // must be a shared_ptr to allow enable_shared_from_this
 
-  template<int N>
-  class Paving : public PavingBase<Paving<N>,N>
-  {
     public:
 
-      explicit Paving(size_t n)
-       : PavingBase<Paving<N>,N>(IntervalVector_<N>(n))
-      { }
+      struct const_iterator : public leaves_container::const_iterator
+      {
+        public:
+          
+          const_iterator(typename leaves_container::const_iterator it) : 
+            leaves_container::const_iterator(it) { }
 
-      explicit Paving(const IntervalVector_<N>& x)
-       : PavingBase<Paving<N>,N>(x)
-      { }
+          const IntervalVector_<N>& operator*() const
+          {
+            return (this->leaves_container::const_iterator::operator*()).lock()->box();
+          }
+      };
+
+      auto begin() const { return const_iterator(_leaves.cbegin()); }
+      auto end() const   { return const_iterator(_leaves.cend()); }
   };
 
+  template<int N_>
+  inline std::ostream& operator<<(std::ostream& os, const Paving<N_>& p)
+  {
+    size_t n = p.nb_leaves();
+    os << "Paving (" << n << " box" << (n > 1 ? "es)" : ")") << flush;
+    return os;
+  }
+
 
 } // namespace codac