|
|
[yoda-svn] r248 - in trunk: . include/YODA src tests
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue Aug 16 10:35:16 BST 2011
Author: mkawalec
Date: Tue Aug 16 10:35:15 2011
New Revision: 248
Log:
(Re)implemented division, fixed equality operators, added come comments, fixed the order of stuff in Axis2D, added merge operator (passes prelimenary tests, but don't treat it as a working product yet!), added a variable determining if a bin is Real or Imaginary. Fixed SVN conflicts.
Modified:
trunk/TODO
trunk/include/YODA/Axis2D.h
trunk/include/YODA/Bin2D.h
trunk/include/YODA/Histo2D.h
trunk/include/YODA/Scatter3D.h
trunk/src/Bin2D.cc
trunk/src/Histo1D.cc
trunk/src/Histo2D.cc
trunk/src/Scatter3D.cc
trunk/tests/TestHisto2D.cc
Modified: trunk/TODO
==============================================================================
--- trunk/TODO Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/TODO Tue Aug 16 10:35:15 2011 (r248)
@@ -16,10 +16,6 @@
exception if binnings are not complete grids. (MICHAL)
(done -- to be reviewed)
-* Remove getBinHash from Axis2D and use private bin hash member for equality
- comparisons (MICHAL)
- AB: Is this done?
-
* Rebinning: merges of n adjacent bins and global rebinning by integer factor
(on widths or on bin groups?) (AB for 1D, MICHAL for 2D)
Modified: trunk/include/YODA/Axis2D.h
==============================================================================
--- trunk/include/YODA/Axis2D.h Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/include/YODA/Axis2D.h Tue Aug 16 10:35:15 2011 (r248)
@@ -16,19 +16,7 @@
using namespace std;
-/// @todo Michal, please work through this file for the @todo markers. They are
-/// hard to see in the raw diff because my editor strips trailing
-/// whitespace. There are lots of code quality and consistency issues to be
-/// addressed, and we need to do it NOW before anyone else is going to use YODA
-/// again...
-
-
-/// A big number for the low/high search:
-/// @todo This isn't good enough! Why would this be guaranteed to be
-/// larger than the axis scale? Use e.g. std::limits<double> and/or determine
-/// the true extent when booking.
-/// @todo Also, the code convention is that _foo is a private member variable/function:
-/// for constants use all-caps, e.g. LARGENUM.
+/// Big and small number for low/high search:
const double LARGENUM = numeric_limits<double>::max();
const double SMALLNUM = numeric_limits<double>::min();
@@ -37,8 +25,6 @@
/// @todo Use the same 2-space indentation scheme as elsewhere
- /// @todo Only use C++ //-type comments, please
-
/// @brief 2D bin container and provider
/// This class handles almost all boiler-plate operations
/// on 2D bins (like creating axis, adding, searching, testing).
@@ -66,581 +52,288 @@
/// Segment, having a beginning and end.
typedef typename std::pair<Point, Point> Segment;
- private:
-
- /// @todo The class layout that we're using everywhere else is to have the
- /// public sections at the top, with the constructors first: please use
- /// the same structure so that we can find what we're looking for.
+ public:
+ /// @name Constructors:
+ //@{
- /// @brief Segment validator function
- /// This a 'dispatcher' function. It checks if the segment in question
- /// is vertical or horizontal and launches the appropriate function
- /// searching for cuts in the prope direction. Since it operates on
- /// a vector of segments it is prepared to act on arbitrarly large sets of edges,
- /// in practice usually being four sides of a rectangular bin.
- ///
- /// Notice that it will never be checked, in the current state, if there is a cut
- /// in edges in the edgeset. This imposes the requirement of provide the program
- /// with non-degenerate bins, until checking in implemented. However, it doesn't seem
- /// to be needed, as edges are not generated by a user.
- ///
- /// This is also a perfect place to paralellize the program, if required.
- bool _validateEdge(vector<Segment>& edges)
+ /// @brief Empty constructor
+ /// Only added because it is required by SWIG.
+ /// It doesn't make much sense to use it.
+ Axis2D()
{
- /// Setting the return variable. True means that no cuts were detected.
- bool ret = true;
-
- /// Looping over all the edges provided
- for (unsigned int i=0; i < edges.size(); ++i) {
- /// If the X coordinate of the starting point is the same
- /// as X coordinate of the ending one, checks if there are cuts
- /// on this vertical segment.
- if (fuzzyEquals(edges[i].first.first, edges[i].second.first)) ret = _findCutsY(edges[i]);
-
- /// Check if the segment is horizontal and is it cutting any bin that already exists
- else if (fuzzyEquals(edges[i].first.second, edges[i].second.second)) ret = _findCutsX(edges[i]);
-
- /// This is a check that discards the bin if it is not a rectangle
- /// composed of vertical and horizontal segments.
- else ret = false;
-
- /// If a cut was detected, say it. There is no point in checking other edges
- /// in the set.
- if (!ret) return false;
- }
- /// If no cuts were detected in any of the edges, tell the launching function about this
- return true;
+ vector<Segment> edges;
+ _mkAxis(edges);
}
- /// @brief A binary search function
- /// This is conceptually the same implementation as in STL
- /// but because it returns the index of an element in a vector,
- /// it is easier to use in our case than the STL implementation
- /// that returns a pointer at the element.
- size_t _binaryS(Utils::cachedvector<EdgeCollection>& toSearch,
- double value, size_t lower, size_t higher)
+ /// Constructor provided with a vector of bin delimiters
+ Axis2D(const vector<Segment>& binLimits)
{
- /// Just a check if such degenerate situation happens
- if(lower == higher) return lower;
+ _mkAxis(binLimits);
+ }
- /// Choose a midpoint that will be our pivot
- size_t where = (higher+lower)/2;
+ ///Most standard constructor, should be self-explanatory
+ Axis2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY)
+ {
+ vector<Segment> binLimits;
+ double coeffX = (upperX - lowerX)/(double)nbinsX;
+ double coeffY = (upperY - lowerX)/(double)nbinsY;
- /// Launch the same procedure on half of the range above the pivot
- if(value >= toSearch[where].first) {
- if(where == toSearch.size() - 1) return where;
- if(value <= toSearch[where+1].first) return where;
- return _binaryS(toSearch, value, where, higher);
+ for(double i=lowerX; i < upperX; i+=coeffX) {
+ for(double j=lowerY; j < upperY; j+=coeffY) {
+ binLimits.push_back(make_pair(make_pair(i, j),
+ make_pair((double)(i+coeffX), (double)(j+coeffY))));
+ }
}
+ _mkAxis(binLimits);
+ }
+ //@}
- /// This is not a redundant check, because
- /// of the nature of int division.
- if (where == 0) return where;
+ /// @name Addition operators:
+ //@{
- /// Check if the value is somewhere inbetween
- /// an element at the position in question and
- /// an element at a lower position. If so, return
- /// an index to the current positon.
- if(value >= toSearch[where-1].first) return where;
+ /// @brief Bin addition operator
+ /// This operator is provided a vector of limiting
+ /// points in the format required by _mkAxis().
+ /// It should be noted that there is nothing special about
+ /// the initiation stage of Axis2D, and the edges can be added
+ /// online if they meet all the requirements of non-degeneracy.
+ /// No merging is supported, and I don't think it should before the support
+ /// for merging for '+' operator (codewise it should be the same thing).
+ void addBin(const vector<Segment>& binLimits)
+ {
+ _mkAxis(binLimits);
+ }
- /// If none of the above occurs, the value must
- /// be smaller that the element at the current
- /// position. In such case, launch the search on
- /// half of the interval below the current position.
- return _binaryS(toSearch, value, lower, where);
+ /// @brief Bin addition operator
+ /// This operator is supplied with whe extreamal coordinates of just
+ /// one bin. It then launches the standard bin addition procedure.
+ void addBin(double lowX, double lowY, double highX, double highY)
+ {
+ vector<Segment> coords;
+ coords.push_back(make_pair(make_pair(lowX, lowY), make_pair(highX, highY)));
+
+ addBin(coords);
}
+ //@}
- /// @brief Function that finds cuts of horizontal edges.
- /// A specialised function that tries to exploit the fact
- /// that edges can cut one another only on right angles
- /// to the highest extent. The inner workings are explained
- /// in the comments placed in the function body.
- bool _findCutsX(Segment& edge) {
- /// Look up the limits of search in the _binHashSparse
- /// structure. We are not interested in the vertical edges
- /// that are before the beginning of our segment or after its end.
- size_t i = _binaryS(_binHashSparse.second, edge.first.first, 0, _binHashSparse.second.size());
- size_t end = _binaryS(_binHashSparse.second, edge.second.first, 0, _binHashSparse.second.size());
+ /// @name Some helper functions:
+ //@{
- for(; i < end; i++) {
+ /// @brief Bin merger
+ /// Try to merge a certain amount of bins
+ void mergeBins(size_t from, size_t to) {
+ HistoBin2D& start = bin(from);
+ HistoBin2D& end = bin(to);
+ HistoBin2D temp = start;
+ start.isReal = false;
+
+ if(start.midpoint().first > end.midpoint().first ||
+ start.midpoint().second > end.midpoint().second)
+ throw GridError("The start/end bins are wrongly defined.");
+
+
+ for(size_t y = (*_binHashSparse.first._cache.lower_bound(start.yMin())).second; y <= (*_binHashSparse.first._cache.lower_bound(end.yMin())).second; ++y)
+ {
+ for(size_t x = 0; x < _binHashSparse.first[y].second.size(); ++x) {
+ if((_binHashSparse.first[y].second[x].second.first > start.xMin() ||
+ fuzzyEquals(_binHashSparse.first[y].second[x].second.first, start.xMin())) &&
+ (_binHashSparse.first[y].second[x].second.second < end.xMax() ||
+ fuzzyEquals(_binHashSparse.first[y].second[x].second.second, end.xMax())) &&
+ bin(_binHashSparse.first[y].second[x].first).isReal)
+ {
+ temp += bin(_binHashSparse.first[y].second[x].first);
+ bin(_binHashSparse.first[y].second[x].first).isReal = false;
+
+ /// @todo The bin that was just added must be changed to virtual!
+ }
+ }
+ }
+ _addEdge(temp.edges(), _binHashSparse, false);
+ _bins.push_back(temp);
+
+ _binHashSparse.first.regenCache();
+ _binHashSparse.second.regenCache();
+ _regenDelimiters();
+ }
+ /// @brief Checks if our bins form a grid.
+ /// This function uses a neat property of _binHashSparse.
+ /// If it is containing a set of edges forming a grid without
+ /// gaps in the middle it will have the same number of edges in the
+ /// inner subcaches and half of this amount in the outer (grid boundary)
+ /// subcaches. This makes isGriddy() a very, very fast function.
+ /// @todo Is the name appropriate?
+ int isGriddy() const
+ {
- /// Scroll through all the vertical segments with a given X coordinate
- /// and see if any of those fulfills the cutting requirement. If it does,
- /// announce it.
- for(unsigned int j = 0; j < _binHashSparse.second[i].second.size(); j++) {
- /// Note that we are not taking into account the edges touching the
- /// segment in question. That's because sides of a bin touch
- if(_binHashSparse.second[i].second[j].second.first < edge.first.second &&
- _binHashSparse.second[i].second[j].second.second > edge.first.second &&
- !fuzzyEquals(_binHashSparse.second[i].second[j].second.first, edge.first.second) &&
- !fuzzyEquals(_binHashSparse.second[i].second[j].second.second, edge.first.second)) {
- return false;
- }
+ /// Check if the number of edges parallel to X axis
+ /// is proper in every subcache.
+ unsigned int sizeX = _binHashSparse.first[0].second.size();
+ for(unsigned int i=1; i < _binHashSparse.first.size(); i++) {
+ if(i == _binHashSparse.first.size() - 1) {
+ if(_binHashSparse.first[i].second.size() != sizeX) {
+ return -1;
}
}
- /// If none of the existing edges is cutting this edge, announce it
- return true;
- }
-
- /// @brief Function that finds cuts of vertical edges
- /// For a detailed descpription, please look into
- /// documentation for _findCutsX().
- bool _findCutsY(Segment& edge) {
- size_t i = _binaryS(_binHashSparse.first, edge.first.second, 0, _binHashSparse.first.size());
- size_t end = _binaryS(_binHashSparse.first, edge.second.second, 0, _binHashSparse.first.size());
- for(; i < end; i++) {
+ else if(_binHashSparse.first[i].second.size() != 2*sizeX) {
+ return -1;
+ }
+ }
- for(unsigned int j = 0; j < _binHashSparse.first[i].second.size(); j++) {
- if(_binHashSparse.first[i].second[j].second.first < edge.first.first &&
- _binHashSparse.first[i].second[j].second.second > edge.first.first &&
- !fuzzyEquals(_binHashSparse.first[i].second[j].second.first, edge.first.first) &&
- !fuzzyEquals(_binHashSparse.first[i].second[j].second.second, edge.first.first)) {
- return false;
- }
+ /// Do the same for edges parallel to Y axis.
+ unsigned int sizeY = _binHashSparse.second[0].second.size();
+ for(unsigned int i=1; i < _binHashSparse.second.size(); i++) {
+ if(i!= _binHashSparse.second.size() - 1) {
+ if(2*sizeY != _binHashSparse.second[i].second.size()) {
+ return -1;
}
}
- return true;
- }
+ else if(_binHashSparse.second[i].second.size() != sizeY) return -1;
+ }
- /// @brief Function executed when a set of edges is dropped.
- /// It does not have any information about which edge in the set
- /// had failed the check. If this is needed such additional information
- /// can be readily implemented.
- void _dropEdge(vector<Segment>& edges) {
- std::cerr << "A set of edges was dropped." << endl;
+ /// If everything is proper, announce it.
+ return 0;
}
- /// @brief Bin adder.
- /// It contains all the commands that need to executed
- /// to properly add a bin. Specifially edges are added to
- /// the edge cache (_binHashSparse) and a bin is created from
- /// those edges.
- void _addEdge(vector<Segment>& edges, pair<Utils::cachedvector<EdgeCollection>,
- Utils::cachedvector<EdgeCollection> >& binHash, bool addBin = true) {
- /// Check if there was no mistake made when adding segments to a vector.
- if(edges.size() != 4) throw Exception("The segments supplied don't describe a full bin!");
+ /// Return a total number of bins in a Histo
+ unsigned int numBinsTotal() const
+ {
+ return _bins.size();
+ }
- /// This is the part in charge of adding each of the segments
- /// to the edge cache. Segments are assumed to be validated
- /// beforehand.
- for(unsigned int j=0; j < edges.size(); j++) {
- /// Association made for convinience.
- Segment edge = edges[j];
+ /// Get inf(X) (non-const version)
+ double lowEdgeX()
+ {
+ return _lowEdgeX;
+ }
- /// Do the following if the edge is vertical.
- /// Those two cases need to be distinguished
- /// because of the way in which edge cache is structured.
- if(edge.first.first == edge.second.first) {
- /// See if our edge has the same X coordinate as any other
- /// edge that is currently in the cache.
+ /// Get sup(X) (non-const version)
+ double highEdgeX()
+ {
+ return _highEdgeX;
+ }
- /// Keeps the status of the search
- bool found = false;
+ /// Get inf(Y) (non-const version)
+ double lowEdgeY()
+ {
+ return _lowEdgeY;
+ }
- /// There is only a certain set of X coordinates that we need to sweep
- /// to check if our segment has the same X coordinate. Find them.
- size_t i = _binaryS(binHash.second, edge.first.first, 0, binHash.second.size())-1;
- if(i < 0) i = 0;
- size_t end = i+3;
+ /// Get sup(Y) (non-const version)
+ double highEdgeY()
+ {
+ return _highEdgeY;
+ }
- /// For the coordinates in range, check if one of them is an X coordinate of
- /// the sement.
- for(; i < binHash.second.size() && i < end ; i++) {
- /// If this is the case, do what is needed to be done.
- if(fuzzyEquals(binHash.second[i].first, edge.first.first)) {
- binHash.second[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.second, edge.second.second)));
- found = true;
- break;
- }
- }
-
- /// If no edge with the same X coordinate exist, create
- /// a new subhash at the X coordinate of a segment.
- if(!found) {
- vector<Edge> temp;
- temp.push_back(make_pair(_bins.size(), make_pair(edge.first.second, edge.second.second)));
- binHash.second.push_back(make_pair(edge.first.first,temp));
- sort(binHash.second.begin(), binHash.second.end());
- }
- }
-
- /// See the vertical case for description of a horizontal one
- else if(edge.first.second == edge.second.second) {
- bool found = false;
- size_t i = _binaryS(binHash.first, edge.first.second, 0, binHash.first.size())-1;
- if(i < 0) i = 0;
- size_t end = i+3;
- for(; i < binHash.first.size() && i < end; i++) {
- if(fuzzyEquals(binHash.first[i].first, edge.first.second)) {
- binHash.first[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.first, edge.second.first)));
- found = true;
- }
- }
- if(!found) {
- vector<Edge> temp;
- temp.push_back(make_pair(_bins.size(), make_pair(edge.first.first, edge.second.first)));
- binHash.first.push_back(make_pair(edge.second.second, temp));
- sort(binHash.first.begin(), binHash.first.end());
- }
- }
- }
- /// Now, create a bin with the edges provided
- if(addBin) _bins.push_back(BIN(edges));
+ /// Get inf(X) (const version)
+ const double lowEdgeX() const
+ {
+ return _lowEdgeX;
}
- /// @brief Orientation fixer
- /// Check if the orientation of an edge is proper
- /// for the rest of the algorithm to work on, and if it is not
- /// fix it.
- void _fixOrientation(Segment& edge) {
- if(fuzzyEquals(edge.first.first, edge.second.first)) {
- if(edge.first.second > edge.second.second) {
- double temp = edge.second.second;
- edge.second.second = edge.first.second;
- edge.first.second = temp;
- }
- }
- else if(edge.first.first > edge.second.first) {
- double temp = edge.first.first;
- edge.first.first = edge.second.first;
- edge.second.first = temp;
- }
+ /// Get sup(X) (const version)
+ const double highEdgeX() const
+ {
+ return _highEdgeX;
}
- /// @brief Axis creator
- /// The top-level function taking part in the process of
- /// adding edges. Creating an axis is the same operation
- /// for it as adding new bins so it can be as well used to
- /// add some custom bins.
- ///
- /// It accepts two extremal points of a rectangle
- /// (top-right and bottom-left) as input.
- void _mkAxis(const vector<Segment>& binLimits) {
-
- /// For each of the rectangles
- for(unsigned int i=0; i < binLimits.size(); i++) {
- /// Produce the segments that a rectangle is composed of
- Segment edge1 =
- make_pair(binLimits[i].first,
- make_pair(binLimits[i].first.first, binLimits[i].second.second));
- Segment edge2 =
- make_pair(make_pair(binLimits[i].first.first, binLimits[i].second.second),
- binLimits[i].second);
- Segment edge3 =
- make_pair(make_pair(binLimits[i].second.first, binLimits[i].first.second),
- binLimits[i].second);
- Segment edge4 =
- make_pair(binLimits[i].first,
- make_pair(binLimits[i].second.first, binLimits[i].first.second));
-
- /// Check if they are made properly
- _fixOrientation(edge1); _fixOrientation(edge2);
- _fixOrientation(edge3); _fixOrientation(edge4);
-
- /// Add all the segments to a vector
- vector<Segment> edges;
- edges.push_back(edge1); edges.push_back(edge2);
- edges.push_back(edge3); edges.push_back(edge4);
-
- /// And check if a bin is a proper one, if it is, add it.
- if(_validateEdge(edges)) _addEdge(edges, _binHashSparse);
- else _dropEdge(edges);
-
- }
-
- /// Setting all the caches
- _binHashSparse.first.regenCache();
- _binHashSparse.second.regenCache();
- _regenDelimiters();
+ /// Get inf(Y)
+ const double lowEdgeY() const
+ {
+ return _lowEdgeY;
}
- /// @brief Plot extrema (re)generator.
- /// Since scrolling through every bin is an expensive
- /// operation to do every time we need the limits of
- /// the plot, there are caches set up. This function
- /// regenerates them. It should be run after any change is made
- /// to bin layout.
- void _regenDelimiters() {
- double highEdgeX = SMALLNUM;
- double highEdgeY = SMALLNUM;
- double lowEdgeX = LARGENUM;
- double lowEdgeY = LARGENUM;
+ ///Get sup(Y)
+ const double highEdgeY() const
+ {
+ return _highEdgeY;
+ }
- /// Scroll through the bins and set the delimiters.
- for(unsigned int i=0; i < _bins.size(); i++) {
- if(_bins[i].xMin() < lowEdgeX) lowEdgeX = _bins[i].xMin();
- if(_bins[i].xMax() > highEdgeX) highEdgeX = _bins[i].xMax();
- if(_bins[i].yMin() < lowEdgeY) lowEdgeY = _bins[i].yMin();
- if(_bins[i].yMax() > highEdgeY) highEdgeY = _bins[i].yMax();
- }
+ /// Get the bins from an Axis (non-const version)
+ Bins& bins()
+ {
+ return _bins;
+ }
- _lowEdgeX = lowEdgeX;
- _highEdgeX = highEdgeX;
- _lowEdgeY = lowEdgeY;
- _highEdgeY = highEdgeY;
+ /// Get the bins from an Axis (const version)
+ const Bins& bins() const
+ {
+ return _bins;
}
- /// @brief BIn index finder
- /// Looks through all the bins to see which
- /// one contains the point of interest.
- int _findBinIndex(double coordX, double coordY) const
+ /// Get a bin with a given index (non-const version)
+ BIN& bin(size_t index)
{
- for(size_t i=0; i < _bins.size(); i++) {
- if(_bins[i].lowEdgeX() <= coordX && _bins[i].highEdgeX() >= coordX &&
- _bins[i].lowEdgeY() <= coordY && _bins[i].highEdgeY() >= coordY) return i;
- }
- return -1;
+ if(index >= _bins.size()) throw RangeError("Bin index out of range.");
+ return _bins[index];
}
+ /// Get a bin with a given index (const version)
+ const BIN& bin(size_t index) const
+ {
+ if(index >= _bins.size()) throw RangeError("Bin index out of range.");
+ return _bins[index];
+ }
- public:
- /// @name Constructors:
- //@{
+ /// Get a bin at given coordinates (non-const version)
+ BIN& binByCoord(double x, double y)
+ {
+ int ret = _findBinIndex(x, y);
+ if(ret != -1) return bin(ret);
+ else throw RangeError("No bin found!!");
+ }
- /// @brief Empty constructor
- /// Only added because it is required by SWIG.
- /// It doesn't make much sense to use it.
- Axis2D()
+ /// Get a bin at given coordinates (const version)
+ const BIN& binByCoord(double x, double y) const
{
- vector<Segment> edges;
- _mkAxis(edges);
+ int ret = _findBinIndex(x, y);
+ if(ret != -1) return bin(ret);
+ else throw RangeError("No bin found!!");
}
- /// Constructor provided with a vector of bin delimiters
- Axis2D(const vector<Segment>& binLimits)
+ /// Get a bin at given coordinates (non-const version)
+ BIN& binByCoord(pair<double, double>& coords)
{
- _mkAxis(binLimits);
+ int ret = _findBinIndex(coords.first, coords.second);
+ if(ret != -1) return bin(ret);
+ else throw RangeError("No bin found!!");
}
- ///Most standard constructor, should be self-explanatory
- Axis2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY)
+ /// Get a bin at given coordinates (const version)
+ const BIN& binByCoord(pair<double, double>& coords) const
{
- vector<Segment> binLimits;
- double coeffX = (upperX - lowerX)/(double)nbinsX;
- double coeffY = (upperY - lowerX)/(double)nbinsY;
+ int ret = _findBinIndex(coords.first, coords.second);
+ if(ret != -1) return bin(ret);
+ else throw RangeError("No bin found!!");
+ }
- for(double i=lowerX; i < upperX; i+=coeffX) {
- for(double j=lowerY; j < upperY; j+=coeffY) {
- binLimits.push_back(make_pair(make_pair(i, j),
- make_pair((double)(i+coeffX), (double)(j+coeffY))));
- }
- }
- _mkAxis(binLimits);
+ /// Get a total distribution (non-const version)
+ Dbn2D& totalDbn()
+ {
+ return _dbn;
}
- //@}
- /// @name Addition operators:
- //@{
+ /// Get a total distribution (const version)
+ const Dbn2D& totalDbn() const
+ {
+ return _dbn;
+ }
- /// @brief Bin addition operator
- /// This operator is provided a vector of limiting
- /// points in the format required by _mkAxis().
- /// It should be noted that there is nothing special about
- /// the initiation stage of Axis2D, and the edges can be added
- /// online if they meet all the requirements of non-degeneracy.
- /// No merging is supported, and I don't think it should before the support
- /// for merging for '+' operator (codewise it should be the same thing).
- void addBin(const vector<Segment>& binLimits)
+ /// Get the overflow distribution (non-const version)
+ Dbn2D& overflow()
{
- _mkAxis(binLimits);
+ return _overflow;
}
- /// @brief Bin addition operator
- /// This operator is supplied with whe extreamal coordinates of just
- /// one bin. It then launches the standard bin addition procedure.
- void addBin(double lowX, double lowY, double highX, double highY)
+ /// Get the overflow distribution (const version)
+ const Dbn2D& overflow() const
{
- vector<Segment> coords;
- coords.push_back(make_pair(make_pair(lowX, lowY), make_pair(highX, highY)));
+ return _overflow;
+ }
- addBin(coords);
- }
- //@}
-
- /// @name Some helper functions:
- //@{
-
- /// @brief Checks if our bins form a grid.
- /// This function uses a neat property of _binHashSparse.
- /// If it is containing a set of edges forming a grid without
- /// gaps in the middle it will have the same number of edges in the
- /// inner subcaches and half of this amount in the outer (grid boundary)
- /// subcaches. This makes isGriddy() a very, very fast function.
- /// @todo Is the name appropriate?
- int isGriddy() const
- {
-
- /// Check if the number of edges parallel to X axis
- /// is proper in every subcache.
- unsigned int sizeX = _binHashSparse.first[0].second.size();
- for(unsigned int i=1; i < _binHashSparse.first.size(); i++) {
- if(i == _binHashSparse.first.size() - 1) {
- if(_binHashSparse.first[i].second.size() != sizeX) {
- return -1;
- }
- }
- else if(_binHashSparse.first[i].second.size() != 2*sizeX) {
- return -1;
- }
- }
-
- /// Do the same for edges parallel to Y axis.
- unsigned int sizeY = _binHashSparse.second[0].second.size();
- for(unsigned int i=1; i < _binHashSparse.second.size(); i++) {
- if(i!= _binHashSparse.second.size() - 1) {
- if(2*sizeY != _binHashSparse.second[i].second.size()) {
- return -1;
- }
- }
- else if(_binHashSparse.second[i].second.size() != sizeY) return -1;
- }
-
- /// If everything is proper, announce it.
- return 0;
- }
-
- /// Return a total number of bins in a Histo
- unsigned int numBinsTotal() const
- {
- return _bins.size();
- }
-
- /// Get inf(X) (non-const version)
- double lowEdgeX()
- {
- return _lowEdgeX;
- }
-
- /// Get sup(X) (non-const version)
- double highEdgeX()
- {
- return _highEdgeX;
- }
-
- /// Get inf(Y) (non-const version)
- double lowEdgeY()
- {
- return _lowEdgeY;
- }
-
- /// Get sup(Y) (non-const version)
- double highEdgeY()
- {
- return _highEdgeY;
- }
-
- /// Get inf(X) (const version)
- const double lowEdgeX() const
- {
- return _lowEdgeX;
- }
-
- /// Get sup(X) (const version)
- const double highEdgeX() const
- {
- return _highEdgeX;
- }
-
- /// Get inf(Y)
- const double lowEdgeY() const
- {
- return _lowEdgeY;
- }
-
- ///Get sup(Y)
- const double highEdgeY() const
- {
- return _highEdgeY;
- }
-
- /// Get the bins from an Axis (non-const version)
- Bins& bins()
- {
- return _bins;
- }
-
- /// Get the bins from an Axis (const version)
- const Bins& bins() const
- {
- return _bins;
- }
-
- /// Get a bin with a given index (non-const version)
- BIN& bin(size_t index)
- {
- /// @todo WTF?! Don't just copy and paste exception warnings. This one is totally misleading: are there more?
- if(index >= _bins.size()) throw RangeError("Bin index out of range.");
- return _bins[index];
- }
-
- /// Get a bin with a given index (const version)
- const BIN& bin(size_t index) const
- {
- /// @todo WTF?! Don't just copy and paste exception warnings. This one is totally misleading: are there more?
- if(index >= _bins.size()) throw RangeError("Bin index out of range.");
- return _bins[index];
- }
-
- /// Get a bin at given coordinates (non-const version)
- BIN& binByCoord(double x, double y)
- {
- int ret = _findBinIndex(x, y);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
- }
-
- /// Get a bin at given coordinates (const version)
- const BIN& binByCoord(double x, double y) const
- {
- int ret = _findBinIndex(x, y);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
- }
-
- /// Get a bin at given coordinates (non-const version)
- BIN& binByCoord(pair<double, double>& coords)
- {
- int ret = _findBinIndex(coords.first, coords.second);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
- }
-
- /// Get a bin at given coordinates (const version)
- const BIN& binByCoord(pair<double, double>& coords) const
- {
- int ret = _findBinIndex(coords.first, coords.second);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
- }
-
- /// Get a total distribution (non-const version)
- Dbn2D& totalDbn()
- {
- return _dbn;
- }
-
- /// Get a total distribution (const version)
- const Dbn2D& totalDbn() const
- {
- return _dbn;
- }
-
- /// Get the overflow distribution (non-const version)
- Dbn2D& overflow()
- {
- return _overflow;
- }
-
- /// Get the overflow distribution (const version)
- const Dbn2D& overflow() const
- {
- return _overflow;
- }
-
- /// Get the underflow distribution (non-const version)
- Dbn2D& underflow()
- {
- return _underflow;
+ /// Get the underflow distribution (non-const version)
+ Dbn2D& underflow()
+ {
+ return _underflow;
}
/// Get the underflow distribution (const version)
@@ -649,18 +342,8 @@
return _underflow;
}
- /// Get bin index from external classes (non-const version)
- /// @todo Change: the YODA method naming convention doesn't have "get"s
- /// @todo Why isn't *this* method const?
- int getBinIndex(double coordX, double coordY)
- {
- return _findBinIndex(coordX, coordY);
- }
-
/// Get bin index from external classes (const version)
- /// @todo Change: the YODA method naming convention doesn't have "get"s
- /// @todo This method is completely pointless: I think that method constancy has been misunderstood :(
- const int getBinIndex(double coordX, double coordY) const
+ int getBinIndex(double coordX, double coordY) const
{
return _findBinIndex(coordX, coordY);
}
@@ -726,7 +409,7 @@
//@{
/// Equality operator
- bool operator == (const Axis2D& other) const
+ bool operator == (const Axis2D& other) const
{
return _binHashSparse == other._binHashSparse;
}
@@ -768,6 +451,330 @@
}
//@}
+
+ private:
+
+ /// @brief Segment validator function
+ /// This a 'dispatcher' function. It checks if the segment in question
+ /// is vertical or horizontal and launches the appropriate function
+ /// searching for cuts in the prope direction. Since it operates on
+ /// a vector of segments it is prepared to act on arbitrarly large sets of edges,
+ /// in practice usually being four sides of a rectangular bin.
+ ///
+ /// Notice that it will never be checked, in the current state, if there is a cut
+ /// in edges in the edgeset. This imposes the requirement of provide the program
+ /// with non-degenerate bins, until checking in implemented. However, it doesn't seem
+ /// to be needed, as edges are not generated by a user.
+ ///
+ /// This is also a perfect place to paralellize the program, if required.
+ bool _validateEdge(vector<Segment>& edges)
+ {
+ /// Setting the return variable. True means that no cuts were detected.
+ bool ret = true;
+
+ /// Looping over all the edges provided
+ for(unsigned int i=0; i < edges.size(); i++) {
+ /// If the X coordinate of the starting point is the same
+ /// as X coordinate of the ending one, checks if there are cuts
+ /// on this vertical segment.
+ if(fuzzyEquals(edges[i].first.first, edges[i].second.first)) ret = _findCutsY(edges[i]);
+
+ /// Check if the segment is horizontal and is it cutting any bin that already exists
+ else if(fuzzyEquals(edges[i].first.second, edges[i].second.second)) ret = _findCutsX(edges[i]);
+
+ /// This is a check that discards the bin if it is not a rectangle
+ /// composed of vertical and horizontal segments.
+ else ret = false;
+
+ /// If a cut was detected, say it. There is no point in checking other edges
+ /// in the set.
+ if(!ret) return false;
+ }
+ /// If no cuts were detected in any of the edges, tell the launching function about this
+ return true;
+ }
+
+ /// @brief A binary search function
+ /// This is conceptually the same implementation as in STL
+ /// but because it returns the index of an element in a vector,
+ /// it is easier to use in our case than the STL implementation
+ /// that returns a pointer at the element.
+ size_t _binaryS(Utils::cachedvector<EdgeCollection>& toSearch,
+ double value, size_t lower, size_t higher)
+ {
+ /// Just a check if such degenerate situation happens
+ if(lower == higher) return lower;
+
+ /// Choose a midpoint that will be our pivot
+ size_t where = (higher+lower)/2;
+
+ /// Launch the same procedure on half of the range above the pivot
+ if(value >= toSearch[where].first) {
+ if(where == toSearch.size() - 1) return where;
+ if(value <= toSearch[where+1].first) return where;
+ return _binaryS(toSearch, value, where, higher);
+ }
+
+ /// This is not a redundant check, because
+ /// of the nature of int division.
+ if (where == 0) return where;
+
+ /// Check if the value is somewhere inbetween
+ /// an element at the position in question and
+ /// an element at a lower position. If so, return
+ /// an index to the current positon.
+ if(value >= toSearch[where-1].first) return where;
+
+ /// If none of the above occurs, the value must
+ /// be smaller that the element at the current
+ /// position. In such case, launch the search on
+ /// half of the interval below the current position.
+ return _binaryS(toSearch, value, lower, where);
+ }
+
+ /// @brief Function that finds cuts of horizontal edges.
+ /// A specialised function that tries to exploit the fact
+ /// that edges can cut one another only on right angles
+ /// to the highest extent. The inner workings are explained
+ /// in the comments placed in the function body.
+ bool _findCutsX(Segment& edge) {
+ /// Look up the limits of search in the _binHashSparse
+ /// structure. We are not interested in the vertical edges
+ /// that are before the beginning of our segment or after its end.
+ size_t i = _binaryS(_binHashSparse.second, edge.first.first, 0, _binHashSparse.second.size());
+ size_t end = _binaryS(_binHashSparse.second, edge.second.first, 0, _binHashSparse.second.size());
+
+ for(; i < end; i++) {
+
+ /// Scroll through all the vertical segments with a given X coordinate
+ /// and see if any of those fulfills the cutting requirement. If it does,
+ /// announce it.
+ for(unsigned int j = 0; j < _binHashSparse.second[i].second.size(); j++) {
+ /// Note that we are not taking into account the edges touching the
+ /// segment in question. That's because sides of a bin touch
+ if(_binHashSparse.second[i].second[j].second.first < edge.first.second &&
+ _binHashSparse.second[i].second[j].second.second > edge.first.second &&
+ !fuzzyEquals(_binHashSparse.second[i].second[j].second.first, edge.first.second) &&
+ !fuzzyEquals(_binHashSparse.second[i].second[j].second.second, edge.first.second)) {
+ return false;
+ }
+ }
+ }
+ /// If none of the existing edges is cutting this edge, announce it
+ return true;
+ }
+
+ /// @brief Function that finds cuts of vertical edges
+ /// For a detailed descpription, please look into
+ /// documentation for _findCutsX().
+ bool _findCutsY(Segment& edge) {
+ size_t i = _binaryS(_binHashSparse.first, edge.first.second, 0, _binHashSparse.first.size());
+ size_t end = _binaryS(_binHashSparse.first, edge.second.second, 0, _binHashSparse.first.size());
+ for(; i < end; i++) {
+
+ for(unsigned int j = 0; j < _binHashSparse.first[i].second.size(); j++) {
+ if(_binHashSparse.first[i].second[j].second.first < edge.first.first &&
+ _binHashSparse.first[i].second[j].second.second > edge.first.first &&
+ !fuzzyEquals(_binHashSparse.first[i].second[j].second.first, edge.first.first) &&
+ !fuzzyEquals(_binHashSparse.first[i].second[j].second.second, edge.first.first)) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+ /// @brief Function executed when a set of edges is dropped.
+ /// It does not have any information about which edge in the set
+ /// had failed the check. If this is needed such additional information
+ /// can be readily implemented.
+ void _dropEdge(vector<Segment>& edges) {
+ std::cerr << "A set of edges was dropped." << endl;
+ }
+
+ /// @brief Bin adder.
+ /// It contains all the commands that need to executed
+ /// to properly add a bin. Specifially edges are added to
+ /// the edge cache (_binHashSparse) and a bin is created from
+ /// those edges.
+ void _addEdge(vector<Segment> edges, pair<Utils::cachedvector<EdgeCollection>,
+ Utils::cachedvector<EdgeCollection> >& binHash, bool addBin = true) {
+ /// Check if there was no mistake made when adding segments to a vector.
+ if(edges.size() != 4) throw Exception("The segments supplied don't describe a full bin!");
+
+ /// This is the part in charge of adding each of the segments
+ /// to the edge cache. Segments are assumed to be validated
+ /// beforehand.
+ for(unsigned int j=0; j < edges.size(); j++) {
+ /// Association made for convinience.
+ Segment edge = edges[j];
+
+ /// Do the following if the edge is vertical.
+ /// Those two cases need to be distinguished
+ /// because of the way in which edge cache is structured.
+ if(edge.first.first == edge.second.first) {
+ /// See if our edge has the same X coordinate as any other
+ /// edge that is currently in the cache.
+
+ /// Keeps the status of the search
+ bool found = false;
+
+ /// There is only a certain set of X coordinates that we need to sweep
+ /// to check if our segment has the same X coordinate. Find them.
+ size_t i = _binaryS(binHash.second, edge.first.first, 0, binHash.second.size())-1;
+ if(i < 0) i = 0;
+ size_t end = i+3;
+
+ /// For the coordinates in range, check if one of them is an X coordinate of
+ /// the sement.
+ for(; i < binHash.second.size() && i < end ; i++) {
+ /// If this is the case, do what is needed to be done.
+ if(fuzzyEquals(binHash.second[i].first, edge.first.first)) {
+ binHash.second[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.second, edge.second.second)));
+ found = true;
+ break;
+ }
+ }
+
+ /// If no edge with the same X coordinate exist, create
+ /// a new subhash at the X coordinate of a segment.
+ if(!found) {
+ vector<Edge> temp;
+ temp.push_back(make_pair(_bins.size(), make_pair(edge.first.second, edge.second.second)));
+ binHash.second.push_back(make_pair(edge.first.first,temp));
+ sort(binHash.second.begin(), binHash.second.end());
+ }
+ }
+
+ /// See the vertical case for description of a horizontal one
+ else if(edge.first.second == edge.second.second) {
+ bool found = false;
+ size_t i = _binaryS(binHash.first, edge.first.second, 0, binHash.first.size())-1;
+ if(i < 0) i = 0;
+ size_t end = i+3;
+ for(; i < binHash.first.size() && i < end; i++) {
+ if(fuzzyEquals(binHash.first[i].first, edge.first.second)) {
+ binHash.first[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.first, edge.second.first)));
+ found = true;
+ }
+ }
+ if(!found) {
+ vector<Edge> temp;
+ temp.push_back(make_pair(_bins.size(), make_pair(edge.first.first, edge.second.first)));
+ binHash.first.push_back(make_pair(edge.second.second, temp));
+ sort(binHash.first.begin(), binHash.first.end());
+ }
+ }
+ }
+ /// Now, create a bin with the edges provided
+ if(addBin) _bins.push_back(BIN(edges));
+ }
+
+ /// @brief Orientation fixer
+ /// Check if the orientation of an edge is proper
+ /// for the rest of the algorithm to work on, and if it is not
+ /// fix it.
+ void _fixOrientation(Segment& edge) {
+ if(fuzzyEquals(edge.first.first, edge.second.first)) {
+ if(edge.first.second > edge.second.second) {
+ double temp = edge.second.second;
+ edge.second.second = edge.first.second;
+ edge.first.second = temp;
+ }
+ }
+ else if(edge.first.first > edge.second.first) {
+ double temp = edge.first.first;
+ edge.first.first = edge.second.first;
+ edge.second.first = temp;
+ }
+ }
+
+ /// @brief Axis creator
+ /// The top-level function taking part in the process of
+ /// adding edges. Creating an axis is the same operation
+ /// for it as adding new bins so it can be as well used to
+ /// add some custom bins.
+ ///
+ /// It accepts two extremal points of a rectangle
+ /// (top-right and bottom-left) as input.
+ void _mkAxis(const vector<Segment>& binLimits) {
+
+ /// For each of the rectangles
+ for(unsigned int i=0; i < binLimits.size(); i++) {
+ /// Produce the segments that a rectangle is composed of
+ Segment edge1 =
+ make_pair(binLimits[i].first,
+ make_pair(binLimits[i].first.first, binLimits[i].second.second));
+ Segment edge2 =
+ make_pair(make_pair(binLimits[i].first.first, binLimits[i].second.second),
+ binLimits[i].second);
+ Segment edge3 =
+ make_pair(make_pair(binLimits[i].second.first, binLimits[i].first.second),
+ binLimits[i].second);
+ Segment edge4 =
+ make_pair(binLimits[i].first,
+ make_pair(binLimits[i].second.first, binLimits[i].first.second));
+
+ /// Check if they are made properly
+ _fixOrientation(edge1); _fixOrientation(edge2);
+ _fixOrientation(edge3); _fixOrientation(edge4);
+
+ /// Add all the segments to a vector
+ vector<Segment> edges;
+ edges.push_back(edge1); edges.push_back(edge2);
+ edges.push_back(edge3); edges.push_back(edge4);
+
+ /// And check if a bin is a proper one, if it is, add it.
+ if(_validateEdge(edges)) _addEdge(edges, _binHashSparse);
+ else _dropEdge(edges);
+
+ }
+
+ /// Setting all the caches
+ _binHashSparse.first.regenCache();
+ _binHashSparse.second.regenCache();
+ _regenDelimiters();
+ }
+
+ /// @brief Plot extrema (re)generator.
+ /// Since scrolling through every bin is an expensive
+ /// operation to do every time we need the limits of
+ /// the plot, there are caches set up. This function
+ /// regenerates them. It should be run after any change is made
+ /// to bin layout.
+ void _regenDelimiters() {
+ double highEdgeX = SMALLNUM;
+ double highEdgeY = SMALLNUM;
+ double lowEdgeX = LARGENUM;
+ double lowEdgeY = LARGENUM;
+
+ /// Scroll through the bins and set the delimiters.
+ for(unsigned int i=0; i < _bins.size(); i++) {
+ if(_bins[i].xMin() < lowEdgeX) lowEdgeX = _bins[i].xMin();
+ if(_bins[i].xMax() > highEdgeX) highEdgeX = _bins[i].xMax();
+ if(_bins[i].yMin() < lowEdgeY) lowEdgeY = _bins[i].yMin();
+ if(_bins[i].yMax() > highEdgeY) highEdgeY = _bins[i].yMax();
+ }
+
+ _lowEdgeX = lowEdgeX;
+ _highEdgeX = highEdgeX;
+ _lowEdgeY = lowEdgeY;
+ _highEdgeY = highEdgeY;
+ }
+
+ /// @brief BIn index finder
+ /// Looks through all the bins to see which
+ /// one contains the point of interest.
+ int _findBinIndex(double coordX, double coordY) const
+ {
+ for(size_t i=0; i < _bins.size(); i++) {
+ if(_bins[i].lowEdgeX() <= coordX && _bins[i].highEdgeX() >= coordX &&
+ _bins[i].lowEdgeY() <= coordY && _bins[i].highEdgeY() >= coordY &&
+ _bins[i].isReal) return i;
+ }
+ return -1;
+ }
+
private:
/// Bins contained in this histogram
Modified: trunk/include/YODA/Bin2D.h
==============================================================================
--- trunk/include/YODA/Bin2D.h Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/include/YODA/Bin2D.h Tue Aug 16 10:35:15 2011 (r248)
@@ -6,12 +6,12 @@
#include <string>
#include <utility>
#include <vector>
+#include <iostream>
#include <cassert>
using namespace std;
namespace YODA {
-
-
+
/// @brief A generic 2D bin type
///
/// This is a generic 2D bin type which supplies the accessors for the two "x"
@@ -23,6 +23,10 @@
class Bin2D : public Bin {
public:
+ /// Convinience typedefs
+ typedef typename std::pair<double, double> Point;
+ typedef typename std::pair<Point, Point> Segment;
+
/// @name Constructors
//@{
@@ -37,7 +41,7 @@
/// Bin slightly faster (this claim is very weakly true). It is not
/// suggested to use it if it is just needed to add few bins to an already
/// created Histo2D.
- Bin2D(std::vector<std::pair<std::pair<double, double>, std::pair<double, double> > > edges);
+ Bin2D(std::vector<Segment> edges);
//@}
@@ -45,6 +49,8 @@
/// @name Modifiers
//@{
+ const vector<Segment> edges() const { return _edges;}
+
/// Reset all bin data
virtual void reset() {
_dbn.reset();
@@ -70,9 +76,15 @@
double lowEdgeY() const {
return _edges[0].first.second;
}
+
/// Synonym for lowEdgeY()
double yMin() const { return lowEdgeY(); }
+ /// A variable that specifies if the bin should be plotted
+ bool isReal;
+
+ ///@name Transformers
+ //@{
/// Get the high x edge of the bin.
double highEdgeX() const {
return _edges[1].second.first;
@@ -87,34 +99,29 @@
/// Synonym for highEdgeY()
double yMax() const { return highEdgeY(); }
+ /// Width of the bin in y
+ double widthY() const {
+ return _edges[0].second.second - _edges[0].first.second;
+ }
+
/// Width of the bin in x
double widthX() const {
return _edges[1].second.first - _edges[0].first.first;
}
+ //@}
- /// Width of the bin in y
- double widthY() const {
- return _edges[0].second.second - _edges[0].first.second;
- }
+ /// @name Distribution statistics
+ //@{
/// Find the geometric midpoint of the bin
- std::pair<double, double> midpoint() const {
- return make_pair(_edges[1].second.first-_edges[0].first.first,
- _edges[0].second.second-_edges[0].first.second);
- }
+ Point midpoint() const;
/// Find the weighted mean point of the bin, or the midpoint if unfilled
- std::pair<double, double> focus() const {
+ Point focus() const {
if (_dbn.sumW() != 0) return make_pair(xMean(), yMean());
return midpoint();
}
- //@}
-
-
- /// @name Distribution statistics
- //@{
-
/// Mean x value
double xMean() const {
return _dbn.xMean();
@@ -206,34 +213,25 @@
//@}
-
protected:
- Bin2D& add(const Bin2D& b) {
- assert(_edges == b._edges);
- _dbn += b._dbn;
- return *this;
- }
-
- Bin2D& subtract(const Bin2D& b) {
- assert(_edges == b._edges);
- _dbn -= b._dbn;
- return *this;
- }
+ Bin2D& add(const Bin2D& b);
+ Bin2D& subtract(const Bin2D& b);
- protected:
-
-
- /// @todo Explain! And then replace it with something sane -- this is really wasteful.
- std::vector<std::pair<std::pair<double,double>,std::pair<double,double> > > _edges;
-
- /// Distribution of fills in this bin.
+ std::vector<Segment> _edges;
Dbn2D _dbn;
+
+ /// Boundaries setter
+ void _setBounds(double xMin, double yMin, double xMax, double yMax) {
+ _edges[0] = make_pair(make_pair(xMin, yMin), make_pair(xMin, yMax));
+ _edges[1] = make_pair(make_pair(xMin, yMax), make_pair(xMax, yMax));
+ _edges[2] = make_pair(make_pair(xMax, yMin), make_pair(xMax, yMax));
+ _edges[3] = make_pair(make_pair(xMin, yMin), make_pair(xMax, yMin));
+ }
};
-
inline Bin2D operator + (const Bin2D& a, const Bin2D& b) {
Bin2D rtn = a;
rtn += b;
Modified: trunk/include/YODA/Histo2D.h
==============================================================================
--- trunk/include/YODA/Histo2D.h Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/include/YODA/Histo2D.h Tue Aug 16 10:35:15 2011 (r248)
@@ -10,7 +10,6 @@
#include "YODA/HistoBin2D.h"
#include "YODA/HistoBin1D.h"
#include "YODA/Axis2D.h"
-#include "YODA/Scatter3D.h"
#include "YODA/Exceptions.h"
#include "YODA/Histo1D.h"
#include <vector>
@@ -19,6 +18,8 @@
namespace YODA {
+ // Forward declaration
+ class Scatter3D;
/// Convenience typedef
typedef Axis2D<HistoBin2D> Histo2DAxis;
@@ -114,7 +115,11 @@
void addBin(double lowX, double lowY, double highX, double highY) {
_axis.addBin(lowX, lowY, highX, highY);
}
-
+
+ /// Merge the bins
+ void mergeBins(size_t from, size_t to) {
+ _axis.mergeBins(from, to);
+ }
//@}
public:
@@ -210,19 +215,6 @@
return _axis.numBinsTotal();
}
- /// Hash returner (non-const version)
- /// @todo This needs a typedef
- /* std::pair<Utils::cachedvector<pair<double,std::vector<pair<size_t, pair<double,double> > > > >,
- Utils::cachedvector<pair<double,std::vector<pair<size_t, pair<double,double> > > > > > getHash() {
- return _axis.getHash();
- }
-
- /// Hash returner (const version)
- /// @todo This needs a typedef
- const std::pair<Utils::cachedvector<pair<double,std::vector<pair<size_t, pair<double,double> > > > >,
- Utils::cachedvector<pair<double,std::vector<pair<size_t, pair<double,double> > > > > > getHash() const {
- return _axis.getHash();
- }*/
//@}
public:
@@ -278,14 +270,14 @@
return *this;
}
- /*bool operator == (const Histo2D& other) {
+ bool operator == (const Histo2D& other) const {
return _axis == other._axis;
}
- bool operator != (const Histo2D& other) {
+ bool operator != (const Histo2D& other) const {
return ! operator == (other);
}
-*/
+
//@}
@@ -300,18 +292,22 @@
Histo1D cutterX(double atY, const std::string& path="", const std::string& title="") {
if (atY < lowEdgeY() || atY > highEdgeY()) throw RangeError("Y is outside the grid");
vector<HistoBin1D> tempBins;
+
/// @todo Make all Bin1D constructions happen in loop, to reduce code duplication
const HistoBin2D& first = binByCoord(lowEdgeX(), atY);
const Dbn1D dbn(first.numEntries(), first.sumW(), first.sumW2(), first.sumWX(), first.sumWX2());
tempBins.push_back(HistoBin1D(first.lowEdgeX(), first.highEdgeX(), dbn));
+
for (double i = first.xMax() + first.widthX()/2; i < highEdgeX(); i += first.widthX()) {
const HistoBin2D& b2 = binByCoord(i, atY);
const Dbn1D dbn2(b2.numEntries(), b2.sumW(), b2.sumW2(), b2.sumWX(), b2.sumWX2());
tempBins.push_back(HistoBin1D(b2.lowEdgeX(), b2.highEdgeX(), dbn2));
}
+
/// @todo Think about the total, underflow and overflow distributions
/// @todo Create total dbn from input bins
return Histo1D(tempBins, Dbn1D(), Dbn1D(), Dbn1D(), path, title);
+
}
@@ -323,15 +319,18 @@
Histo1D cutterY(double atX, const std::string& path="", const std::string& title="") {
if (atX < lowEdgeX() || atX > highEdgeX()) throw RangeError("X is outside the grid");
vector<HistoBin1D> tempBins;
+
/// @todo Make all Bin1D constructions happen in loop, to reduce code duplication
const HistoBin2D& first = binByCoord(atX, lowEdgeY());
const Dbn1D dbn(first.numEntries(), first.sumW(), first.sumW2(), first.sumWX(), first.sumWX2());
tempBins.push_back(HistoBin1D(first.lowEdgeY(), first.highEdgeY(), dbn));
+
for (double i = first.yMax() + first.widthY()/2; i < highEdgeY(); i += first.widthY()) {
const HistoBin2D& b2 = binByCoord(atX, i);
const Dbn1D dbn2(b2.numEntries(), b2.sumW(), b2.sumW2(), b2.sumWX(), b2.sumWX2());
tempBins.push_back(HistoBin1D(b2.lowEdgeY(), b2.highEdgeY(), dbn2));
}
+
/// @todo Think about the total, underflow and overflow distributions
/// @todo Create total dbn from input bins
return Histo1D(tempBins, Dbn1D(), Dbn1D(), Dbn1D(), path, title);
@@ -375,7 +374,7 @@
return tmp;
}
- //Scatter3D operator / (const Histo2D& numer, const Histo2D& denom);
+ Scatter3D operator / (const Histo2D& numer, const Histo2D& denom);
//@}
Modified: trunk/include/YODA/Scatter3D.h
==============================================================================
--- trunk/include/YODA/Scatter3D.h Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/include/YODA/Scatter3D.h Tue Aug 16 10:35:15 2011 (r248)
@@ -84,7 +84,8 @@
addPoint(Point3D(x[i], exminus[i], explus[i], y[i], eyminus[i], eyplus[i], z[i], ezminus[i], ezplus[i]));
}
}
-
+ /// Copy constructor
+ Scatter3D(const Scatter3D&, const std::string& path="");
//@}
Modified: trunk/src/Bin2D.cc
==============================================================================
--- trunk/src/Bin2D.cc Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/src/Bin2D.cc Tue Aug 16 10:35:15 2011 (r248)
@@ -10,10 +10,9 @@
namespace YODA {
-
Bin2D::Bin2D(double lowedgeX, double lowedgeY, double highedgeX, double highedgeY) {
assert(lowedgeX <= highedgeX && lowedgeY <= highedgeY);
-
+
/// @todo Why store with so much redundancy?
pair<pair<double,double>, pair<double,double> > edge1 =
make_pair(make_pair(lowedgeX, lowedgeY), make_pair(lowedgeX, highedgeY));
@@ -28,8 +27,9 @@
_edges.push_back(edge2);
_edges.push_back(edge3);
_edges.push_back(edge4);
- }
+ isReal = true;
+ }
Bin2D::Bin2D(std::vector<std::pair<std::pair<double,double>,
std::pair<double,double> > > edges) {
@@ -38,6 +38,8 @@
_edges.push_back(edges[1]);
_edges.push_back(edges[2]);
_edges.push_back(edges[3]);
+
+ isReal = true;
}
@@ -51,5 +53,31 @@
_dbn.scaleXY(scaleX, scaleY);
}
+ std::pair<double,double> Bin2D::midpoint() const {
+ return make_pair((double)(xMax() - xMin())/2 + xMin(), (double)(yMax() - yMin())/2 + yMin());
+ }
+
+ Bin2D& Bin2D::subtract(const Bin2D& b) {
+ /// Automatically resize if adding a bin that does not have the same location
+ /// this way merging the bins works perfectly
+ if(_edges != b._edges) {
+ if (b.xMax() > xMax()) _setBounds(xMin(), yMin(), b.xMax(), yMax());
+ if (b.yMax() > yMax()) _setBounds(xMin(), yMin(), xMax(), b.yMax());
+ if (b.xMin() < xMin()) _setBounds(b.xMin(), yMin(), xMax(), yMax());
+ if (b.yMin() < yMin()) _setBounds(xMin(), b.yMin(), xMax(), yMax());
+ }
+ _dbn -= b._dbn;
+ return *this;
+ }
+ Bin2D& Bin2D::add(const Bin2D& b) {
+ if(_edges != b._edges) {
+ if (b.highEdgeX() > highEdgeX());// _setBounds(xMin(), yMin(), b.xMax(), yMax());
+ if (b.yMax() > yMax()) _setBounds(xMin(), yMin(), xMax(), b.yMax());
+ if (b.xMin() < xMin()) _setBounds(b.xMin(), yMin(), xMax(), yMax());
+ if (b.yMin() < yMin()) _setBounds(xMin(), b.yMin(), xMax(), yMax());
+ }
+ _dbn += b._dbn;
+ return *this;
+ }
}
Modified: trunk/src/Histo1D.cc
==============================================================================
--- trunk/src/Histo1D.cc Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/src/Histo1D.cc Tue Aug 16 10:35:15 2011 (r248)
@@ -117,8 +117,11 @@
for (size_t i = 0; i < numer.numBins(); ++i) {
const HistoBin1D& b1 = numer.bin(i);
const HistoBin1D& b2 = denom.bin(i);
- assert(fuzzyEquals(b1.focus(), b2.focus()));
- const double x = b1.focus();
+ const HistoBin1D& bL = b1 + b2;
+
+ assert(fuzzyEquals(b1.midpoint(), b2.midpoint()));
+
+ const double x = bL.focus();
assert(fuzzyEquals(b1.xMin(), b2.xMin()));
assert(fuzzyEquals(b1.xMax(), b2.xMax()));
const double exminus = x - b1.xMin();
Modified: trunk/src/Histo2D.cc
==============================================================================
--- trunk/src/Histo2D.cc Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/src/Histo2D.cc Tue Aug 16 10:35:15 2011 (r248)
@@ -16,20 +16,21 @@
typedef vector<HistoBin2D> Bins;
- //TODO: It should return/throw something if no bin exist.
- //So, what is the concept of under/over flow in 2D case?
+ /// @brief Fill function
+ /// It relies on Axis2D for searching procedures and
+ /// complies loudly if no bin was found.
int Histo2D::fill(double x, double y, double weight) {
- // Fill the normal bins
+ /// Fill the normal bins
int index = findBinIndex(x, y);
- //cout << index << endl;
if(index != -1) {
HistoBin2D& b = bin(index);
- // Fill the underflow and overflow nicely
+ /// Fill the underflow and overflow nicely
_axis.totalDbn().fill(x, y, weight);
b.fill(x, y, weight);
}
+
else if (x < _axis.lowEdgeX()) { _axis.underflow().fill(x, y, weight); }
else if (x >= _axis.highEdgeX()) { _axis.overflow().fill(x, y, weight); }
else throw GridError("You are trying to fill an empty space on a grid!");
@@ -114,35 +115,27 @@
///////////////////////////////////////
- /// Divide two histograms
- /* Scatter3D operator / (Histo2D& numer, const Histo2D& denom) {
- if(numer!=denom) throw GridError("The two Histos are not the same!");
+ /// @brief Divide two histograms
+ /// This uses the midpoint instead of the focus
+ Scatter3D operator / (const Histo2D& numer, const Histo2D& denom) {
+ if(numer != denom) throw GridError("The two Histos are not the same!");
Scatter3D tmp;
for (size_t i = 0; i < numer.numBinsTotal(); ++i) {
const HistoBin2D& b1 = numer.bin(i);
const HistoBin2D& b2 = denom.bin(i);
- assert(fuzzyEquals(b1.focus().first, b2.focus().first));
- assert(fuzzyEquals(b1.focus().second, b2.focus().second));
+ const HistoBin2D& bL = b1 + b2;
+ assert(fuzzyEquals(b1.midpoint().first, b2.midpoint().first));
+ assert(fuzzyEquals(b1.midpoint().second, b2.midpoint().second));
- const double x = b1.focus().first;
- const double y = b1.focus().second;
-
- assert(fuzzyEquals(b1.xMin(), b2.xMin()));
- assert(fuzzyEquals(b1.xMax(), b2.xMax()));
-
- assert(fuzzyEquals(b1.yMin(), b2.yMin()));
- assert(fuzzyEquals(b1.yMax(), b2.yMax()));
+ const double x = bL.focus().first;
+ const double y = bL.focus().second;
- const double exminus = x - b1.xMin();
- const double explus = b1.xMax() - x;
+ const double exminus = x - bL.xMin();
+ const double explus = bL.xMax() - x;
const double eyminus = y - b1.yMin();
- const double eyplus = b1.yMax() - y;
-
- assert(exminus >= 0);
- assert(explus >= 0);
- assert(eyminus >= 0);
- assert(eyplus >= 0);
+ const double eyplus = bL.yMax() - y;
+ //cout << b1.xMin() << " " << b1.xMax() << " " << b1.yMin() << " " << b1.yMax() << " EMinus: " << exminus << " " << eyminus << " Focus: " << x << " " << y << endl;
const double z = b1.height() / b2.height();
const double ez = z * sqrt( sqr(b1.heightErr()/b1.height()) + sqr(b2.heightErr()/b2.height()) );
@@ -151,5 +144,5 @@
assert(tmp.numPoints() == numer.numBinsTotal());
return tmp;
}
-*/
+
}
Modified: trunk/src/Scatter3D.cc
==============================================================================
--- trunk/src/Scatter3D.cc Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/src/Scatter3D.cc Tue Aug 16 10:35:15 2011 (r248)
@@ -29,6 +29,12 @@
*/
////////////////////////////////////////
+ /// A copy constructor:
+ Scatter3D::Scatter3D(const Scatter3D& s, const std::string& path)
+ : AnalysisObject("Scatter3D", (path.size() == 0) ? s.path() : path, s, s.title())
+ {
+ _points = s._points;
+ }
/// Subtract two scatters
Scatter3D operator + (const Scatter3D& first, const Scatter3D& second) {
Modified: trunk/tests/TestHisto2D.cc
==============================================================================
--- trunk/tests/TestHisto2D.cc Mon Aug 15 19:33:59 2011 (r247)
+++ trunk/tests/TestHisto2D.cc Tue Aug 16 10:35:15 2011 (r248)
@@ -1,4 +1,5 @@
#include "YODA/Histo2D.h"
+#include "YODA/Scatter3D.h"
#include <cmath>
#include <iostream>
#include <unistd.h>
@@ -38,12 +39,15 @@
gettimeofday(&startTime, NULL);
Histo2D h(200, 0, 100, 200, 0, 100);
gettimeofday(&endTime, NULL);
+ cout << h.binByCoord(1,1).isReal << endl;
double tS = (startTime.tv_sec*1000000 + startTime.tv_usec)/(double)1000000;
double tE = (endTime.tv_sec*1000000 + endTime.tv_usec)/(double)1000000;
cout << "Time to create 40K bins: " << tE - tS << "s" << endl;
printStats(h);
+ h.mergeBins(1, 100);
+
cout << h.numBinsTotal() << endl;
h.addBin(0.1, 0.1, 0.2, 0.2);
h.addBin(110, 0, 200, 12.100);
@@ -160,16 +164,17 @@
/// Addition/Subtraction:
Histo2D first(100, 0, 100, 100, 0, 100);
first.fill(1,1,1);
+ first.fill(23,1,1);
Histo2D second(100, 0, 100, 100, 0, 100);
second.fill(1,1,1);
Histo2D added(first+second);
Histo2D subtracted(first-second);
- //Histo2D divided(first / second);
+ Scatter3D divided(first/second);
printStats(added);
printStats(subtracted);
- //printStats(divided);
+
///And now, test cuts:
More information about the yoda-svn
mailing list
|