InterpolationTables/src/GridAxis.cc

0001 #include <cmath>
0002 #include <algorithm>
0003 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0004
0005 #include "Alignment/Geners/interface/binaryIO.hh"
0006 #include "Alignment/Geners/interface/IOException.hh"
0007
0008 #include "JetMETCorrections/InterpolationTables/interface/GridAxis.h"
0009 #include "JetMETCorrections/InterpolationTables/interface/closeWithinTolerance.h"
0010
0011 namespace npstat {
0012   void GridAxis::initialize() {
0013     if (npt_ <= 1U)
0014       throw npstat::NpstatInvalidArgument(
0015           "In npstat::GridAxis::initialize: insufficient number "
0016           "of coordinates (at least 2 are required)");
0017     std::sort(coords_.begin(), coords_.end());
0018     const double* c = &coords_[0];
0019     if (useLogSpace_) {
0020       logs_.reserve(npt_);
0021       for (unsigned i = 0; i < npt_; ++i) {
0022         if (c[i] <= 0.0)
0023           throw npstat::NpstatInvalidArgument(
0024               "In npstat::GridAxis::initialize: in log space"
0025               "all coordinates must be positive");
0026         logs_.push_back(log(c[i]));
0027       }
0028     }
0029
0030     // Can't have duplicate coordinates
0031     for (unsigned i = 1; i < npt_; ++i)
0032       if (c[i] <= c[i - 1U])
0033         throw npstat::NpstatInvalidArgument(
0034             "In npstat::GridAxis::initialize: "
0035             "all coordinates must be distinct");
0036   }
0037
0038   GridAxis::GridAxis(const std::vector<double>& coords, const bool useLogSpace)
0039       : coords_(coords), npt_(coords_.size()), useLogSpace_(useLogSpace) {
0040     initialize();
0041   }
0042
0043   GridAxis::GridAxis(const std::vector<double>& coords, const char* label, const bool useLogSpace)
0044       : coords_(coords), label_(label ? label : ""), npt_(coords_.size()), useLogSpace_(useLogSpace) {
0045     initialize();
0046   }
0047
0048   std::pair<unsigned, double> GridAxis::getInterval(const double x) const {
0049     if (useLogSpace_)
0050       if (x <= 0.0)
0051         throw npstat::NpstatDomainError("In GridAxis::getInterval: argument must be positive");
0052     const double* c = &coords_[0];
0053     if (x <= c[0])
0054       return std::pair<unsigned, double>(0U, 1.0);
0055     else if (x >= c[npt_ - 1U])
0056       return std::pair<unsigned, double>(npt_ - 2U, 0.0);
0057     else {
0058       unsigned ibnd = lower_bound(coords_.begin(), coords_.end(), x) - coords_.begin();
0059       --ibnd;
0060       if (useLogSpace_) {
0061         const double* l = &logs_[0];
0062         const double w = 1.0 - (log(x) - l[ibnd]) / (l[ibnd + 1U] - l[ibnd]);
0063         return std::pair<unsigned, double>(ibnd, w);
0064       } else {
0065         const double w = 1.0 - (x - c[ibnd]) / (c[ibnd + 1U] - c[ibnd]);
0066         return std::pair<unsigned, double>(ibnd, w);
0067       }
0068     }
0069   }
0070
0071   std::pair<unsigned, double> GridAxis::linearInterval(const double x) const {
0072     if (useLogSpace_)
0073       if (x <= 0.0)
0074         throw npstat::NpstatDomainError("In GridAxis::linearInterval: argument must be positive");
0075     const double* c = &coords_[0];
0076     if (x <= c[0]) {
0077       if (useLogSpace_) {
0078         const double* l = &logs_[0];
0079         const double bw = l[1] - l[0];
0080         return std::pair<unsigned, double>(0U, 1.0 - (log(x) - l[0]) / bw);
0081       } else {
0082         const double bw = c[1] - c[0];
0083         return std::pair<unsigned, double>(0U, 1.0 - (x - c[0]) / bw);
0084       }
0085     } else if (x >= c[npt_ - 1U]) {
0086       if (useLogSpace_) {
0087         const double* l = &logs_[0];
0088         const double bw = l[npt_ - 1U] - l[npt_ - 2U];
0089         return std::pair<unsigned, double>(npt_ - 2U, (l[npt_ - 1U] - log(x)) / bw);
0090       } else {
0091         const double bw = c[npt_ - 1U] - c[npt_ - 2U];
0092         return std::pair<unsigned, double>(npt_ - 2U, (c[npt_ - 1U] - x) / bw);
0093       }
0094     } else {
0095       unsigned ibnd = lower_bound(coords_.begin(), coords_.end(), x) - coords_.begin();
0096       --ibnd;
0097       if (useLogSpace_) {
0098         const double* l = &logs_[0];
0099         const double w = 1.0 - (log(x) - l[ibnd]) / (l[ibnd + 1U] - l[ibnd]);
0100         return std::pair<unsigned, double>(ibnd, w);
0101       } else {
0102         const double w = 1.0 - (x - c[ibnd]) / (c[ibnd + 1U] - c[ibnd]);
0103         return std::pair<unsigned, double>(ibnd, w);
0104       }
0105     }
0106   }
0107
0108   bool GridAxis::write(std::ostream& of) const {
0109     // It is unlikely that this object will be written in isolation.
0110     // So, do not bother with too many checks.
0111     gs::write_pod_vector(of, coords_);
0112     gs::write_pod(of, label_);
0113     gs::write_pod(of, useLogSpace_);
0114     return !of.fail();
0115   }
0116
0117   GridAxis* GridAxis::read(const gs::ClassId& id, std::istream& in) {
0118     static const gs::ClassId current(gs::ClassId::makeId<GridAxis>());
0119     current.ensureSameId(id);
0120
0121     std::vector<double> coords;
0122     gs::read_pod_vector(in, &coords);
0123
0124     std::string label;
0125     gs::read_pod(in, &label);
0126
0127     bool useLogSpace;
0128     gs::read_pod(in, &useLogSpace);
0129
0130     if (in.fail())
0131       throw gs::IOReadFailure(
0132           "In npstat::GridAxis::read: "
0133           "input stream failure");
0134     return new GridAxis(coords, label.c_str(), useLogSpace);
0135   }
0136
0137   bool GridAxis::operator==(const GridAxis& r) const {
0138     return useLogSpace_ == r.useLogSpace_ && coords_ == r.coords_ && label_ == r.label_;
0139   }
0140
0141   bool GridAxis::isClose(const GridAxis& r, const double tol) const {
0142     if (!(useLogSpace_ == r.useLogSpace_ && label_ == r.label_))
0143       return false;
0144     const unsigned long n = coords_.size();
0145     if (n != r.coords_.size())
0146       return false;
0147     for (unsigned long i = 0; i < n; ++i)
0148       if (!closeWithinTolerance(coords_[i], r.coords_[i], tol))
0149         return false;
0150     return true;
0151   }
0152 }  // namespace npstat