InterpolationTables/interface/Interval.h

0001 #ifndef NPSTAT_INTERVAL_HH_
0002 #define NPSTAT_INTERVAL_HH_
0003
0004 /*!
0005 // \file Interval.h
0006 //
0007 // \brief Template to represent intervals in one dimension
0008 //
0009 // Author: I. Volobouev
0010 //
0011 // March 2010
0012 */
0013
0014 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0015 #include <algorithm>
0016
0017 namespace npstat {
0018   /**
0019     // Representation of 1-d intervals. The following invariant
0020     // is maintained: min() will not exceed max().
0021     //
0022     // See BoxND class for rectangles, boxes, and hyperboxes.
0023     */
0024   template <typename Numeric>
0025   class Interval {
0026   public:
0027     /** Both lower and upper interval bounds are set to Numeric() */
0028     inline Interval() : min_(Numeric()), max_(Numeric()) {}
0029
0030     /**
0031         // Minimum is set to Numeric(), maximum to the given argument.
0032         // An exception is thrown if the argument is smaller than Numeric().
0033         */
0034     inline explicit Interval(const Numeric max) : min_(Numeric()), max_(max) {
0035       if (min_ > max_)
0036         throw npstat::NpstatInvalidArgument("In npstat::Interval constructor: invalid limits");
0037     }
0038
0039     /**
0040         // Constructor from both bounds. Set "swapIfOutOfOrder" argument
0041         // to "true" if the minimum can be larger than the maximum (in this
0042         // case the bounds will be swapped internally).
0043         */
0044     inline Interval(const Numeric min, const Numeric max, const bool swapIfOutOfOrder = false) : min_(min), max_(max) {
0045       if (min_ > max_) {
0046         if (swapIfOutOfOrder)
0047           std::swap(min_, max_);
0048         else
0049           throw npstat::NpstatInvalidArgument("In npstat::Interval constructor: invalid limits");
0050       }
0051     }
0052
0053     /** Set the lower interval bound */
0054     inline void setMin(const Numeric value) {
0055       if (value > max_)
0056         throw npstat::NpstatInvalidArgument("In npstat::Interval::setMin: argument above max");
0057       min_ = value;
0058     }
0059
0060     /** Set the upper interval bound */
0061     inline void setMax(const Numeric value) {
0062       if (value < min_)
0063         throw npstat::NpstatInvalidArgument("In npstat::Interval::setMax: argument below min");
0064       max_ = value;
0065     }
0066
0067     /** Set both interval bounds */
0068     inline void setBounds(const Numeric minval, const Numeric maxval, const bool swapIfOutOfOrder = false) {
0069       if (maxval < minval && !swapIfOutOfOrder)
0070         throw npstat::NpstatInvalidArgument("In npstat::Interval::setBounds: invalid limits");
0071       min_ = minval;
0072       max_ = maxval;
0073       if (swapIfOutOfOrder && min_ > max_)
0074         std::swap(min_, max_);
0075     }
0076
0077     /** Return the lower bound */
0078     inline const Numeric min() const { return min_; }
0079
0080     /** Return the upper bound */
0081     inline const Numeric max() const { return max_; }
0082
0083     /** Return both bounds */
0084     inline void getBounds(Numeric* pmin, Numeric* pmax) const {
0085       *pmin = min_;
0086       *pmax = max_;
0087     }
0088
0089     /** Interval length */
0090     inline Numeric length() const { return max_ - min_; }
0091
0092     /** The middle point of the interval */
0093     inline Numeric midpoint() const { return static_cast<Numeric>((max_ + min_) * 0.5); }
0094
0095     /** Is the point inside the interval or on the lower boundary? */
0096     inline bool isInsideLower(const Numeric value) const { return value >= min_ && value < max_; }
0097
0098     /** Is the point inside the interval or on the upper boundary? */
0099     inline bool isInsideUpper(const Numeric value) const { return value > min_ && value <= max_; }
0100
0101     /** Is the point inside the interval or on one of the boundaries? */
0102     inline bool isInsideWithBounds(const Numeric value) const { return value >= min_ && value <= max_; }
0103
0104     /**
0105         // Is the point completely inside the interval
0106         // (and does not coincide with any bound)?
0107         */
0108     inline bool isInside(const Numeric value) const { return value > min_ && value < max_; }
0109
0110     //@{
0111     /**
0112         // Scaling of both the minimum and the maximum by a constant.
0113         // Minimum and maximum will be swapped internally in case the
0114         // constant is negative.
0115         */
0116     Interval& operator*=(double r);
0117     Interval& operator/=(double r);
0118     //@}
0119
0120     //@{
0121     /** Shift both interval bounds by a constant */
0122     Interval& operator+=(const Numeric value);
0123     Interval& operator-=(const Numeric value);
0124     //@}
0125
0126     /** Move the interval so that the midpoint ends up at 0 */
0127     Interval& moveMidpointTo0();
0128
0129     //@{
0130     /**
0131         // Scaling the bounds by a constant in such a way
0132         // that the midpoint remains unchanged
0133         */
0134     Interval& expand(double r);
0135     //@}
0136
0137     /**
0138         // The following function returns default-constructed empty interval
0139         // in case this interval and the argument interval do not overlap
0140         */
0141     Interval overlap(const Interval& r) const;
0142
0143     /** Same as overlap.length() but a tad faster */
0144     Numeric overlapLength(const Interval& r) const;
0145
0146     /** Same as overlapLength(r)/length() but a tad faster */
0147     double overlapFraction(const Interval& r) const;
0148
0149     /**
0150         // Derive the coefficients a and b such that the linear
0151         // mapping y = a*x + b maps the lower limit of this interval
0152         // into the lower limit of the argument interval and the
0153         // upper limit of this interval into the upper limit of the
0154         // argument interval
0155         */
0156     template <typename Num2>
0157     void linearMap(const Interval<Num2>& r, double* a, double* b) const;
0158
0159   private:
0160     Numeric min_;
0161     Numeric max_;
0162   };
0163 }  // namespace npstat
0164
0165 //@{
0166 /** Binary comparison for equality */
0167 template <typename Numeric>
0168 bool operator==(const npstat::Interval<Numeric>& l, const npstat::Interval<Numeric>& r);
0169
0170 template <typename Numeric>
0171 bool operator!=(const npstat::Interval<Numeric>& l, const npstat::Interval<Numeric>& r);
0172 //@}
0173
0174 #include <cmath>
0175 #include <cassert>
0176
0177 namespace npstat {
0178   template <typename Numeric>
0179   inline Interval<Numeric> Interval<Numeric>::overlap(const Interval<Numeric>& r) const {
0180     Interval<Numeric> over;
0181     if (max_ == r.min_)
0182       over.setBounds(max_, max_);
0183     else if (r.max_ == min_)
0184       over.setBounds(min_, min_);
0185     else if (max_ > r.min_ && r.max_ > min_) {
0186       over.min_ = min_ < r.min_ ? r.min_ : min_;
0187       over.max_ = max_ < r.max_ ? max_ : r.max_;
0188     }
0189     return over;
0190   }
0191
0192   template <typename Numeric>
0193   inline Numeric Interval<Numeric>::overlapLength(const Interval& r) const {
0194     if (max_ > r.min_ && r.max_ > min_) {
0195       const Numeric mn = min_ < r.min_ ? r.min_ : min_;
0196       const Numeric mx = max_ < r.max_ ? max_ : r.max_;
0197       return mx - mn;
0198     } else
0199       return Numeric();
0200   }
0201
0202   template <typename Numeric>
0203   inline double Interval<Numeric>::overlapFraction(const Interval& r) const {
0204     if (max_ > r.min_ && r.max_ > min_) {
0205       const Numeric mn = min_ < r.min_ ? r.min_ : min_;
0206       const Numeric mx = max_ < r.max_ ? max_ : r.max_;
0207       return (mx - mn) * 1.0 / (max_ - min_);
0208     } else
0209       return 0.0;
0210   }
0211
0212   template <typename Numeric>
0213   inline Interval<Numeric>& Interval<Numeric>::operator*=(const double r) {
0214     min_ *= r;
0215     max_ *= r;
0216     if (max_ < min_)
0217       std::swap(min_, max_);
0218     return *this;
0219   }
0220
0221   template <typename Numeric>
0222   inline Interval<Numeric>& Interval<Numeric>::moveMidpointTo0() {
0223     const Numeric len = max_ - min_;
0224     max_ = len / static_cast<Numeric>(2);
0225     min_ = -max_;
0226     return *this;
0227   }
0228
0229   template <typename Numeric>
0230   inline Interval<Numeric>& Interval<Numeric>::expand(const double ir) {
0231     const double r = fabs(ir);
0232     if (r != 1.0) {
0233       const Numeric center(static_cast<Numeric>((max_ + min_) * 0.5));
0234       min_ = center + (min_ - center) * r;
0235       max_ = center + (max_ - center) * r;
0236     }
0237     return *this;
0238   }
0239
0240   template <typename Numeric>
0241   inline Interval<Numeric>& Interval<Numeric>::operator/=(const double r) {
0242     if (!r)
0243       throw npstat::NpstatDomainError("In npstat::Interval::operator/=: division by zero");
0244     min_ /= r;
0245     max_ /= r;
0246     if (max_ < min_)
0247       std::swap(min_, max_);
0248     return *this;
0249   }
0250
0251   template <typename Numeric>
0252   inline Interval<Numeric>& Interval<Numeric>::operator+=(const Numeric r) {
0253     min_ += r;
0254     max_ += r;
0255     return *this;
0256   }
0257
0258   template <typename Numeric>
0259   inline Interval<Numeric>& Interval<Numeric>::operator-=(const Numeric r) {
0260     min_ -= r;
0261     max_ -= r;
0262     return *this;
0263   }
0264
0265   template <typename Numeric>
0266   template <typename Num2>
0267   void Interval<Numeric>::linearMap(const Interval<Num2>& r, double* a, double* b) const {
0268     if (max_ == min_)
0269       throw npstat::NpstatDomainError("In npstat::Interval::linearMap: zero length interval");
0270     assert(a);
0271     assert(b);
0272     const Num2 rmax(r.max());
0273     const Num2 rmin(r.min());
0274     *a = static_cast<double>((rmax - rmin) * 1.0 / (max_ - min_));
0275     *b = static_cast<double>((rmax + rmin) - *a * (max_ + min_)) / 2.0;
0276   }
0277 }  // namespace npstat
0278
0279 template <typename Numeric>
0280 bool operator==(const npstat::Interval<Numeric>& l, const npstat::Interval<Numeric>& r) {
0281   return r.min() == l.min() && r.max() == l.max();
0282 }
0283
0284 template <typename Numeric>
0285 bool operator!=(const npstat::Interval<Numeric>& l, const npstat::Interval<Numeric>& r) {
0286   return !(l == r);
0287 }
0288
0289 #endif  // NPSTAT_INTERVAL_HH_