[yoda-svn] r444 - in trunk: . include/YODA pyext/yoda/include

[blackhole at projects.hepforge.org]

Author: buckley
Date: Wed May  2 12:06:24 2012
New Revision: 444

Log:
Improvements (I hope) to the binary search in Axis1D, and providing an experimental default constructor for Histo1D.

Modified:
   trunk/ChangeLog
   trunk/include/YODA/Axis1D.h
   trunk/include/YODA/Histo1D.h
   trunk/pyext/yoda/include/40-Histo1D.pyx
   trunk/pyext/yoda/include/40-Profile1D.pyx

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Wed May  2 11:30:48 2012	(r443)
+++ trunk/ChangeLog	Wed May  2 12:06:24 2012	(r444)
@@ -1,3 +1,8 @@
+2012-05-02  Andy Buckley  <andy.buckley at cern.ch>
+
+	* Improvements (I hope) to the binary search in Axis1D, and
+	providing an experimental default constructor for Histo1D.
+
 2011-12-08  Hendrik Hoeth  <hendrik.hoeth at cern.ch>
 
 	* ReaderYODA can now parse Histo1D and Profile1D flat files

Modified: trunk/include/YODA/Axis1D.h
==============================================================================
--- trunk/include/YODA/Axis1D.h	Wed May  2 11:30:48 2012	(r443)
+++ trunk/include/YODA/Axis1D.h	Wed May  2 12:06:24 2012	(r444)
@@ -26,7 +26,7 @@
     /// only for bins storage.
     typedef typename std::vector<BIN1D> Bins;
 
-    /// A standard edge container. First number defines location
+    /// A bin edge container. First number defines location
     /// of the edge, the second one a bin number the edge is a member of.
     typedef typename std::pair<double, size_t> Edge;
 
@@ -167,18 +167,11 @@
 
     /// Returns an index of a bin at a given coord, -1 if no bin.
     int getBinIndex(double coord) const {
-      // This is NOT a magic number, it is merely a 'trick' to yield a
-      // correct answer in the case of coord pointing exactly on an edge.
-      // This will never change the answer, since it is *much* smaller than
-      // the fuzzyEquals() tolerance!
-      /// @todo No it isn't: the fuzzyEquals tolerance is relative. You don't know that the typical scale of coord isn't 10e-10. This has to go.
-      coord += 0.00000000001;
-
       // Search for an edge
       size_t index = _binaryS(coord, 0, _binHashSparse.size());
 
       // If lower bound is a last edge, it means that we are 'above' the
-      // axis and no bin canbe found in such case. Therefore, announce it.
+      // axis and no bin can be found in such case. Therefore, announce it.
       if (index == _binHashSparse.size() - 1) return -1;
 
       // If on the left to the point in consideration there is an edge that is
@@ -389,20 +382,28 @@
     }
 
 
-    // Binary search algorithm. For an in-depth explanation
-    // please see the documentation of an analogous one in Axis2D
+    /// Binary search algorithm
+    ///
+    /// We are searching for a containing bin for 'value', located in
+    /// the range of bins from indices 'lower' to 'higher'. Failure to
+    /// find 'value' in the provided range will result in an assertion
+    /// failure. A value exactly on a bin edge will be accepted into the
+    /// bin on the RHS of the edge, unless it is the rightmost edge.
     size_t _binaryS(double value, size_t lower, size_t higher) const {
+      assert(_binHashSparse[lower].first <= value && _binHashSparse[higher].first >= value);
       if (lower == higher) return lower;
-      size_t where = (higher+lower)/2;
+      size_t where = (higher+lower)/2; // binary split point
+
+      /// @todo Explicitly throw if coord is exactly equal to an edge value?
 
       if (value >= _binHashSparse[where].first) {
         if (where == _binHashSparse.size() - 1) return where;
-        if (value <= _binHashSparse[where+1].first) return where;
+        if (value < _binHashSparse[where+1].first) return where;
         return _binaryS(value, where, higher);
+      } else {
+        if (where == 0) return where;
+        return _binaryS(value, lower, where);
       }
-
-      if (where == 0) return where;
-      return _binaryS(value, lower, where);
     }
 
 

Modified: trunk/include/YODA/Histo1D.h
==============================================================================
--- trunk/include/YODA/Histo1D.h	Wed May  2 11:30:48 2012	(r443)
+++ trunk/include/YODA/Histo1D.h	Wed May  2 12:06:24 2012	(r444)
@@ -37,6 +37,14 @@
     /// @name Constructors
     //@{
 
+
+    /// Default constructor of an invalid histo
+    Histo1D()
+      : AnalysisObject("Histo1D", "", ""),
+        _axis()
+    { }
+
+
     /// Constructor giving range and number of bins.
     Histo1D(size_t nbins, double lower, double upper,
             const std::string& path="", const std::string& title="")

Modified: trunk/pyext/yoda/include/40-Histo1D.pyx
==============================================================================
--- trunk/pyext/yoda/include/40-Histo1D.pyx	Wed May  2 11:30:48 2012	(r443)
+++ trunk/pyext/yoda/include/40-Histo1D.pyx	Wed May  2 12:06:24 2012	(r444)
@@ -4,6 +4,7 @@
     #cScatter2D operator / (cHisto1D &, cHisto1D &)"""
 
     cdef cppclass cHisto1D "YODA::Histo1D"(cAnalysisObject):
+        cHisto1D()
         cHisto1D(size_t nbins, double lower, double upper, string &path, string &title)
         cHisto1D(vector[double] &binedges, string &path, string &title)
         cHisto1D(vector[double] &binedges)

Modified: trunk/pyext/yoda/include/40-Profile1D.pyx
==============================================================================
--- trunk/pyext/yoda/include/40-Profile1D.pyx	Wed May  2 11:30:48 2012	(r443)
+++ trunk/pyext/yoda/include/40-Profile1D.pyx	Wed May  2 12:06:24 2012	(r444)
@@ -55,7 +55,6 @@
 
         if len(args) == 3:
             nbins, lower, upper = args[0], args[1], args[2]
-
             self.setptr(
                 new cProfile1D(nbins, lower, upper, string(path), string(title))
             )