Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:10:39

0001 #ifndef NPSTAT_BOXND_HH_
0002 #define NPSTAT_BOXND_HH_
0003 
0004 /*!
0005 // \file BoxND.h
0006 //
0007 // \brief Template to represent rectangles, boxes, and hyperboxes
0008 //
0009 // Author: I. Volobouev
0010 //
0011 // March 2010
0012 */
0013 
0014 #include <vector>
0015 
0016 #include "Alignment/Geners/interface/ClassId.hh"
0017 #include "JetMETCorrections/InterpolationTables/interface/Interval.h"
0018 
0019 namespace npstat {
0020   /**
0021     // Class to represent rectangles, boxes, and hyperboxes
0022     */
0023   template <typename Numeric>
0024   struct BoxND : public std::vector<Interval<Numeric> > {
0025     /** Default constructor makes a 0-dimensional box */
0026     inline BoxND() {}
0027 
0028     /** Interval in each dimension is made by its default constructor */
0029     inline explicit BoxND(const unsigned long dim) : std::vector<Interval<Numeric> >(dim) {}
0030 
0031     /** Use the same interval in each dimension */
0032     inline BoxND(const unsigned long dim, const Interval<Numeric>& v) : std::vector<Interval<Numeric> >(dim, v) {}
0033 
0034     /**
0035         // Constructor where one of the limits will be 0 and the other
0036         // will be generated from the given vector (which also determines
0037         // the dimensionality)
0038         */
0039     template <typename Num2>
0040     explicit BoxND(const std::vector<Num2>& limits);
0041 
0042     /** Converting constructor */
0043     template <typename Num2>
0044     explicit BoxND(const BoxND<Num2>& r);
0045 
0046     /**
0047         // Get the data from a box of a different type. This method
0048         // works essentially as a converting assignment operator.
0049         */
0050     template <typename Num2>
0051     BoxND& copyFrom(const BoxND<Num2>& r);
0052 
0053     /** Box dimensionality */
0054     inline unsigned long dim() const { return this->size(); }
0055 
0056     /** Box volume */
0057     Numeric volume() const;
0058 
0059     /**
0060         // Midpoint for every coordinate. The size of the "coord"
0061         // array should be at least as large as the box dimensionality.
0062         */
0063     void getMidpoint(Numeric* coord, unsigned long coordLen) const;
0064 
0065     //@{
0066     /**
0067         // This method return "true" if the corresponding function
0068         // of the Interval returns "true" for every coordinate.
0069         // There must be an automatic conversion from Num2 type into Numeric.
0070         */
0071     template <typename Num2>
0072     bool isInsideLower(const Num2* coord, unsigned long coordLen) const;
0073     template <typename Num2>
0074     bool isInsideUpper(const Num2* coord, unsigned long coordLen) const;
0075     template <typename Num2>
0076     bool isInsideWithBounds(const Num2* coord, unsigned long coordLen) const;
0077     template <typename Num2>
0078     bool isInside(const Num2* coord, unsigned long coordLen) const;
0079     //@}
0080 
0081     //@{
0082     /** Scaling of all limits by a constant */
0083     BoxND& operator*=(double r);
0084     BoxND& operator/=(double r);
0085     //@}
0086 
0087     //@{
0088     /** Scaling by a different constant in each dimension */
0089     BoxND& operator*=(const std::vector<double>& scales);
0090     BoxND& operator/=(const std::vector<double>& scales);
0091     //@}
0092 
0093     /**
0094         // Scaling of all limits by a constant in such a way that the midpoint
0095         // remains unchanged
0096         */
0097     BoxND& expand(double r);
0098 
0099     //@{
0100     /**
0101         // Scaling of all limits in such a way that the midpoint
0102         // remains unchanged, using a different scaling factor
0103         // in each dimension
0104         */
0105     BoxND& expand(const std::vector<double>& scales);
0106     BoxND& expand(const double* scales, unsigned long lenScales);
0107     //@}
0108 
0109     //@{
0110     /** Shifting this object */
0111     template <typename Num2>
0112     BoxND& operator+=(const std::vector<Num2>& shifts);
0113     template <typename Num2>
0114     BoxND& operator-=(const std::vector<Num2>& shifts);
0115     template <typename Num2>
0116     BoxND& shift(const Num2* shifts, unsigned long lenShifts);
0117     //@}
0118 
0119     /** Moving this object so that the midpoint is (0, 0, ..., 0) */
0120     BoxND& moveToOrigin();
0121 
0122     /** Overlap volume with another box */
0123     Numeric overlapVolume(const BoxND& r) const;
0124 
0125     /** A faster way to calculate overlapVolume(r)/volume() */
0126     double overlapFraction(const BoxND& r) const;
0127 
0128     /** Box with lower limit 0 and upper limit 1 in all coordinates */
0129     static BoxND unitBox(unsigned long ndim);
0130 
0131     /**
0132         // Box with lower limit -1 and upper limit 1 in all coordinates.
0133         // Note that this will produce nonsense in case the Numeric type
0134         // is unsigned.
0135         */
0136     static BoxND sizeTwoBox(unsigned long ndim);
0137 
0138     /**
0139         // Box with all upper limits set to maximum possible Numeric
0140         // number and with lower limits set to negative maximum (this
0141         // will not work with unsigned long types)
0142         */
0143     static BoxND allSpace(unsigned long ndim);
0144 
0145     //@{
0146     /** Methods related to I/O */
0147     inline gs::ClassId classId() const { return gs::ClassId(*this); }
0148     bool write(std::ostream& of) const;
0149     //@}
0150 
0151     static const char* classname();
0152     static inline unsigned version() { return 1; }
0153     static void restore(const gs::ClassId& id, std::istream& in, BoxND* box);
0154   };
0155 }  // namespace npstat
0156 
0157 //@{
0158 /** Binary comparison for equality */
0159 template <typename Numeric>
0160 bool operator==(const npstat::BoxND<Numeric>& l, const npstat::BoxND<Numeric>& r);
0161 
0162 template <typename Numeric>
0163 bool operator!=(const npstat::BoxND<Numeric>& l, const npstat::BoxND<Numeric>& r);
0164 //@}
0165 
0166 #include <limits>
0167 #include <cassert>
0168 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0169 
0170 #include "Alignment/Geners/interface/GenericIO.hh"
0171 #include "Alignment/Geners/interface/IOIsUnsigned.hh"
0172 
0173 namespace npstat {
0174   template <typename Numeric>
0175   template <typename Num2>
0176   BoxND<Numeric>::BoxND(const BoxND<Num2>& r) {
0177     const unsigned long dim = r.size();
0178     if (dim) {
0179       this->reserve(dim);
0180       for (unsigned long i = 0; i < dim; ++i) {
0181         const Interval<Num2>& ri(r[i]);
0182         this->push_back(Interval<Numeric>(ri.min(), ri.max()));
0183       }
0184     }
0185   }
0186 
0187   template <typename Numeric>
0188   template <typename Num2>
0189   BoxND<Numeric>::BoxND(const std::vector<Num2>& limits) {
0190     const unsigned long dim = limits.size();
0191     if (dim) {
0192       this->reserve(dim);
0193       Numeric zero = Numeric();
0194       for (unsigned long i = 0; i < dim; ++i) {
0195         const Numeric value(static_cast<Numeric>(limits[i]));
0196         if (value >= zero)
0197           this->push_back(Interval<Numeric>(zero, value));
0198         else
0199           this->push_back(Interval<Numeric>(value, zero));
0200       }
0201     }
0202   }
0203 
0204   template <typename Numeric>
0205   template <typename Num2>
0206   BoxND<Numeric>& BoxND<Numeric>::copyFrom(const BoxND<Num2>& r) {
0207     if ((void*)this == (void*)(&r))
0208       return *this;
0209     const unsigned long n = r.size();
0210     this->clear();
0211     this->reserve(n);
0212     for (unsigned long i = 0; i < n; ++i) {
0213       const Interval<Num2>& ir(r[i]);
0214       this->push_back(Interval<Numeric>(ir.min(), ir.max()));
0215     }
0216     return *this;
0217   }
0218 
0219   template <typename Numeric>
0220   Numeric BoxND<Numeric>::volume() const {
0221     Numeric v(static_cast<Numeric>(1));
0222     const unsigned long mydim = this->size();
0223     for (unsigned long i = 0U; i < mydim; ++i)
0224       v *= (*this)[i].length();
0225     return v;
0226   }
0227 
0228   template <typename Numeric>
0229   Numeric BoxND<Numeric>::overlapVolume(const BoxND& r) const {
0230     const unsigned long mydim = this->size();
0231     if (mydim == r.size()) {
0232       Numeric v(static_cast<Numeric>(1));
0233       for (unsigned long i = 0U; i < mydim; ++i)
0234         v *= (*this)[i].overlapLength(r[i]);
0235       return v;
0236     } else
0237       return static_cast<Numeric>(0);
0238   }
0239 
0240   template <typename Numeric>
0241   double BoxND<Numeric>::overlapFraction(const BoxND& r) const {
0242     const unsigned long mydim = this->size();
0243     if (mydim == r.size()) {
0244       double f = 1.0;
0245       for (unsigned long i = 0U; i < mydim; ++i)
0246         f *= (*this)[i].overlapFraction(r[i]);
0247       return f;
0248     } else
0249       return 0.0;
0250   }
0251 
0252   template <typename Numeric>
0253   void BoxND<Numeric>::getMidpoint(Numeric* coord, const unsigned long coordLen) const {
0254     const unsigned long mydim = this->size();
0255     if (coordLen < mydim)
0256       throw npstat::NpstatInvalidArgument("In npstat::BoxND::getMidpoint: insufficient output buffer length");
0257     if (mydim) {
0258       assert(coord);
0259       for (unsigned long i = 0U; i < mydim; ++i)
0260         coord[i] = (*this)[i].midpoint();
0261     }
0262   }
0263 
0264   template <typename Numeric>
0265   template <typename Num2>
0266   bool BoxND<Numeric>::isInsideLower(const Num2* coords, const unsigned long coordLen) const {
0267     if (coordLen != this->size())
0268       throw npstat::NpstatInvalidArgument(
0269           "In npstat::BoxND::isInsideLower: "
0270           "incompatible point dimensionality");
0271     const Interval<Numeric>* myptr = &(*this)[0];
0272     for (unsigned long i = 0; i < coordLen; ++i)
0273       if (!myptr[i].isInsideLower(coords[i]))
0274         return false;
0275     return true;
0276   }
0277 
0278   template <typename Numeric>
0279   template <typename Num2>
0280   bool BoxND<Numeric>::isInsideUpper(const Num2* coords, const unsigned long coordLen) const {
0281     if (coordLen != this->size())
0282       throw npstat::NpstatInvalidArgument(
0283           "In npstat::BoxND::isInsideUpper: "
0284           "incompatible point dimensionality");
0285     const Interval<Numeric>* myptr = &(*this)[0];
0286     for (unsigned long i = 0; i < coordLen; ++i)
0287       if (!myptr[i].isInsideUpper(coords[i]))
0288         return false;
0289     return true;
0290   }
0291 
0292   template <typename Numeric>
0293   template <typename Num2>
0294   bool BoxND<Numeric>::isInsideWithBounds(const Num2* coords, const unsigned long coordLen) const {
0295     if (coordLen != this->size())
0296       throw npstat::NpstatInvalidArgument(
0297           "In npstat::BoxND::isInsideWithBounds: "
0298           "incompatible point dimensionality");
0299     const Interval<Numeric>* myptr = &(*this)[0];
0300     for (unsigned long i = 0; i < coordLen; ++i)
0301       if (!myptr[i].isInsideWithBounds(coords[i]))
0302         return false;
0303     return true;
0304   }
0305 
0306   template <typename Numeric>
0307   template <typename Num2>
0308   bool BoxND<Numeric>::isInside(const Num2* coords, const unsigned long coordLen) const {
0309     if (coordLen != this->size())
0310       throw npstat::NpstatInvalidArgument("In npstat::BoxND::isInside: incompatible point dimensionality");
0311     const Interval<Numeric>* myptr = &(*this)[0];
0312     for (unsigned long i = 0; i < coordLen; ++i)
0313       if (!myptr[i].isInside(coords[i]))
0314         return false;
0315     return true;
0316   }
0317 
0318   template <typename Numeric>
0319   BoxND<Numeric>& BoxND<Numeric>::operator*=(const double r) {
0320     const unsigned long mydim = this->size();
0321     for (unsigned long i = 0; i < mydim; ++i)
0322       (*this)[i] *= r;
0323     return *this;
0324   }
0325 
0326   template <typename Numeric>
0327   BoxND<Numeric>& BoxND<Numeric>::moveToOrigin() {
0328     const unsigned long mydim = this->size();
0329     for (unsigned long i = 0; i < mydim; ++i)
0330       (*this)[i].moveMidpointTo0();
0331     return *this;
0332   }
0333 
0334   template <typename Numeric>
0335   BoxND<Numeric>& BoxND<Numeric>::expand(const double r) {
0336     const unsigned long mydim = this->size();
0337     for (unsigned long i = 0; i < mydim; ++i)
0338       (*this)[i].expand(r);
0339     return *this;
0340   }
0341 
0342   template <typename Numeric>
0343   BoxND<Numeric>& BoxND<Numeric>::operator*=(const std::vector<double>& scales) {
0344     const unsigned long mydim = this->size();
0345     if (mydim != scales.size())
0346       throw npstat::NpstatInvalidArgument(
0347           "In npstat::BoxND::operator*=: "
0348           "incompatible argument dimensionality");
0349     for (unsigned long i = 0; i < mydim; ++i)
0350       (*this)[i] *= scales[i];
0351     return *this;
0352   }
0353 
0354   template <typename Numeric>
0355   BoxND<Numeric>& BoxND<Numeric>::expand(const std::vector<double>& scales) {
0356     const unsigned long mydim = this->size();
0357     if (mydim != scales.size())
0358       throw npstat::NpstatInvalidArgument("In npstat::BoxND::expand: incompatible argument dimensionality");
0359     for (unsigned long i = 0; i < mydim; ++i)
0360       (*this)[i].expand(scales[i]);
0361     return *this;
0362   }
0363 
0364   template <typename Numeric>
0365   BoxND<Numeric>& BoxND<Numeric>::expand(const double* scales, const unsigned long lenScales) {
0366     const unsigned long mydim = this->size();
0367     if (mydim != lenScales)
0368       throw npstat::NpstatInvalidArgument("In npstat::BoxND::expand: incompatible argument dimensionality");
0369     if (mydim) {
0370       assert(scales);
0371       for (unsigned long i = 0; i < mydim; ++i)
0372         (*this)[i].expand(scales[i]);
0373     }
0374     return *this;
0375   }
0376 
0377   template <typename Numeric>
0378   BoxND<Numeric>& BoxND<Numeric>::operator/=(const double r) {
0379     const unsigned long mydim = this->size();
0380     for (unsigned long i = 0; i < mydim; ++i)
0381       (*this)[i] /= r;
0382     return *this;
0383   }
0384 
0385   template <typename Numeric>
0386   BoxND<Numeric>& BoxND<Numeric>::operator/=(const std::vector<double>& scales) {
0387     const unsigned long mydim = this->size();
0388     if (mydim != scales.size())
0389       throw npstat::NpstatInvalidArgument(
0390           "In npstat::BoxND::operator/=: "
0391           "incompatible argument dimensionality");
0392     for (unsigned long i = 0; i < mydim; ++i)
0393       (*this)[i] /= scales[i];
0394     return *this;
0395   }
0396 
0397   template <typename Numeric>
0398   template <typename Num2>
0399   BoxND<Numeric>& BoxND<Numeric>::operator+=(const std::vector<Num2>& shifts) {
0400     const unsigned long mydim = this->size();
0401     if (mydim != shifts.size())
0402       throw npstat::NpstatInvalidArgument(
0403           "In npstat::BoxND::operator+=: "
0404           "incompatible argument dimensionality");
0405     for (unsigned long i = 0; i < mydim; ++i)
0406       (*this)[i] += static_cast<Numeric>(shifts[i]);
0407     return *this;
0408   }
0409 
0410   template <typename Numeric>
0411   template <typename Num2>
0412   BoxND<Numeric>& BoxND<Numeric>::shift(const Num2* shifts, const unsigned long shiftsLen) {
0413     const unsigned long mydim = this->size();
0414     if (mydim != shiftsLen)
0415       throw npstat::NpstatInvalidArgument("In npstat::BoxND::shift: incompatible argument dimensionality");
0416     if (mydim) {
0417       assert(shifts);
0418       for (unsigned long i = 0; i < mydim; ++i)
0419         (*this)[i] += static_cast<Numeric>(shifts[i]);
0420     }
0421     return *this;
0422   }
0423 
0424   template <typename Numeric>
0425   template <typename Num2>
0426   BoxND<Numeric>& BoxND<Numeric>::operator-=(const std::vector<Num2>& shifts) {
0427     const unsigned long mydim = this->size();
0428     if (mydim != shifts.size())
0429       throw npstat::NpstatInvalidArgument(
0430           "In npstat::BoxND::operator-=: "
0431           "incompatible argument dimensionality");
0432     for (unsigned long i = 0; i < mydim; ++i)
0433       (*this)[i] -= static_cast<Numeric>(shifts[i]);
0434     return *this;
0435   }
0436 
0437   template <typename Numeric>
0438   BoxND<Numeric> BoxND<Numeric>::unitBox(const unsigned long ndim) {
0439     Interval<Numeric> unit(static_cast<Numeric>(0), static_cast<Numeric>(1));
0440     return BoxND<Numeric>(ndim, unit);
0441   }
0442 
0443   template <typename Numeric>
0444   BoxND<Numeric> BoxND<Numeric>::sizeTwoBox(const unsigned long ndim) {
0445     const Numeric one = static_cast<Numeric>(1);
0446     Interval<Numeric> i(-one, one);
0447     return BoxND<Numeric>(ndim, i);
0448   }
0449 
0450   template <typename Numeric>
0451   BoxND<Numeric> BoxND<Numeric>::allSpace(const unsigned long ndim) {
0452     const Numeric maxval = std::numeric_limits<Numeric>::max();
0453     Interval<Numeric> i(static_cast<Numeric>(0), maxval);
0454     if (!gs::IOIsUnsigned<Numeric>::value)
0455       i.setMin(-maxval);
0456     return BoxND<Numeric>(ndim, i);
0457   }
0458 
0459   template <typename Numeric>
0460   const char* BoxND<Numeric>::classname() {
0461     static const std::string na(gs::template_class_name<Numeric>("npstat::BoxND"));
0462     return na.c_str();
0463   }
0464 
0465   template <typename Numeric>
0466   bool BoxND<Numeric>::write(std::ostream& of) const {
0467     const unsigned long mydim = this->size();
0468     std::vector<Numeric> limits;
0469     limits.reserve(2UL * mydim);
0470     for (unsigned long i = 0; i < mydim; ++i) {
0471       limits.push_back((*this)[i].min());
0472       limits.push_back((*this)[i].max());
0473     }
0474     return gs::write_item(of, limits);
0475   }
0476 
0477   template <typename Numeric>
0478   void BoxND<Numeric>::restore(const gs::ClassId& id, std::istream& in, BoxND* b) {
0479     static const gs::ClassId current(gs::ClassId::makeId<BoxND<Numeric> >());
0480     current.ensureSameId(id);
0481 
0482     std::vector<Numeric> limits;
0483     gs::restore_item(in, &limits);
0484     if (in.fail())
0485       throw gs::IOReadFailure("In npstat::BoxND::restore: input stream failure");
0486     const unsigned long nlimits = limits.size();
0487     if (nlimits % 2UL)
0488       throw gs::IOInvalidData("In npstat::BoxND::restore: bad limits");
0489     assert(b);
0490     b->clear();
0491     b->reserve(nlimits / 2UL);
0492     for (unsigned long i = 0; i < nlimits / 2UL; ++i)
0493       b->push_back(npstat::Interval<Numeric>(limits[2U * i], limits[2U * i + 1U]));
0494   }
0495 }  // namespace npstat
0496 
0497 template <typename Numeric>
0498 bool operator==(const npstat::BoxND<Numeric>& l, const npstat::BoxND<Numeric>& r) {
0499   const unsigned long dim = l.size();
0500   if (dim != r.size())
0501     return false;
0502   for (unsigned long i = 0; i < dim; ++i)
0503     if (l[i] != r[i])
0504       return false;
0505   return true;
0506 }
0507 
0508 template <typename Numeric>
0509 bool operator!=(const npstat::BoxND<Numeric>& l, const npstat::BoxND<Numeric>& r) {
0510   return !(l == r);
0511 }
0512 
0513 #endif  // NPSTAT_BOXND_HH_