CarpetLib: Optimise bbox and bboxset classes for speed

Use std::list instead of std::set to store the bboxes internally.

Define some functions as inline.

Perform some checking only when NDEBUG is not defined.

Optimise the algorithms for normalising bboxsets and for calculating
the set difference.

darcs-hash:20070203172717-dae7b-1e77f0a810f786913cd2a1aaed1ea1a5fde604cf.gz

4 files changed, 141 insertions, 100 deletions

diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index e26daf9c7..488b33598 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -13,42 +13,13 @@ using namespace std;
 
 
 
-// Constructors
-template<class T, int D>
-bbox<T,D>::bbox (): _lower(T(1)), _upper(T(0)), _stride(T(1)) { }
-
-template<class T, int D>
-bbox<T,D>::bbox (const bbox& b)
-  : _lower(b._lower), _upper(b._upper), _stride(b._stride)
-{ }
-
-template<class T, int D>
-bbox<T,D>& bbox<T,D>::operator= (const bbox& b) {
-  _lower=b._lower; _upper=b._upper; _stride=b._stride;
-  return *this;
-}
-
-template<class T, int D>
-bbox<T,D>::bbox (const vect<T,D>& lower_, const vect<T,D>& upper_,
-		 const vect<T,D>& stride_)
-  : _lower(lower_), _upper(upper_), _stride(stride_)
-{
-  assert (all(_stride>T(0)));
-  assert (all((_upper-_lower)%_stride == T(0)));
-  if (numeric_limits<T>::is_integer && numeric_limits<T>::is_signed) {
-    // prevent accidental wrap-around
-    assert (all(_lower < numeric_limits<T>::max() / 2));
-    assert (all(_lower > numeric_limits<T>::min() / 2));
-    assert (all(_upper < numeric_limits<T>::max() / 2));
-    assert (all(_upper > numeric_limits<T>::min() / 2));
-  }
-}
-
 // Accessors
 template<class T, int D>
 T bbox<T,D>::size () const {
   if (empty()) return 0;
-//   return prod((shape()+stride()-1)/stride());
+#ifdef NDEBUG
+  return prod(shape()/stride()+1);
+#else
   const vect<T,D> sh((shape()+stride()-T(1))/stride());
   T sz = 1, max = numeric_limits<T>::max();
   for (int d=0; d<D; ++d) {
@@ -57,6 +28,7 @@ T bbox<T,D>::size () const {
     max /= sh[d];
   }
   return sz;
+#endif
 }
 
 // Queries
@@ -114,15 +86,6 @@ bool bbox<T,D>::operator>= (const bbox& b) const {
   return b <= *this;
 }
 
-// Intersection
-template<class T, int D>
-bbox<T,D> bbox<T,D>::operator& (const bbox& b) const {
-  assert (all(stride()==b.stride()));
-  vect<T,D> lo = max(lower(),b.lower());
-  vect<T,D> up = min(upper(),b.upper());
-  return bbox(lo,up,stride());
-}
-
 // Containment
 template<class T, int D>
 bool bbox<T,D>::is_contained_in (const bbox& b) const {
@@ -131,6 +94,13 @@ bool bbox<T,D>::is_contained_in (const bbox& b) const {
   return all(lower()>=b.lower() && upper()<=b.upper());
 }
 
+// Intersection
+template<class T, int D>
+bool bbox<T,D>::intersects (const bbox& b) const {
+  // no alignment check
+  return all(upper()>=b.lower() && lower()<=b.upper());
+}
+
 // Alignment check
 template<class T, int D>
 bool bbox<T,D>::is_aligned_with (const bbox& b) const {
diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh
index 8d0f8f54e..92563651b 100644
--- a/Carpet/CarpetLib/src/bbox.hh
+++ b/Carpet/CarpetLib/src/bbox.hh
@@ -1,7 +1,9 @@
 #ifndef BBOX_HH
 #define BBOX_HH
 
+#include <cassert>
 #include <iostream>
+#include <limits>
 
 #include "defs.hh"
 #include "vect.hh"
@@ -38,17 +40,40 @@ public:
   // Constructors
   
   /** Construct an empty bbox.  */
-  bbox ();
+  bbox ()
+  : _lower(T(1)), _upper(T(0)), _stride(T(1))
+  {
+  }
   
   /** Copy constructor.  */
-  bbox (const bbox& b);
+  bbox (const bbox& b)
+    : _lower(b._lower), _upper(b._upper), _stride(b._stride)
+  {
+  }
   
   /** Assignment operator.  */
-  bbox& operator= (const bbox& b);
+  bbox& operator= (const bbox& b)
+  {
+    _lower=b._lower; _upper=b._upper; _stride=b._stride;
+    return *this;
+  }
   
   /** Create a bbox from bounds and stride.  */
-  bbox (const vect<T,D>& lower, const vect<T,D>& upper,
-	const vect<T,D>& stride);
+  bbox (const vect<T,D>& lower_,
+        const vect<T,D>& upper_,
+	const vect<T,D>& stride_)
+    : _lower(lower_), _upper(upper_), _stride(stride_)
+  {
+    assert (all(_stride>T(0)));
+    assert (all((_upper-_lower)%_stride == T(0)));
+    if (numeric_limits<T>::is_integer and numeric_limits<T>::is_signed) {
+      // prevent accidental wrap-around
+      assert (all(_lower < numeric_limits<T>::max() / 2));
+      assert (all(_lower > numeric_limits<T>::min() / 2));
+      assert (all(_upper < numeric_limits<T>::max() / 2));
+      assert (all(_upper > numeric_limits<T>::min() / 2));
+    }
+}
   
   // Accessors
   // (Don't return references; *this might be a temporary)
@@ -88,11 +113,20 @@ public:
   
   /** Calculate the intersection (the set of common points) with the
       bbox b.  */
-  bbox operator& (const bbox& b) const;
+  bbox operator& (const bbox& b) const
+  {
+    assert (all(stride()==b.stride()));
+    vect<T,D> lo = max(lower(),b.lower());
+    vect<T,D> up = min(upper(),b.upper());
+    return bbox(lo,up,stride());
+  }
   
   /** Find out whether this bbox is contained in the bbox b.  */
   bool is_contained_in (const bbox& b) const;
   
+  /** Find out whether this bbox intersects the bbox b.  */
+  bool intersects (const bbox& b) const;
+  
   /** Find out whether this bbox is aligned with the bbox b.
       ("aligned" means that both bboxes have the same stride and that
       their boundaries are commesurate.)  */
diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index 7725e899b..c36e8a524 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -21,7 +21,8 @@ bboxset<T,D>::bboxset () {
 
 template<class T, int D>
 bboxset<T,D>::bboxset (const box& b) {
-  if (!b.empty()) bs.insert(b);
+  //S if (!b.empty()) bs.insert(b);
+  if (!b.empty()) bs.push_back(b);
   assert (invariant());
 }
 
@@ -41,13 +42,13 @@ bboxset<T,D>::bboxset (const bset& bs_): bs(bs_) {
 template<class T, int D>
 bool bboxset<T,D>::invariant () const {
 // This is very slow when there are many bboxes
-#if 0
+#if 0 && ! defined(CARPET_OPTIMISE)
   for (const_iterator bi=begin(); bi!=end(); ++bi) {
     if ((*bi).empty()) return false;
-    if (! (*bi).is_aligned_with(*bs.begin())) return false;
+    if (not (*bi).is_aligned_with(*bs.begin())) return false;
     // check for overlap (quadratic -- expensive)
     for (const_iterator bi2=begin(); bi2!=bi; ++bi2) {
-      if (! ((*bi2) & (*bi)).empty()) return false;
+      if ((*bi2).intersects(*bi)) return false;
     }
   }
 #endif
@@ -69,15 +70,23 @@ void bboxset<T,D>::normalize ()
   // boundaries aligned.
   for (int d=0; d<D; ++d) {
     // Find all boundaries
-    set<T> sbnds;
+    //S typedef set<T> buf;
+    typedef vector<T> buf;
+    buf sbnds;
+    sbnds.reserve (2 * bs.size());
     for (typename bset::const_iterator si = bs.begin(); si != bs.end(); ++ si) {
       box const & b = * si;
       int const bstr = b.stride()[d];
       int const blo = b.lower()[d];
       int const bhi = b.upper()[d] + bstr;
-      sbnds.insert (blo);
-      sbnds.insert (bhi);
+      //S sbnds.insert (blo);
+      //S sbnds.insert (bhi);
+      sbnds.push_back (blo);
+      sbnds.push_back (bhi);
     }
+    sort (sbnds.begin(), sbnds.end());
+    typename buf::iterator const last = unique (sbnds.begin(), sbnds.end());
+    sbnds.resize (last - sbnds.begin());
     // Split bboxes
     bset nbs;
     for (typename bset::const_iterator si = bs.begin(); si != bs.end(); ++ si) {
@@ -85,24 +94,30 @@ void bboxset<T,D>::normalize ()
       int const bstr = b.stride()[d];
       int const blo = b.lower()[d];
       int const bhi = b.upper()[d] + bstr;
-      typename set<T>::const_iterator const ilo
-        = find (sbnds.begin(), sbnds.end(), blo);
-      typename set<T>::const_iterator const ihi
-        = find (sbnds.begin(), sbnds.end(), bhi);
+//       typename buf::const_iterator const ilo
+//         = find (sbnds.begin(), sbnds.end(), blo);
+//       typename buf::const_iterator const ihi
+//         = find (sbnds.begin(), sbnds.end(), bhi);
+      typename buf::const_iterator const ilo
+        = lower_bound (sbnds.begin(), sbnds.end(), blo);
+      typename buf::const_iterator const ihi
+        = lower_bound (sbnds.begin(), sbnds.end(), bhi);
       assert (ilo != sbnds.end());
       assert (ihi != sbnds.end());
       assert (* ilo == blo);
       assert (* ihi == bhi);
       // Split one bbox
-      for (typename set<T>::const_iterator curr = ilo; curr != ihi; ++ curr) {
-        typename set<T>::const_iterator next = curr;
+      for (typename buf::const_iterator curr = ilo; curr != ihi; ++ curr) {
+        typename buf::const_iterator next = curr;
         advance (next, 1);
         int const nblo = * curr;
         int const nbhi = * next;
+        assert (nbhi > nblo);   // ensure that the set remains sorted
         box nb (b.lower().replace(d, nblo),
                 b.upper().replace(d, nbhi - bstr),
                 b.stride());
-        nbs.insert (nb);
+        //S nbs.insert (nb);
+        nbs.push_back (nb);
       }
     }
     // Replace old set
@@ -142,7 +157,8 @@ void bboxset<T,D>::normalize ()
                           b.stride());
             bs.erase (si);
             nbs.erase (nsi);
-            bs.insert (cb);
+            //S bs.insert (cb);
+            bs.push_back (cb);
             goto done;
           } else if (bhi == nblo) {
             // Combine boxes, b < nb
@@ -151,13 +167,15 @@ void bboxset<T,D>::normalize ()
                           b.stride());
             bs.erase (si);
             nbs.erase (nsi);
-            bs.insert (cb);
+            //S bs.insert (cb);
+            bs.push_back (cb);
             goto done;
           }
         }
       }
       bs.erase (si);
-      nbs.insert (b);
+      //S nbs.insert (b);
+      nbs.push_back (b);
     done:;
     }
     bs.swap (nbs);
@@ -188,22 +206,17 @@ T bboxset<T,D>::size () const {
 
 // Add (bboxes that don't overlap)
 template<class T, int D>
-bboxset<T,D>& bboxset<T,D>::operator+= (const box& b) {
-  if (b.empty()) return *this;
-  // check for overlap
-  for (const_iterator bi=begin(); bi!=end(); ++bi) {
-    assert ((*bi & b).empty());
+bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s) {
+  for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) {
+    *this += *bi;
   }
-  bs.insert(b);
   assert (invariant());
   return *this;
 }
 
 template<class T, int D>
-bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s) {
-  for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) {
-    *this += *bi;
-  }
+bboxset<T,D>& bboxset<T,D>::add_transfer (bboxset& s) {
+  bs.splice (bs.end(), s.bs);
   assert (invariant());
   return *this;
 }
@@ -239,14 +252,16 @@ bboxset<T,D> bboxset<T,D>::plus (const bbox<T,D>& b, const bboxset<T,D>& s) {
 // Union
 template<class T, int D>
 bboxset<T,D>& bboxset<T,D>::operator|= (const box& b) {
-  *this += b - *this;
+  bboxset tmp = b - *this;
+  add_transfer (tmp);
   assert (invariant());
   return *this;
 }
 
 template<class T, int D>
 bboxset<T,D>& bboxset<T,D>::operator|= (const bboxset& s) {
-  *this += s - *this;
+  bboxset tmp = s - *this;
+  add_transfer (tmp);
   assert (invariant());
   return *this;
 }
@@ -290,7 +305,8 @@ bboxset<T,D> bboxset<T,D>::operator& (const bboxset& s) const {
   // walk all the bboxes
   for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) {
     // insert the intersection with this bbox
-    r += *this & *bi;
+    bboxset tmp = *this & *bi;
+    r.add_transfer (tmp);
   }
   assert (r.invariant());
   return r;
@@ -317,8 +333,12 @@ template<class T, int D>
 bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2) {
   assert (b1.is_aligned_with(b2));
   if (b1.empty()) return bboxset<T,D>();
-  if (b2.empty()) return bboxset<T,D>(b1);
-  const vect<T,D> str = b1.stride();
+  if (not b1.intersects(b2)) return bboxset<T,D>(b1);
+  vect<T,D> const b2lo = max (b1.lower(), b2.lower());
+  vect<T,D> const b2up = min (b1.upper(), b2.upper());
+  vect<T,D> const & b1lo = b1.lower();
+  vect<T,D> const & b1up = b1.upper();
+  vect<T,D> const & str = b1.stride();
   bboxset<T,D> r;
   for (int d=0; d<D; ++d) {
     // make resulting bboxes as large as possible in x-direction (for
@@ -327,20 +347,20 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2) {
     bbox<T,D> b;
     for (int dd=0; dd<D; ++dd) {
       if (dd<d) {
-	lb[dd] = b2.lower()[dd];
-	ub[dd] = b2.upper()[dd];
+	lb[dd] = b2lo[dd];
+	ub[dd] = b2up[dd];
       } else if (dd>d) {
-	lb[dd] = b1.lower()[dd];
-	ub[dd] = b1.upper()[dd];
+	lb[dd] = b1lo[dd];
+	ub[dd] = b1up[dd];
       }
     }
-    lb[d] = b1.lower()[d];
-    ub[d] = b2.lower()[d] - str[d];
-    b = bbox<T,D>(lb,ub,str) & b1;
+    lb[d] = b1lo[d];
+    ub[d] = b2lo[d] - str[d];
+    b = bbox<T,D>(lb,ub,str);
     r += b;
-    lb[d] = b2.upper()[d] + str[d];
-    ub[d] = b1.upper()[d];
-    b = bbox<T,D>(lb,ub,str) & b1;
+    lb[d] = b2up[d] + str[d];
+    ub[d] = b1up[d];
+    b = bbox<T,D>(lb,ub,str);
     r += b;
   }
   assert (r.invariant());
@@ -354,18 +374,12 @@ bboxset<T,D> bboxset<T,D>::operator- (const box& b) const {
   // walk all my elements
   for (const_iterator bi=begin(); bi!=end(); ++bi) {
     // insert the difference with the bbox
-    r += *bi - b;
+    bboxset tmp = *bi - b;
+    r.add_transfer (tmp);
   }
   assert (r.invariant());
   return r;
 }
-
-template<class T, int D>
-bboxset<T,D>& bboxset<T,D>::operator-= (const box& b) {
-  *this = *this - b;
-  assert (invariant());
-  return *this;
-}
   
 template<class T, int D>
 bboxset<T,D>& bboxset<T,D>::operator-= (const bboxset& s) {
diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh
index 44eee3c0a..5ad2dea94 100644
--- a/Carpet/CarpetLib/src/bboxset.hh
+++ b/Carpet/CarpetLib/src/bboxset.hh
@@ -3,6 +3,7 @@
 
 #include <cassert>
 #include <iostream>
+#include <list>
 #include <set>
 
 #include "bbox.hh"
@@ -38,7 +39,8 @@ class bboxset {
   
   // Types
   typedef bbox<T,D> box;
-  typedef set<box> bset;
+  //S typedef set<box> bset;
+  typedef list<box> bset;
   
   // Fields
   bset bs;
@@ -67,8 +69,24 @@ public:
   int setsize () const { return bs.size(); }
   
   // Add (bboxes that don't overlap)
-  bboxset& operator+= (const box& b);
+  bboxset& operator+= (const box& b)
+  {
+    if (b.empty()) return *this;
+    // This is very slow when there are many bboxes
+#if 0 && ! defined(CARPET_OPTIMISE)
+    // check for overlap
+    for (const_iterator bi=begin(); bi!=end(); ++bi) {
+      assert (not (*bi).intersects(b));
+    }
+#endif
+    //S bs.insert(b);
+    bs.push_back(b);
+    assert (invariant());
+    return *this;
+  }
+  
   bboxset& operator+= (const bboxset& s);
+  bboxset& add_transfer (bboxset& s);
   bboxset operator+ (const box& b) const;
   bboxset operator+ (const bboxset& s) const;
   static bboxset plus (const box& b1, const box& b2);
@@ -90,7 +108,12 @@ public:
   // friend bboxset operator- <T,D>(const box& b1, const box& b2);
   static bboxset minus (const box& b1, const box& b2);
   bboxset operator- (const box& b) const;
-  bboxset& operator-= (const box& b);
+  bboxset& operator-= (const box& b)
+  {
+    *this = *this - b;
+    assert (invariant());
+    return *this;
+  }
   bboxset& operator-= (const bboxset& s);
   bboxset operator- (const bboxset& s) const;
   // friend bboxset operator- <T,D>(const box& b, const bboxset& s);