CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoJets/FFTJetAlgorithms/interface/LinInterpolatedTable1D.h

Go to the documentation of this file.
00001 //=========================================================================
00002 // LinInterpolatedTable1D.h
00003 //
00004 // This class can be used to linearly interpolate data in one dimension.
00005 // It differs from similar facilities in the fftjet package by its handling
00006 // of the extrapolation beyound the grid limits.
00007 //
00008 // I. Volobouev
00009 // April 2011
00010 //=========================================================================
00011 
00012 #ifndef RecoJets_FFTJetAlgorithms_LinInterpolatedTable1D_h
00013 #define RecoJets_FFTJetAlgorithms_LinInterpolatedTable1D_h
00014 
00015 #include <utility>
00016 #include <vector>
00017 #include <memory>
00018 #include <algorithm>
00019 
00020 #include "fftjet/SimpleFunctors.hh"
00021 #include "FWCore/Utilities/interface/Exception.h"
00022 
00023 namespace fftjetcms {
00024     class LinInterpolatedTable1D : public fftjet::Functor1<double,double>
00025     {
00026     public:
00027         // Constructor from a regularly spaced data. Extrapolation
00028         // from the edge to infinity can be either linear or constant.
00029         // "npoints" must be larger than 1.
00030         template <typename RealN>
00031         LinInterpolatedTable1D(const RealN* data, unsigned npoints,
00032                                double x_min, double x_max,
00033                                bool leftExtrapolationLinear,
00034                                bool rightExtrapolationLinear);
00035 
00036         // Constructor from a list of points, not necessarily regularly
00037         // spaced (but must be sorted in the increasing order). The first
00038         // member of the pair is the x coordinate, the second is the
00039         // tabulated function value. The input list will be interpolated
00040         // to "npoints" internal points linearly.
00041         template <typename RealN>
00042         LinInterpolatedTable1D(const std::vector<std::pair<RealN,RealN> >& v,
00043                                unsigned npoints,
00044                                bool leftExtrapolationLinear,
00045                                bool rightExtrapolationLinear);
00046 
00047         // Constructor which builds a function returning the given constant
00048         explicit LinInterpolatedTable1D(double c);
00049 
00050         inline virtual ~LinInterpolatedTable1D() {}
00051 
00052         // Main calculations are performed inside the following operator
00053         virtual double operator()(const double& x) const;
00054 
00055         // Comparisons (useful for testing)
00056         bool operator==(const LinInterpolatedTable1D& r) const;
00057         inline bool operator!=(const LinInterpolatedTable1D& r) const
00058             {return !(*this == r);}
00059 
00060         // Various simple inspectors
00061         inline double xmin() const {return xmin_;}
00062         inline double xmax() const {return xmax_;}
00063         inline unsigned npoints() const {return npoints_;}
00064         inline bool leftExtrapolationLinear() const
00065             {return leftExtrapolationLinear_;}
00066         inline bool rightExtrapolationLinear() const
00067             {return rightExtrapolationLinear_;}
00068         inline const double* data() const {return &data_[0];}
00069 
00070         // The following checks whether the table is monotonous
00071         // (and, therefore, can be inverted). Possible flat regions
00072         // at the edges are not taken into account.
00073         bool isMonotonous() const;
00074 
00075         // Generate the inverse lookup table. Note that it is only
00076         // possible if the original table is monotonous (not taking
00077         // into account possible flat regions at the edges). If the
00078         // inversion is not possible, NULL pointer will be returned.
00079         //
00080         // The new table will have "npoints" points. The parameters
00081         // "leftExtrapolationLinear" and "rightExtrapolationLinear"
00082         // refer to the inverted table (note that left and right will
00083         // exchange places if the original table is decreasing).
00084         //
00085         std::auto_ptr<LinInterpolatedTable1D> inverse(
00086             unsigned npoints, bool leftExtrapolationLinear,
00087             bool rightExtrapolationLinear) const;
00088 
00089     private:
00090         static inline double interpolateSimple(
00091             const double x0, const double x1,
00092             const double y0, const double y1,
00093             const double x)
00094         {
00095             return y0 + (y1 - y0)*((x - x0)/(x1 - x0));
00096         }
00097 
00098         std::vector<double> data_;
00099         double xmin_;
00100         double xmax_;
00101         double binwidth_;
00102         unsigned npoints_;
00103         bool leftExtrapolationLinear_;
00104         bool rightExtrapolationLinear_;
00105         mutable bool monotonous_;
00106         mutable bool monotonicityKnown_;
00107     };
00108 }
00109 
00110 
00111 // Implementation of the templated constructors
00112 namespace fftjetcms {
00113     template <typename RealN>
00114     inline LinInterpolatedTable1D::LinInterpolatedTable1D(
00115         const RealN* data, const unsigned npoints,
00116         const double x_min, const double x_max,
00117         const bool leftExtrapolationLinear,
00118         const bool rightExtrapolationLinear)
00119         : data_(data, data+npoints),
00120           xmin_(x_min),
00121           xmax_(x_max),
00122           binwidth_((x_max - x_min)/(npoints - 1U)),
00123           npoints_(npoints),
00124           leftExtrapolationLinear_(leftExtrapolationLinear),
00125           rightExtrapolationLinear_(rightExtrapolationLinear),
00126           monotonous_(false),
00127           monotonicityKnown_(false)
00128     {
00129         if (!data)
00130             throw cms::Exception("FFTJetBadConfig")
00131                 << "No data configured" << std::endl;
00132         if (npoints <= 1U)
00133             throw cms::Exception("FFTJetBadConfig")
00134                 << "Not enough data points" << std::endl;
00135     }
00136 
00137     template <typename RealN>
00138     LinInterpolatedTable1D::LinInterpolatedTable1D(
00139         const std::vector<std::pair<RealN,RealN> >& v,
00140         const unsigned npoints,
00141         const bool leftExtrapolationLinear,
00142         const bool rightExtrapolationLinear)
00143         : xmin_(v[0].first),
00144           xmax_(v[v.size() - 1U].first),
00145           binwidth_((xmax_ - xmin_)/(npoints - 1U)),
00146           npoints_(npoints),
00147           leftExtrapolationLinear_(leftExtrapolationLinear),
00148           rightExtrapolationLinear_(rightExtrapolationLinear),
00149           monotonous_(false),
00150           monotonicityKnown_(false)
00151     {
00152         const unsigned len = v.size();
00153         if (len <= 1U)
00154             throw cms::Exception("FFTJetBadConfig")
00155                 << "Not enough data for interpolation"
00156                 << std::endl;
00157 
00158         if (npoints <= 1U)
00159             throw cms::Exception("FFTJetBadConfig")
00160                 << "Not enough interpolation table entries"
00161                 << std::endl;
00162 
00163         const std::pair<RealN,RealN>* vdata = &v[0];
00164         for (unsigned i=1; i<len; ++i)
00165             if (vdata[i-1U].first >= vdata[i].first)
00166                 throw cms::Exception("FFTJetBadConfig")
00167                     << "Input data is not sorted properly"
00168                     << std::endl;
00169 
00170         unsigned shift = 0U;
00171         if (leftExtrapolationLinear)
00172         {
00173             ++npoints_;
00174             xmin_ -= binwidth_;
00175             shift = 1U;
00176         }
00177         if (rightExtrapolationLinear)
00178         {
00179             ++npoints_;
00180             xmax_ += binwidth_;
00181         }
00182 
00183         data_.resize(npoints_);
00184 
00185         if (leftExtrapolationLinear)
00186         {
00187             data_[0] = interpolateSimple(
00188                 vdata[0].first, vdata[1].first,
00189                 vdata[0].second, vdata[1].second, xmin_);
00190         }
00191         if (rightExtrapolationLinear)
00192         {
00193             data_[npoints_-1U] = interpolateSimple(
00194                 vdata[len - 2U].first, vdata[len - 1U].first,
00195                 vdata[len - 2U].second, vdata[len - 1U].second, xmax_);
00196         }
00197 
00198         data_[shift] = vdata[0].second;
00199         data_[npoints - 1U + shift] = vdata[len - 1U].second;
00200         unsigned ibelow = 0, iabove = 1;
00201         for (unsigned i=1; i<npoints-1; ++i)
00202         {
00203             const double x = xmin_ + (i + shift)*binwidth_;
00204             while (static_cast<double>(v[iabove].first) <= x)
00205             {
00206                 ++ibelow;
00207                 ++iabove;
00208             }
00209             if (v[ibelow].first == v[iabove].first)
00210                 data_[i + shift] = (v[ibelow].second + v[iabove].second)/2.0;
00211             else
00212                 data_[i + shift] = interpolateSimple(
00213                     v[ibelow].first, v[iabove].first,
00214                     v[ibelow].second, v[iabove].second, x);
00215         }
00216     }
00217 }
00218 
00219 #endif // RecoJets_FFTJetAlgorithms_LinInterpolatedTable1D_h