[yoda-svn] r428 - in trunk: . include/YODA include/YODA/Utils src

blackhole at projects.hepforge.org blackhole at projects.hepforge.org
Thu Dec 8 18:15:30 GMT 2011


Author: buckley
Date: Thu Dec  8 18:15:30 2011
New Revision: 428

Log:
Adding a Utils::ndarray object and using it to implement a general Scatter<N> system, with generalised Point<N> and Error<N> to boot.

Added:
   trunk/include/YODA/ErrorND.h
   trunk/include/YODA/PointND.h
   trunk/include/YODA/ScatterND.h
   trunk/include/YODA/Utils/ndarray.h
   trunk/src/Scatter1D.cc
Modified:
   trunk/ChangeLog
   trunk/src/Makefile.am

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Thu Dec  8 18:11:55 2011	(r427)
+++ trunk/ChangeLog	Thu Dec  8 18:15:30 2011	(r428)
@@ -1,3 +1,9 @@
+2011-12-08  Andy Buckley  <andy.buckley at cern.ch>
+
+	* Adding a Utils::ndarray object and using it to implement a
+	general Scatter<N> system, with generalised Point<N> and Error<N>
+	to boot.
+
 2011-12-07  Andy Buckley  <andy.buckley at cern.ch>
 
 	* Mapping the Dbn1D and Dbn2D classes into Python.

Added: trunk/include/YODA/ErrorND.h
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/include/YODA/ErrorND.h	Thu Dec  8 18:15:30 2011	(r428)
@@ -0,0 +1,220 @@
+// -*- C++ -*-
+//
+// This file is part of YODA -- Yet more Objects for Data Analysis
+// Copyright (C) 2008-2011 The YODA collaboration (see AUTHORS for details)
+//
+#ifndef YODA_ERRORND_H
+#define YODA_ERRORND_H
+
+#include "YODA/Exceptions.h"
+#include "YODA/Utils/MathUtils.h"
+#include "YODA/Utils/sortedvector.h"
+#include "YODA/Utils/ndarray.h"
+#include <utility>
+#include <algorithm>
+
+namespace YODA {
+
+
+  /// An N-dimensional error to be contained in a Point<N>
+  template<int N>
+  class Error {
+  public:
+
+    typedef Utils::ndarray<double, N> NdVal;
+    typedef Utils::ndarray<std::pair<double,double>, N> NdValPair;
+
+
+    /// @name Constructors
+    //@{
+
+    // Default constructor
+    Error(const std::string& name="")
+      : _name(name)
+    {
+      clear();
+    }
+
+
+    /// Constructor from ND array of error pairs
+    Error(const NdValPair& err, const std::string& name="")
+      : _name(name), _val(err)
+    {    }
+
+
+    /// Constructor of a symmetric error from one ND array of values
+    Error(const NdVal& errsymm, const std::string& name="")
+      : _name(name)
+    {
+      for (size_t i = 0; i < N; ++i) {
+        _val[i] = std::make_pair(errsymm[i], errsymm[i]);
+      }
+    }
+
+
+    /// Constructor of an asymmetric error from two ND arrays of values
+    Error(const NdVal& errminus, const NdVal& errplus, const std::string& name="")
+      : _name(name)
+    {
+      for (size_t i = 0; i < N; ++i) {
+        _val[i] = std::make_pair(errminus[i], errplus[i]);
+      }
+    }
+
+    //@}
+
+
+    /// @name Modifiers and accessors
+    //@{
+
+    /// Clear the point values and errors
+    const std::string& name() const {
+      return _name;
+    }
+
+    /// Clear the point values and errors
+    void setName(const std::string& name) {
+      _name = name;
+    }
+
+    /// Clear the point values and errors
+    void clear() {
+      _val.clear();
+    }
+
+    // /// Get the error pair array
+    // NdValPair& errs() {
+    //   return _val;
+    // }
+
+    // /// Get the error pair array (const version)
+    // const NdValPair& errs() const {
+    //   return _val;
+    // }
+
+    /// Access the error pair in dimension @a dim
+    std::pair<double,double>& err(size_t dim) {
+      return _val[dim];
+    }
+
+    /// Get the error pair in dimension @a dim (const version)
+    const std::pair<double,double>& err(size_t dim) const {
+      return _val[dim];
+    }
+
+    /// Access the error pair in dimension @a dim
+    std::pair<double,double>& operator[](size_t dim) {
+      return _val[dim];
+    }
+
+    /// Get the error pair in dimension @a dim (const version)
+    const std::pair<double,double>& operator[](size_t dim) const {
+      return _val[dim];
+    }
+
+    /// Get the minus error in dimension @a dim
+    double errMinus(size_t dim) const {
+      return _val[dim].first;
+    }
+
+    /// Get the plus error in dimension @a dim
+    double errPlus(size_t dim) const {
+      return _val[dim].second;
+    }
+
+    /// Get the mean error in dimension @a dim
+    double errAvg(size_t dim) const {
+      return (_val[dim].first + _val[dim].second)/2.0;
+    }
+
+    //@}
+
+
+
+    /// @name Scaling and transformations
+    //@{
+
+    /// Uniform scaling
+    void scale(const NdVal& scales) {
+      for (size_t i = 0; i < N; ++i) {
+        _val[i].first *= scales[i];
+        _val[i].second *= scales[i];
+      }
+    }
+
+    /// @todo Generic trf functor support -- need abs position of error bar
+
+    //@}
+
+
+  protected:
+
+    /// @name Value and error variables
+    //@{
+
+    std::string _name;
+    NdValPair _val;
+
+    //@}
+
+  };
+
+
+
+  /// @name Comparison operators
+  //@{
+
+  /// Equality test
+  template <int N>
+  inline bool operator==(const Error<N>& a, const Error<N>& b) {
+    if (a.name() != b.name()) return false;
+    for (size_t i = 0; i < N; ++i) {
+      if ( !fuzzyEquals(a.err()[i], b.err()[i]) ) return false;
+    }
+    return true;
+  }
+
+  /// Inequality test
+  template <int N>
+  inline bool operator!=(const Error<N>& a, const Error<N>& b) {
+    return !(a == b);
+  }
+
+
+  /// Less-than operator used to sort errors
+  template <int N>
+  inline bool operator<(const Error<N>& a, const Error<N>& b) {
+    #define LT_IF_NOT_EQ(a,b) { if (!fuzzyEquals(a, b)) return a < b; }
+    for (size_t i = 0; i < N; ++i) {
+      LT_IF_NOT_EQ(a.err(i).first, b.err(i).first);
+      LT_IF_NOT_EQ(a.err(i).second, b.err(i).second);
+    }
+    #undef LT_IF_NOT_EQ
+    return a.name() < b.name();
+  }
+
+  /// Less-than-or-equals operator used to sort errors
+  template <int N>
+  inline bool operator<=(const Error<N>& a, const Error<N>& b) {
+    if (a == b) return true;
+    return a < b;
+  }
+
+  /// Greater-than operator used to sort errors
+  template <int N>
+  inline bool operator>(const Error<N>& a, const Error<N>& b) {
+    return !(a <= b);
+  }
+
+  /// Greater-than-or-equals operator used to sort errors
+  template <int N>
+  inline bool operator>=(const Error<N>& a, const Error<N>& b) {
+    return !(a < b);
+  }
+
+  //@}
+
+
+}
+
+#endif

Added: trunk/include/YODA/PointND.h
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/include/YODA/PointND.h	Thu Dec  8 18:15:30 2011	(r428)
@@ -0,0 +1,235 @@
+// -*- C++ -*-
+//
+// This file is part of YODA -- Yet more Objects for Data Analysis
+// Copyright (C) 2008-2011 The YODA collaboration (see AUTHORS for details)
+//
+#ifndef YODA_POINTND_H
+#define YODA_POINTND_H
+
+#include "YODA/Exceptions.h"
+#include "YODA/ErrorND.h"
+#include "YODA/Utils/MathUtils.h"
+#include "YODA/Utils/sortedvector.h"
+#include "YODA/Utils/ndarray.h"
+#include <utility>
+#include <algorithm>
+
+namespace YODA {
+
+
+  /// An N-dimensional data point to be contained in a Scatter<N>
+  template<int N>
+  class Point {
+  public:
+
+    // Typedefs
+    typedef Utils::ndarray<double, N> NdVal;
+    typedef Utils::ndarray<std::pair<double,double>, N> NdValPair;
+    typedef Utils::sortedvector< Error<N> > Errors;
+
+
+    /// @name Constructors
+    //@{
+
+    // Default constructor
+    Point() {
+      clear();
+    }
+
+
+    /// Constructor from position values without errors
+    Point(const NdVal& pos)
+      : _pos(pos)
+    {    }
+
+
+    /// Constructor from values with a single set of symmetric errors
+    Point(const NdVal& pos, const NdVal& errs)
+      : _pos(pos)
+    {
+      _errs.push_back(Error<N>(errs));
+    }
+
+
+    /// Constructor from values with a single set of asymmetric errors
+    Point(const NdVal& pos, const NdValPair& errs)
+      : _pos(pos)
+    {
+      _errs.push_back(Error<N>(errs));
+    }
+
+
+    /// Constructor from values and a single Error object
+    Point(const NdVal& pos, const Error<N>& err)
+      : _pos(pos)
+    {
+      _errs.push_back(err);
+    }
+
+
+    /// Constructor from values and a collection of Error objects
+    Point(const std::vector<double>& pos, const std::vector< Error<N> >& errs)
+      : _pos(pos), _errs(errs)
+    {    }
+
+    //@}
+
+
+    /// @name Modifiers
+    //@{
+
+    /// Clear the point values and errors
+    void clear() {
+      for (size_t i = 0; i < N; ++i) _pos[i] = 0;
+      _errs.clear();
+    }
+
+    /// @todo addError, addErrors, setErrors
+
+    //@}
+
+
+  public:
+
+    /// @name Coordinate accessors
+    //@{
+
+    /// Get the coordinate vector
+    NdVal& pos() { return _pos; }
+
+    /// Get the coordinate vector (const version)
+    const NdVal& pos() const { return _pos; }
+
+    /// Set the coordinate vector
+    void setPos(const NdVal& pos) {
+      _pos = pos;
+    }
+
+    //@}
+
+
+    /// @name Error accessors
+    //@{
+
+    /// Get error values
+    Errors& errs() {
+      return _errs;
+    }
+
+    /// Get error values (const version)
+    const Errors& errs() const {
+      return _errs;
+    }
+
+    /// Set the error values
+    void setErrs(const Errors& errs) {
+      _errs = errs;
+    }
+
+    //@}
+
+
+    /// @name Scaling and transformations
+    //@{
+
+    /// Uniform scaling
+    void scale(const NdVal& scales) {
+      for (size_t i = 0; i < N; ++i) {
+        _pos[i] *= scales[i];
+      }
+      foreach (Error<N>& e, errs()) {
+        e.scale(scales);
+      }
+    }
+
+
+    // /// Generalised transformations with functors
+    // void scale(const Trf<N>& trf) {
+    //   for (size_t i = 0; i < N; ++i) {
+    //     _pos = trf.transform(_pos);
+    //   }
+    //   foreach (Error e, errs()) {
+    //     rf.transformErrs(_pos, e);
+    //   }
+    // }
+
+    //@}
+
+
+  protected:
+
+    /// @name Value and error variables
+    //@{
+
+    NdVal _pos;
+    Errors _errs;
+
+    //@}
+
+  };
+
+
+
+  /// @name Comparison operators
+  //@{
+
+  /// Equality test
+  template <int N>
+  inline bool operator==(const Point<N>& a, const Point<N>& b) {
+    // Compare positions
+    for (size_t i = 0; i < N; ++i) {
+      if ( !fuzzyEquals(a.pos()[i], b.pos()[i]) ) return false;
+    }
+    // Compare number of errors and then (sorted) error equality
+    if (a.errs().size() != b.errs().size()) return false;
+    for (size_t i = 0; i < a.errs().size(); ++i) {
+      if (a.errs()[i] != b.errs()[i]) return false;
+    }
+    return true;
+  }
+
+  /// Inequality test
+  template <int N>
+  inline bool operator!=(const Point<N>& a, const Point<N>& b) {
+    return !(a == b);
+  }
+
+
+  /// Less-than operator used to sort points
+  template <int N>
+  inline bool operator<(const Point<N>& a, const Point<N>& b) {
+    #define LT_IF_NOT_EQ(a,b) { if (!fuzzyEquals(a, b)) return a < b; }
+    for (size_t i = 0; i < N; ++i) LT_IF_NOT_EQ(a.pos()[i], b.pos()[i]);
+    LT_IF_NOT_EQ(a.errs().size(), b.errs().size());
+    for (size_t i = 0; i < a.errs().size(); ++i) {
+      if (a.errs()[i] != b.errs()[i]) return a.errs()[i] < b.errs()[i];
+    }
+    #undef LT_IF_NOT_EQ
+    return false;
+  }
+
+  /// Less-than-or-equals operator used to sort points
+  template <int N>
+  inline bool operator<=(const Point<N>& a, const Point<N>& b) {
+    if (a == b) return true;
+    return a < b;
+  }
+
+  /// Greater-than operator used to sort points
+  template <int N>
+  inline bool operator>(const Point<N>& a, const Point<N>& b) {
+    return !(a <= b);
+  }
+
+  /// Greater-than-or-equals operator used to sort points
+  template <int N>
+  inline bool operator>=(const Point<N>& a, const Point<N>& b) {
+    return !(a < b);
+  }
+
+  //@}
+
+
+}
+
+#endif

Added: trunk/include/YODA/ScatterND.h
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/include/YODA/ScatterND.h	Thu Dec  8 18:15:30 2011	(r428)
@@ -0,0 +1,288 @@
+// -*- C++ -*-
+//
+// This file is part of YODA -- Yet more Objects for Data Analysis
+// Copyright (C) 2008-2011 The YODA collaboration (see AUTHORS for details)
+//
+#ifndef YODA_SCATTERND_H
+#define YODA_SCATTERND_H
+
+#include "YODA/AnalysisObject.h"
+#include "YODA/PointND.h"
+#include "YODA/Utils/sortedvector.h"
+#include "YODA/Utils/ndarray.h"
+#include <vector>
+#include <set>
+#include <string>
+#include <utility>
+
+namespace YODA {
+
+
+  /// A very generic data type which is just a collection of ND data points with errors
+  template <int N>
+  class Scatter : public AnalysisObject {
+  public:
+
+    // Typedefs
+    typedef Utils::ndarray<double, N> NdVal;
+    typedef Utils::ndarray<std::pair<double,double>, N> NdValPair;
+    typedef Utils::sortedvector< Point<N> > Points;
+
+
+    /// @name Constructors
+    //@{
+
+    /// Empty constructor
+    Scatter(const std::string& path="", const std::string& title="")
+      : AnalysisObject("Scatter", path, title)
+    {  }
+
+
+    /// Constructor from a set of points
+    Scatter(const Points& points,
+            const std::string& path="", const std::string& title="")
+      : AnalysisObject("Scatter", path, title),
+        _points(points)
+    {  }
+
+
+    /// Constructor from a vector of position values with no errors
+    Scatter(const std::vector<NdVal>& positions,
+            const std::string& path="", const std::string& title="")
+      : AnalysisObject("Scatter", path, title)
+    {
+      for (size_t i = 0; i < positions.size(); ++i) {
+        addPoint(Point<N>(positions[i]));
+      }
+    }
+
+
+    /// Constructor from vectors of values for positions and a single set of symmetric errors
+    Scatter(const std::vector<NdVal>& positions,
+            const std::vector<NdVal>& errors,
+            const std::string& path="", const std::string& title="")
+      : AnalysisObject("Scatter", path, title)
+    {
+      assert(positions.size() == errors.size());
+      for (size_t i = 0; i < positions.size(); ++i) {
+        addPoint(Point<N>(positions[i], errors[i]));
+      }
+    }
+
+
+    // /// Constructor from values with completely explicit asymmetric errors
+    // Scatter(const std::vector<double>& x, const std::vector<double>& y,
+    //         const std::vector<double>& exminus,
+    //         const std::vector<double>& explus,
+    //         const std::vector<double>& eyminus,
+    //         const std::vector<double>& eyplus,
+    //         const std::string& path="", const std::string& title="")
+    //   : AnalysisObject("Scatter", path, title)
+    // {
+    //   assert(x.size() == y.size() &&
+    //          x.size() == exminus.size() && x.size() == explus.size() &&
+    //          x.size() == eyminus.size() && x.size() == eyplus.size());
+    //   for (size_t i = 0; i < x.size(); ++i) {
+    //     addPoint(Point<N>2D(x[i], exminus[i], explus[i], y[i], eyminus[i], eyplus[i]));
+    //   }
+    // }
+
+
+    /// Copy constructor with optional new path
+    Scatter(const Scatter<N>& s, const std::string& path="")
+      : AnalysisObject("Scatter", (path.size() == 0) ? s.path() : path, s, s.title()),
+        _points(s._points)
+    {  }
+
+
+    /// Assignment operator
+    Scatter<N>& operator = (const Scatter<N>& s) {
+      setPath(s.path());
+      setTitle(s.title());
+      _points = s._points;
+      return *this;
+    }
+
+    //@}
+
+
+    /// @name Modifiers
+    //@{
+
+    /// Clear all points
+    void reset() {
+      _points.clear();
+    }
+
+    /// Scaling
+    void scale(const NdVal& scales) {
+      foreach (Point<N>& p, _points) {
+        p.scale(scales);
+      }
+    }
+
+    //@}
+
+
+    ///////////////////////////////////////////////////
+
+
+    /// @name Point accessors
+    //@{
+
+    /// Number of points in the scatter
+    size_t numPoints() const {
+      return _points.size();
+    }
+
+
+    /// Get the collection of points
+    Points& points() {
+      return _points;
+    }
+
+
+    /// Get the collection of points (const version)
+    const Points& points() const {
+      return _points;
+    }
+
+
+    /// Get a reference to the point with index @a index
+    Point<N>& point(size_t index) {
+      assert(index < numPoints());
+      return _points.at(index);
+    }
+
+
+    /// Get the point with index @a index (const version)
+    const Point<N>& point(size_t index) const {
+      assert(index < numPoints());
+      return _points.at(index);
+    }
+
+    //@}
+
+
+    /// @name Point inserters
+    //@{
+
+    /// Insert a new point
+    Scatter<N>& addPoint(const Point<N>& pt) {
+      _points.insert(pt);
+      return *this;
+    }
+
+    /// Insert a new point, from a position array
+    Scatter<N>& addPoint(const NdVal& pos) {
+      _points.insert(Point<N>(pos));
+      return *this;
+    }
+
+    /// @todo More addPoint combinations with arrays for errors
+
+    /// Insert a collection of new points
+    Scatter<N>& addPoints(Points pts) {
+      foreach (const Point<N>& pt, pts) {
+        addPoint(pt);
+      }
+      return *this;
+    }
+
+    //@}
+
+
+    /// @name Combining sets of scatter points
+    //@{
+
+    /// @todo Better name?
+    Scatter<N>& combineWith(const Scatter<N>& other) {
+      addPoints(other.points());
+      return *this;
+    }
+
+    /// @todo Better name?
+    Scatter<N>& combineWith(const std::vector< Scatter<N> >& others) {
+      foreach (const Scatter<N>& s, others) {
+        combineWith(s);
+      }
+      return *this;
+    }
+
+    //@}
+
+
+  private:
+
+    Points _points;
+
+  };
+
+
+  /// @name Combining scatters by merging sets of points
+  //@{
+
+  template <int N>
+  inline Scatter<N> combine(const Scatter<N>& a, const Scatter<N>& b) {
+    Scatter<N> rtn = a;
+    rtn.combineWith(b);
+    return rtn;
+  }
+
+  template <int N>
+  inline Scatter<N> combine(const std::vector< Scatter<N> >& scatters) {
+    Scatter<N> rtn;
+    rtn.combineWith(scatters);
+    return rtn;
+  }
+
+  //@}
+
+
+  //////////////////////////////////
+
+
+  // /// @name Combining scatters: global operators, assuming aligned points
+  // //@{
+
+  // /// Add two scatters
+  // template <int N>
+  // Scatter add(const Scatter& first, const Scatter& second);
+
+
+  // /// Add two scatters
+  // template <int N>
+  // inline Scatter operator + (const Scatter& first, const Scatter& second) {
+  //   return add(first, second);
+  // }
+
+
+  // /// Subtract two scatters
+  // template <int N>
+  // Scatter subtract(const Scatter& first, const Scatter& second);
+
+
+  // /// Subtract two scatters
+  // template <int N>
+  // inline Scatter operator - (const Scatter& first, const Scatter& second) {
+  //   return subtract(first, second);
+  // }
+
+
+  // /// Divide two scatters
+  // template <int N>
+  // Scatter divide(const Scatter& numer, const Scatter& denom);
+
+
+  // /// Divide two scatters
+  // template <int N>
+  // inline Scatter operator / (const Scatter& numer, const Scatter& denom) {
+  //   return divide(numer, denom);
+  // }
+
+  // //@}
+
+
+}
+
+#endif

Added: trunk/include/YODA/Utils/ndarray.h
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/include/YODA/Utils/ndarray.h	Thu Dec  8 18:15:30 2011	(r428)
@@ -0,0 +1,88 @@
+#ifndef YODA_ARRAY_H
+#define YODA_ARRAY_H
+
+#include <vector>
+#include <algorithm>
+#include <sstream>
+#include "boost/array.hpp"
+
+namespace YODA {
+  namespace Utils {
+
+
+    /// @brief Fixed-size array with constructors from standard array types
+    template <typename T, int N>
+    class ndarray {
+    public:
+
+      /// @name Constructors
+      //@{
+
+      /// Default constructor
+      ndarray() {
+        clear();
+      }
+
+      /// Conversion from std::vector
+      ndarray(const std::vector<T>& vec) {
+        if (vec.size() != N) {
+          std::stringstream msg;
+          msg << "Value vector of wrong size supplied to a " << N << " dimensional array";
+          throw RangeError(msg.str());
+        }
+        std::copy(vec, vec+N, _val);
+      }
+
+
+      /// Conversion from Boost array
+      ndarray(const boost::array<T,N> arr) {
+        _val = arr;
+      }
+
+
+      /// Conversion from C array
+      ndarray(const T arr[5]) {
+        try {
+          std::copy(arr, arr+N, _val);
+        } catch (...) {
+          std::stringstream msg;
+          msg << "Value vector of wrong size supplied to a " << N << " dimensional array";
+          throw RangeError(msg.str());
+        }
+      }
+
+      //@}
+
+
+      /// @name Accessors and modifiers
+      //@{
+
+      /// Clear the point values and errors
+      void clear() {
+        for (size_t i = 0; i < N; ++i) _val[i] = T();
+      }
+
+      /// Array-like accessor
+      T& operator[](size_t i) {
+        return _val[i];
+      }
+
+      /// Array-like accessor (const version)
+      const T& operator[](size_t i) const {
+        return _val[i];
+      }
+
+      //@}
+
+
+    private:
+
+      boost::array<T,N> _val;
+
+    };
+
+
+  }
+}
+
+#endif

Modified: trunk/src/Makefile.am
==============================================================================
--- trunk/src/Makefile.am	Thu Dec  8 18:11:55 2011	(r427)
+++ trunk/src/Makefile.am	Thu Dec  8 18:15:30 2011	(r428)
@@ -8,6 +8,7 @@
     Histo2D.cc \
     Profile1D.cc \
     Profile2D.cc \
+    Scatter1D.cc \
     Scatter2D.cc \
     Scatter3D.cc\
     Writer.cc \

Added: trunk/src/Scatter1D.cc
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ trunk/src/Scatter1D.cc	Thu Dec  8 18:15:30 2011	(r428)
@@ -0,0 +1,11 @@
+#include "YODA/ScatterND.h"
+
+namespace YODA {
+
+
+  void test() {
+    Scatter<2> s2;
+  }
+
+
+}


More information about the yoda-svn mailing list