Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:25:22

0001 //=========================================================================
0002 // LinInterpolatedTable1D.h
0003 //
0004 // This class can be used to linearly interpolate data in one dimension.
0005 // It differs from similar facilities in the fftjet package by its handling
0006 // of the extrapolation beyound the grid limits.
0007 //
0008 // I. Volobouev
0009 // April 2011
0010 //=========================================================================
0011 
0012 #ifndef RecoJets_FFTJetAlgorithms_LinInterpolatedTable1D_h
0013 #define RecoJets_FFTJetAlgorithms_LinInterpolatedTable1D_h
0014 
0015 #include <utility>
0016 #include <vector>
0017 #include <memory>
0018 #include <algorithm>
0019 
0020 #include "fftjet/SimpleFunctors.hh"
0021 #include "FWCore/Utilities/interface/Exception.h"
0022 
0023 namespace fftjetcms {
0024   class LinInterpolatedTable1D : public fftjet::Functor1<double, double> {
0025   public:
0026     // Constructor from a regularly spaced data. Extrapolation
0027     // from the edge to infinity can be either linear or constant.
0028     // "npoints" must be larger than 1.
0029     template <typename RealN>
0030     LinInterpolatedTable1D(const RealN* data,
0031                            unsigned npoints,
0032                            double x_min,
0033                            double x_max,
0034                            bool leftExtrapolationLinear,
0035                            bool rightExtrapolationLinear);
0036 
0037     // Constructor from a list of points, not necessarily regularly
0038     // spaced (but must be sorted in the increasing order). The first
0039     // member of the pair is the x coordinate, the second is the
0040     // tabulated function value. The input list will be interpolated
0041     // to "npoints" internal points linearly.
0042     template <typename RealN>
0043     LinInterpolatedTable1D(const std::vector<std::pair<RealN, RealN> >& v,
0044                            unsigned npoints,
0045                            bool leftExtrapolationLinear,
0046                            bool rightExtrapolationLinear);
0047 
0048     // Constructor which builds a function returning the given constant
0049     explicit LinInterpolatedTable1D(double c);
0050 
0051     inline ~LinInterpolatedTable1D() override {}
0052 
0053     // Main calculations are performed inside the following operator
0054     double operator()(const double& x) const override;
0055 
0056     // Comparisons (useful for testing)
0057     bool operator==(const LinInterpolatedTable1D& r) const;
0058     inline bool operator!=(const LinInterpolatedTable1D& r) const { return !(*this == r); }
0059 
0060     // Various simple inspectors
0061     inline double xmin() const { return xmin_; }
0062     inline double xmax() const { return xmax_; }
0063     inline unsigned npoints() const { return npoints_; }
0064     inline bool leftExtrapolationLinear() const { return leftExtrapolationLinear_; }
0065     inline bool rightExtrapolationLinear() const { return rightExtrapolationLinear_; }
0066     inline const double* data() const { return &data_[0]; }
0067 
0068     // The following checks whether the table is monotonous
0069     // (and, therefore, can be inverted). Possible flat regions
0070     // at the edges are not taken into account.
0071     bool isMonotonous() const;
0072 
0073     // Generate the inverse lookup table. Note that it is only
0074     // possible if the original table is monotonous (not taking
0075     // into account possible flat regions at the edges). If the
0076     // inversion is not possible, NULL pointer will be returned.
0077     //
0078     // The new table will have "npoints" points. The parameters
0079     // "leftExtrapolationLinear" and "rightExtrapolationLinear"
0080     // refer to the inverted table (note that left and right will
0081     // exchange places if the original table is decreasing).
0082     //
0083     std::unique_ptr<LinInterpolatedTable1D> inverse(unsigned npoints,
0084                                                     bool leftExtrapolationLinear,
0085                                                     bool rightExtrapolationLinear) const;
0086 
0087   private:
0088     static inline double interpolateSimple(
0089         const double x0, const double x1, const double y0, const double y1, const double x) {
0090       return y0 + (y1 - y0) * ((x - x0) / (x1 - x0));
0091     }
0092 
0093     std::vector<double> data_;
0094     double xmin_;
0095     double xmax_;
0096     double binwidth_;
0097     unsigned npoints_;
0098     bool leftExtrapolationLinear_;
0099     bool rightExtrapolationLinear_;
0100     mutable bool monotonous_;
0101     mutable bool monotonicityKnown_;
0102   };
0103 }  // namespace fftjetcms
0104 
0105 // Implementation of the templated constructors
0106 namespace fftjetcms {
0107   template <typename RealN>
0108   inline LinInterpolatedTable1D::LinInterpolatedTable1D(const RealN* data,
0109                                                         const unsigned npoints,
0110                                                         const double x_min,
0111                                                         const double x_max,
0112                                                         const bool leftExtrapolationLinear,
0113                                                         const bool rightExtrapolationLinear)
0114       : data_(data, data + npoints),
0115         xmin_(x_min),
0116         xmax_(x_max),
0117         binwidth_((x_max - x_min) / (npoints - 1U)),
0118         npoints_(npoints),
0119         leftExtrapolationLinear_(leftExtrapolationLinear),
0120         rightExtrapolationLinear_(rightExtrapolationLinear),
0121         monotonous_(false),
0122         monotonicityKnown_(false) {
0123     if (!data)
0124       throw cms::Exception("FFTJetBadConfig") << "No data configured" << std::endl;
0125     if (npoints <= 1U)
0126       throw cms::Exception("FFTJetBadConfig") << "Not enough data points" << std::endl;
0127   }
0128 
0129   template <typename RealN>
0130   LinInterpolatedTable1D::LinInterpolatedTable1D(const std::vector<std::pair<RealN, RealN> >& v,
0131                                                  const unsigned npoints,
0132                                                  const bool leftExtrapolationLinear,
0133                                                  const bool rightExtrapolationLinear)
0134       : xmin_(v[0].first),
0135         xmax_(v[v.size() - 1U].first),
0136         binwidth_((xmax_ - xmin_) / (npoints - 1U)),
0137         npoints_(npoints),
0138         leftExtrapolationLinear_(leftExtrapolationLinear),
0139         rightExtrapolationLinear_(rightExtrapolationLinear),
0140         monotonous_(false),
0141         monotonicityKnown_(false) {
0142     const unsigned len = v.size();
0143     if (len <= 1U)
0144       throw cms::Exception("FFTJetBadConfig") << "Not enough data for interpolation" << std::endl;
0145 
0146     if (npoints <= 1U)
0147       throw cms::Exception("FFTJetBadConfig") << "Not enough interpolation table entries" << std::endl;
0148 
0149     const std::pair<RealN, RealN>* vdata = &v[0];
0150     for (unsigned i = 1; i < len; ++i)
0151       if (vdata[i - 1U].first >= vdata[i].first)
0152         throw cms::Exception("FFTJetBadConfig") << "Input data is not sorted properly" << std::endl;
0153 
0154     unsigned shift = 0U;
0155     if (leftExtrapolationLinear) {
0156       ++npoints_;
0157       xmin_ -= binwidth_;
0158       shift = 1U;
0159     }
0160     if (rightExtrapolationLinear) {
0161       ++npoints_;
0162       xmax_ += binwidth_;
0163     }
0164 
0165     data_.resize(npoints_);
0166 
0167     if (leftExtrapolationLinear) {
0168       data_[0] = interpolateSimple(vdata[0].first, vdata[1].first, vdata[0].second, vdata[1].second, xmin_);
0169     }
0170     if (rightExtrapolationLinear) {
0171       data_[npoints_ - 1U] = interpolateSimple(
0172           vdata[len - 2U].first, vdata[len - 1U].first, vdata[len - 2U].second, vdata[len - 1U].second, xmax_);
0173     }
0174 
0175     data_[shift] = vdata[0].second;
0176     data_[npoints - 1U + shift] = vdata[len - 1U].second;
0177     unsigned ibelow = 0, iabove = 1;
0178     for (unsigned i = 1; i < npoints - 1; ++i) {
0179       const double x = xmin_ + (i + shift) * binwidth_;
0180       while (static_cast<double>(v[iabove].first) <= x) {
0181         ++ibelow;
0182         ++iabove;
0183       }
0184       if (v[ibelow].first == v[iabove].first)
0185         data_[i + shift] = (v[ibelow].second + v[iabove].second) / 2.0;
0186       else
0187         data_[i + shift] = interpolateSimple(v[ibelow].first, v[iabove].first, v[ibelow].second, v[iabove].second, x);
0188     }
0189   }
0190 }  // namespace fftjetcms
0191 
0192 #endif  // RecoJets_FFTJetAlgorithms_LinInterpolatedTable1D_h