File indexing completed on 2024-04-06 12:25:22
0001
0002
0003
0004
0005
0006
0007
0008
0009
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
0027
0028
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
0038
0039
0040
0041
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
0049 explicit LinInterpolatedTable1D(double c);
0050
0051 inline ~LinInterpolatedTable1D() override {}
0052
0053
0054 double operator()(const double& x) const override;
0055
0056
0057 bool operator==(const LinInterpolatedTable1D& r) const;
0058 inline bool operator!=(const LinInterpolatedTable1D& r) const { return !(*this == r); }
0059
0060
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
0069
0070
0071 bool isMonotonous() const;
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
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 }
0104
0105
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 }
0191
0192 #endif