Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #ifndef NPSTAT_LININTERPOLATEDTABLEND_HH_
0002 #define NPSTAT_LININTERPOLATEDTABLEND_HH_
0003 
0004 /**
0005 // \file LinInterpolatedTableND.h
0006 //
0007 // \brief Multilinear interpolation/extrapolation on rectangular grids
0008 //
0009 // Author: I. Volobouev
0010 //
0011 // June 2012
0012 */
0013 
0014 #include <climits>
0015 #include <vector>
0016 #include <utility>
0017 
0018 #include <memory>
0019 
0020 #include "JetMETCorrections/InterpolationTables/interface/ArrayND.h"
0021 #include "JetMETCorrections/InterpolationTables/interface/UniformAxis.h"
0022 
0023 namespace npstat {
0024   /** 
0025     // Template for multilinear interpolation/extrapolation of values provided
0026     // on a rectangular coordinate grid. "Numeric" is the type stored in
0027     // the table. "Axis" should be one of GridAxis, UniformAxis, or DualAxis
0028     // classes or a user-provided class with a similar set of methods.
0029     */
0030   template <class Numeric, class Axis = UniformAxis>
0031   class LinInterpolatedTableND {
0032     template <typename Num2, typename Axis2>
0033     friend class LinInterpolatedTableND;
0034 
0035   public:
0036     typedef Numeric value_type;
0037     typedef Axis axis_type;
0038 
0039     /** 
0040         // Main constructor for arbitrary-dimensional interpolators.
0041         //
0042         // "axes" are the axes of the rectangular coordinate grid.
0043         //
0044         // In each pair provided by the "extrapolationType" argument,
0045         // the first element of the pair specifies whether the extrapolation
0046         // to negative infinity should be linear (if "true") or constant
0047         // (if "false"). The second element of the pair specifies whether
0048         // to extrapolate linearly to positive infinity.
0049         //
0050         // "functionLabel" is an arbitrary string which can potentially
0051         // be used by plotting programs and such.
0052         */
0053     LinInterpolatedTableND(const std::vector<Axis>& axes,
0054                            const std::vector<std::pair<bool, bool> >& extrapolationType,
0055                            const char* functionLabel = nullptr);
0056 
0057     /** Convenience constructor for 1-d interpolators */
0058     LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX, const char* functionLabel = nullptr);
0059 
0060     /** Convenience constructor for 2-d interpolators */
0061     LinInterpolatedTableND(const Axis& xAxis,
0062                            bool leftX,
0063                            bool rightX,
0064                            const Axis& yAxis,
0065                            bool leftY,
0066                            bool rightY,
0067                            const char* functionLabel = nullptr);
0068 
0069     /** Convenience constructor for 3-d interpolators */
0070     LinInterpolatedTableND(const Axis& xAxis,
0071                            bool leftX,
0072                            bool rightX,
0073                            const Axis& yAxis,
0074                            bool leftY,
0075                            bool rightY,
0076                            const Axis& zAxis,
0077                            bool leftZ,
0078                            bool rightZ,
0079                            const char* functionLabel = nullptr);
0080 
0081     /** Convenience constructor for 4-d interpolators */
0082     LinInterpolatedTableND(const Axis& xAxis,
0083                            bool leftX,
0084                            bool rightX,
0085                            const Axis& yAxis,
0086                            bool leftY,
0087                            bool rightY,
0088                            const Axis& zAxis,
0089                            bool leftZ,
0090                            bool rightZ,
0091                            const Axis& tAxis,
0092                            bool leftT,
0093                            bool rightT,
0094                            const char* functionLabel = nullptr);
0095 
0096     /** Convenience constructor for 5-d interpolators */
0097     LinInterpolatedTableND(const Axis& xAxis,
0098                            bool leftX,
0099                            bool rightX,
0100                            const Axis& yAxis,
0101                            bool leftY,
0102                            bool rightY,
0103                            const Axis& zAxis,
0104                            bool leftZ,
0105                            bool rightZ,
0106                            const Axis& tAxis,
0107                            bool leftT,
0108                            bool rightT,
0109                            const Axis& vAxis,
0110                            bool leftV,
0111                            bool rightV,
0112                            const char* functionLabel = nullptr);
0113 
0114     /**
0115         // Converting copy constructor from interpolator
0116         // with another storage type
0117         */
0118     template <class Num2>
0119     LinInterpolatedTableND(const LinInterpolatedTableND<Num2, Axis>&);
0120 
0121     /** Default constructor not implemented  **/
0122     LinInterpolatedTableND() = delete;
0123 
0124     /**
0125         // Basic interpolation result. Argument point dimensionality must be
0126         // compatible with the interpolator dimensionality.
0127         */
0128     Numeric operator()(const double* point, unsigned dim) const;
0129 
0130     //@{
0131     /** Convenience function for low-dimensional interpolators */
0132     Numeric operator()(const double& x0) const;
0133     Numeric operator()(const double& x0, const double& x1) const;
0134     Numeric operator()(const double& x0, const double& x1, const double& x2) const;
0135     Numeric operator()(const double& x0, const double& x1, const double& x2, const double& x3) const;
0136     Numeric operator()(const double& x0, const double& x1, const double& x2, const double& x3, const double& x4) const;
0137     //@}
0138 
0139     //@{
0140     /** Examine interpolator contents */
0141     inline unsigned dim() const { return dim_; }
0142     inline const std::vector<Axis>& axes() const { return axes_; }
0143     inline const Axis& axis(const unsigned i) const { return axes_.at(i); }
0144     inline unsigned long length() const { return data_.length(); }
0145     bool leftInterpolationLinear(unsigned i) const;
0146     bool rightInterpolationLinear(unsigned i) const;
0147     std::vector<std::pair<bool, bool> > interpolationType() const;
0148     inline const std::string& functionLabel() const { return functionLabel_; }
0149     //@}
0150 
0151     //@{
0152     /** Access the interpolator table data */
0153     inline const ArrayND<Numeric>& table() const { return data_; }
0154     inline ArrayND<Numeric>& table() { return data_; }
0155     //@}
0156 
0157     /** Convenience function for getting coordinates of the grid points */
0158     void getCoords(unsigned long linearIndex, double* coords, unsigned coordsBufferSize) const;
0159 
0160     /** 
0161         // This method returns "true" if the method isUniform()
0162         // of each interpolator axis returns "true" 
0163         */
0164     bool isUniformlyBinned() const;
0165 
0166     /**
0167         // This method will return "true" if the point
0168         // is inside the grid limits or on the boundary
0169         */
0170     bool isWithinLimits(const double* point, unsigned dim) const;
0171 
0172     /** Modifier for the function label */
0173     inline void setFunctionLabel(const char* newlabel) { functionLabel_ = newlabel ? newlabel : ""; }
0174 
0175     /**
0176         // Invert the function w.r.t. the variable represented by one of
0177         // the axes. The function values must change monotonously along
0178         // the chosen axis. Note that this operation is meaningful only
0179         // in case "Numeric" type is either float or double.
0180         */
0181     template <typename ConvertibleToUnsigned>
0182     std::unique_ptr<LinInterpolatedTableND> invertWRTAxis(ConvertibleToUnsigned axisNumber,
0183                                                           const Axis& replacementAxis,
0184                                                           bool newAxisLeftLinear,
0185                                                           bool newAxisRightLinear,
0186                                                           const char* functionLabel = nullptr) const;
0187 
0188     /**
0189         // This method inverts the ratio response.
0190         // That is, we assume that the table encodes r(x) for
0191         // some function f(x) = x*r(x). We also assume that the
0192         // original axis does not represent x directly -- instead,
0193         // axis coordinates are given by g(x) (in practice, g is
0194         // often the natural log). We will also assume that the new
0195         // axis is not going to represent f(x) directly -- it
0196         // will be h(f(x)) instead. The functors "invg" and "invh"
0197         // below must do the inverse: taking the axes coordinates
0198         // to the actual values of x and f(x). Both "invg" and "invh"
0199         // must be monotonously increasing functions. The code assumes
0200         // that x*r(x) -> 0 when x->0 (that is, r(x) is bounded at 0)
0201         // and it also assumes (but does not check explicitly)
0202         // that x*r(x) is monotonously increasing with x.
0203         //
0204         // The returned interpolation table encodes the values
0205         // of x/f(x). Of course, they are just 1/r(x), but the trick
0206         // is to be able to look them up quickly as a function of
0207         // h(f(x)). Naturally, this whole operation is meaningful
0208         // only in case "Numeric" type is either float or double.
0209         */
0210     template <class Functor1, class Functor2>
0211     std::unique_ptr<LinInterpolatedTableND> invertRatioResponse(unsigned axisNumber,
0212                                                                 const Axis& replacementAxis,
0213                                                                 bool newAxisLeftLinear,
0214                                                                 bool newAxisRightLinear,
0215                                                                 Functor1 invg,
0216                                                                 Functor2 invh,
0217                                                                 const char* functionLabel = nullptr) const;
0218 
0219     /** Comparison for equality */
0220     bool operator==(const LinInterpolatedTableND&) const;
0221 
0222     /** Logical negation of operator== */
0223     inline bool operator!=(const LinInterpolatedTableND& r) const { return !(*this == r); }
0224 
0225     //@{
0226     // Method related to "geners" I/O
0227     inline gs::ClassId classId() const { return gs::ClassId(*this); }
0228     bool write(std::ostream& of) const;
0229     //@}
0230 
0231     static const char* classname();
0232     static inline unsigned version() { return 1; }
0233     static LinInterpolatedTableND* read(const gs::ClassId& id, std::istream& in);
0234 
0235   private:
0236     LinInterpolatedTableND(const ArrayND<Numeric>& data,
0237                            const std::vector<Axis>& axes,
0238                            const char* leftInterpolation,
0239                            const char* rightInterpolation,
0240                            const std::string& label);
0241 
0242     bool allConstInterpolated() const;
0243 
0244     ArrayND<Numeric> data_;
0245     std::vector<Axis> axes_;
0246     std::string functionLabel_;
0247     char leftInterpolationLinear_[CHAR_BIT * sizeof(unsigned long)];
0248     char rightInterpolationLinear_[CHAR_BIT * sizeof(unsigned long)];
0249     unsigned dim_;
0250     bool allConstInterpolated_;
0251 
0252     template <class Functor1>
0253     static double solveForRatioArg(double xmin, double xmax, double rmin, double rmax, double fval, Functor1 invg);
0254 
0255     template <class Functor1>
0256     static void invert1DResponse(const ArrayND<Numeric>& fromSlice,
0257                                  const Axis& fromAxis,
0258                                  const Axis& toAxis,
0259                                  bool newLeftLinear,
0260                                  bool newRightLinear,
0261                                  Functor1 invg,
0262                                  const double* rawx,
0263                                  const double* rawf,
0264                                  double* workspace,
0265                                  ArrayND<Numeric>* toSlice);
0266   };
0267 }  // namespace npstat
0268 
0269 #include <cmath>
0270 #include <cfloat>
0271 #include <cassert>
0272 #include <algorithm>
0273 #include <functional>
0274 
0275 #include "Alignment/Geners/interface/binaryIO.hh"
0276 
0277 #include "JetMETCorrections/InterpolationTables/interface/ArrayNDScanner.h"
0278 #include "JetMETCorrections/InterpolationTables/interface/isMonotonous.h"
0279 
0280 namespace npstat {
0281   namespace Private {
0282     template <class Axis>
0283     ArrayShape makeTableShape(const std::vector<Axis>& axes) {
0284       const unsigned n = axes.size();
0285       ArrayShape result;
0286       result.reserve(n);
0287       for (unsigned i = 0; i < n; ++i)
0288         result.push_back(axes[i].nCoords());
0289       return result;
0290     }
0291 
0292     template <class Axis>
0293     ArrayShape makeTableShape(const Axis& xAxis) {
0294       ArrayShape result;
0295       result.reserve(1U);
0296       result.push_back(xAxis.nCoords());
0297       return result;
0298     }
0299 
0300     template <class Axis>
0301     ArrayShape makeTableShape(const Axis& xAxis, const Axis& yAxis) {
0302       ArrayShape result;
0303       result.reserve(2U);
0304       result.push_back(xAxis.nCoords());
0305       result.push_back(yAxis.nCoords());
0306       return result;
0307     }
0308 
0309     template <class Axis>
0310     ArrayShape makeTableShape(const Axis& xAxis, const Axis& yAxis, const Axis& zAxis) {
0311       ArrayShape result;
0312       result.reserve(3U);
0313       result.push_back(xAxis.nCoords());
0314       result.push_back(yAxis.nCoords());
0315       result.push_back(zAxis.nCoords());
0316       return result;
0317     }
0318 
0319     template <class Axis>
0320     ArrayShape makeTableShape(const Axis& xAxis, const Axis& yAxis, const Axis& zAxis, const Axis& tAxis) {
0321       ArrayShape result;
0322       result.reserve(4U);
0323       result.push_back(xAxis.nCoords());
0324       result.push_back(yAxis.nCoords());
0325       result.push_back(zAxis.nCoords());
0326       result.push_back(tAxis.nCoords());
0327       return result;
0328     }
0329 
0330     template <class Axis>
0331     ArrayShape makeTableShape(
0332         const Axis& xAxis, const Axis& yAxis, const Axis& zAxis, const Axis& tAxis, const Axis& vAxis) {
0333       ArrayShape result;
0334       result.reserve(5U);
0335       result.push_back(xAxis.nCoords());
0336       result.push_back(yAxis.nCoords());
0337       result.push_back(zAxis.nCoords());
0338       result.push_back(tAxis.nCoords());
0339       result.push_back(vAxis.nCoords());
0340       return result;
0341     }
0342 
0343     inline double lind_interpolateSimple(
0344         const double x0, const double x1, const double y0, const double y1, const double x) {
0345       return y0 + (y1 - y0) * ((x - x0) / (x1 - x0));
0346     }
0347 
0348     template <typename Numeric, class Axis>
0349     void lind_invert1DSlice(const ArrayND<Numeric>& fromSlice,
0350                             const Axis& fromAxis,
0351                             const Axis& toAxis,
0352                             const bool leftLinear,
0353                             const bool rightLinear,
0354                             ArrayND<Numeric>* toSlice) {
0355       assert(toSlice);
0356       assert(fromSlice.rank() == 1U);
0357       assert(toSlice->rank() == 1U);
0358 
0359       const Numeric* fromData = fromSlice.data();
0360       const unsigned fromLen = fromSlice.length();
0361       assert(fromLen > 1U);
0362       assert(fromLen == fromAxis.nCoords());
0363       const Numeric* fromDataEnd = fromData + fromLen;
0364       if (!isStrictlyMonotonous(fromData, fromDataEnd))
0365         throw npstat::NpstatInvalidArgument(
0366             "In npstat::Private::lind_invert1DSlice: "
0367             "slice data is not monotonous and can not be inverted");
0368 
0369       const Numeric yfirst = fromData[0];
0370       const Numeric ylast = fromData[fromLen - 1U];
0371       const bool increasing = yfirst < ylast;
0372 
0373       Numeric* toD = const_cast<Numeric*>(toSlice->data());
0374       const unsigned nAxisPoints = toAxis.nCoords();
0375       assert(toSlice->length() == nAxisPoints);
0376 
0377       for (unsigned ipt = 0; ipt < nAxisPoints; ++ipt) {
0378         const Numeric y = static_cast<Numeric>(toAxis.coordinate(ipt));
0379         if (increasing) {
0380           if (y <= yfirst) {
0381             if (leftLinear)
0382               toD[ipt] = Private::lind_interpolateSimple(
0383                   yfirst, fromData[1], fromAxis.coordinate(0), fromAxis.coordinate(1), y);
0384             else
0385               toD[ipt] = fromAxis.coordinate(0);
0386           } else if (y >= ylast) {
0387             if (rightLinear)
0388               toD[ipt] = Private::lind_interpolateSimple(ylast,
0389                                                          fromData[fromLen - 2U],
0390                                                          fromAxis.coordinate(fromLen - 1U),
0391                                                          fromAxis.coordinate(fromLen - 2U),
0392                                                          y);
0393             else
0394               toD[ipt] = fromAxis.coordinate(fromLen - 1U);
0395           } else {
0396             const unsigned i = std::lower_bound(fromData, fromDataEnd, y) - fromData;
0397             toD[ipt] = Private::lind_interpolateSimple(
0398                 fromData[i - 1U], fromData[i], fromAxis.coordinate(i - 1U), fromAxis.coordinate(i), y);
0399           }
0400         } else {
0401           // The role of left and right are exchanged
0402           // with respect to first and last point
0403           if (y <= ylast) {
0404             if (leftLinear)
0405               toD[ipt] = Private::lind_interpolateSimple(ylast,
0406                                                          fromData[fromLen - 2U],
0407                                                          fromAxis.coordinate(fromLen - 1U),
0408                                                          fromAxis.coordinate(fromLen - 2U),
0409                                                          y);
0410             else
0411               toD[ipt] = fromAxis.coordinate(fromLen - 1U);
0412           } else if (y >= yfirst) {
0413             if (rightLinear)
0414               toD[ipt] = Private::lind_interpolateSimple(
0415                   yfirst, fromData[1], fromAxis.coordinate(0), fromAxis.coordinate(1), y);
0416             else
0417               toD[ipt] = fromAxis.coordinate(0);
0418           } else {
0419             const unsigned i = std::lower_bound(fromData, fromDataEnd, y, std::greater<Numeric>()) - fromData;
0420             toD[ipt] = Private::lind_interpolateSimple(
0421                 fromData[i - 1U], fromData[i], fromAxis.coordinate(i - 1U), fromAxis.coordinate(i), y);
0422           }
0423         }
0424       }
0425     }
0426   }  // namespace Private
0427 
0428   template <class Numeric, class Axis>
0429   bool LinInterpolatedTableND<Numeric, Axis>::allConstInterpolated() const {
0430     for (unsigned i = 0; i < dim_; ++i)
0431       if (leftInterpolationLinear_[i] || rightInterpolationLinear_[i])
0432         return false;
0433     return true;
0434   }
0435 
0436   template <class Numeric, class Axis>
0437   bool LinInterpolatedTableND<Numeric, Axis>::operator==(const LinInterpolatedTableND& r) const {
0438     if (dim_ != r.dim_)
0439       return false;
0440     for (unsigned i = 0; i < dim_; ++i) {
0441       if (leftInterpolationLinear_[i] != r.leftInterpolationLinear_[i])
0442         return false;
0443       if (rightInterpolationLinear_[i] != r.rightInterpolationLinear_[i])
0444         return false;
0445     }
0446     return data_ == r.data_ && axes_ == r.axes_ && functionLabel_ == r.functionLabel_;
0447   }
0448 
0449   template <typename Numeric, class Axis>
0450   const char* LinInterpolatedTableND<Numeric, Axis>::classname() {
0451     static const std::string myClass(gs::template_class_name<Numeric, Axis>("npstat::LinInterpolatedTableND"));
0452     return myClass.c_str();
0453   }
0454 
0455   template <typename Numeric, class Axis>
0456   bool LinInterpolatedTableND<Numeric, Axis>::write(std::ostream& of) const {
0457     const bool status = data_.classId().write(of) && data_.write(of) && gs::write_obj_vector(of, axes_);
0458     if (status) {
0459       gs::write_pod_array(of, leftInterpolationLinear_, dim_);
0460       gs::write_pod_array(of, rightInterpolationLinear_, dim_);
0461       gs::write_pod(of, functionLabel_);
0462     }
0463     return status && !of.fail();
0464   }
0465 
0466   template <typename Numeric, class Axis>
0467   LinInterpolatedTableND<Numeric, Axis>* LinInterpolatedTableND<Numeric, Axis>::read(const gs::ClassId& id,
0468                                                                                      std::istream& in) {
0469     static const gs::ClassId current(gs::ClassId::makeId<LinInterpolatedTableND<Numeric, Axis> >());
0470     current.ensureSameId(id);
0471 
0472     gs::ClassId ida(in, 1);
0473     ArrayND<Numeric> data;
0474     ArrayND<Numeric>::restore(ida, in, &data);
0475     std::vector<Axis> axes;
0476     gs::read_heap_obj_vector_as_placed(in, &axes);
0477     const unsigned dim = axes.size();
0478     if (dim > CHAR_BIT * sizeof(unsigned long) || data.rank() != dim)
0479       throw gs::IOInvalidData(
0480           "In npstat::LinInterpolatedTableND::read: "
0481           "read back invalid dimensionality");
0482     char leftInterpolation[CHAR_BIT * sizeof(unsigned long)];
0483     gs::read_pod_array(in, leftInterpolation, dim);
0484     char rightInterpolation[CHAR_BIT * sizeof(unsigned long)];
0485     gs::read_pod_array(in, rightInterpolation, dim);
0486     std::string label;
0487     gs::read_pod(in, &label);
0488     if (in.fail())
0489       throw gs::IOReadFailure("In npstat::LinInterpolatedTableND::read: input stream failure");
0490     return new LinInterpolatedTableND(data, axes, leftInterpolation, rightInterpolation, label);
0491   }
0492 
0493   template <typename Numeric, class Axis>
0494   bool LinInterpolatedTableND<Numeric, Axis>::leftInterpolationLinear(const unsigned i) const {
0495     if (i >= dim_)
0496       throw npstat::NpstatOutOfRange(
0497           "In npstat::LinInterpolatedTableND::leftInterpolationLinear: "
0498           "index out of range");
0499     return leftInterpolationLinear_[i];
0500   }
0501 
0502   template <typename Numeric, class Axis>
0503   bool LinInterpolatedTableND<Numeric, Axis>::rightInterpolationLinear(const unsigned i) const {
0504     if (i >= dim_)
0505       throw npstat::NpstatOutOfRange(
0506           "In npstat::LinInterpolatedTableND::rightInterpolationLinear: "
0507           "index out of range");
0508     return rightInterpolationLinear_[i];
0509   }
0510 
0511   template <typename Numeric, class Axis>
0512   bool LinInterpolatedTableND<Numeric, Axis>::isUniformlyBinned() const {
0513     for (unsigned i = 0; i < dim_; ++i)
0514       if (!axes_[i].isUniform())
0515         return false;
0516     return true;
0517   }
0518 
0519   template <typename Numeric, class Axis>
0520   std::vector<std::pair<bool, bool> > LinInterpolatedTableND<Numeric, Axis>::interpolationType() const {
0521     std::vector<std::pair<bool, bool> > vec;
0522     vec.reserve(dim_);
0523     for (unsigned i = 0; i < dim_; ++i)
0524       vec.push_back(std::pair<bool, bool>(leftInterpolationLinear_[i], rightInterpolationLinear_[i]));
0525     return vec;
0526   }
0527 
0528   template <typename Numeric, class Axis>
0529   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(
0530       const std::vector<Axis>& axes, const std::vector<std::pair<bool, bool> >& interpolationType, const char* label)
0531       : data_(Private::makeTableShape(axes)), axes_(axes), functionLabel_(label ? label : ""), dim_(axes.size()) {
0532     if (dim_ == 0 || dim_ >= CHAR_BIT * sizeof(unsigned long))
0533       throw npstat::NpstatInvalidArgument(
0534           "In npstat::LinInterpolatedTableND constructor: requested "
0535           "table dimensionality is not supported");
0536     if (dim_ != interpolationType.size())
0537       throw npstat::NpstatInvalidArgument(
0538           "In npstat::LinInterpolatedTableND constructor: "
0539           "incompatible number of interpolation specifications");
0540     for (unsigned i = 0; i < dim_; ++i) {
0541       const std::pair<bool, bool>& pair(interpolationType[i]);
0542       leftInterpolationLinear_[i] = pair.first;
0543       rightInterpolationLinear_[i] = pair.second;
0544     }
0545 
0546     allConstInterpolated_ = allConstInterpolated();
0547   }
0548 
0549   template <typename Numeric, class Axis>
0550   template <class Num2>
0551   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(const LinInterpolatedTableND<Num2, Axis>& r)
0552       : data_(r.data_),
0553         axes_(r.axes_),
0554         functionLabel_(r.functionLabel_),
0555         dim_(r.dim_),
0556         allConstInterpolated_(r.allConstInterpolated_) {
0557     for (unsigned i = 0; i < dim_; ++i) {
0558       leftInterpolationLinear_[i] = r.leftInterpolationLinear_[i];
0559       rightInterpolationLinear_[i] = r.rightInterpolationLinear_[i];
0560     }
0561   }
0562 
0563   template <typename Numeric, class Axis>
0564   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(const ArrayND<Numeric>& data,
0565                                                                 const std::vector<Axis>& axes,
0566                                                                 const char* leftInterpolation,
0567                                                                 const char* rightInterpolation,
0568                                                                 const std::string& label)
0569       : data_(data), axes_(axes), functionLabel_(label), dim_(data.rank()) {
0570     for (unsigned i = 0; i < dim_; ++i) {
0571       leftInterpolationLinear_[i] = leftInterpolation[i];
0572       rightInterpolationLinear_[i] = rightInterpolation[i];
0573     }
0574     allConstInterpolated_ = allConstInterpolated();
0575   }
0576 
0577   template <typename Numeric, class Axis>
0578   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(const Axis& xAxis,
0579                                                                 bool leftX,
0580                                                                 bool rightX,
0581                                                                 const Axis& yAxis,
0582                                                                 bool leftY,
0583                                                                 bool rightY,
0584                                                                 const Axis& zAxis,
0585                                                                 bool leftZ,
0586                                                                 bool rightZ,
0587                                                                 const Axis& tAxis,
0588                                                                 bool leftT,
0589                                                                 bool rightT,
0590                                                                 const Axis& vAxis,
0591                                                                 bool leftV,
0592                                                                 bool rightV,
0593                                                                 const char* label)
0594       : data_(Private::makeTableShape(xAxis, yAxis, zAxis, tAxis, vAxis)),
0595         functionLabel_(label ? label : ""),
0596         dim_(5U) {
0597     axes_.reserve(dim_);
0598     axes_.push_back(xAxis);
0599     axes_.push_back(yAxis);
0600     axes_.push_back(zAxis);
0601     axes_.push_back(tAxis);
0602     axes_.push_back(vAxis);
0603 
0604     unsigned i = 0;
0605     leftInterpolationLinear_[i] = leftX;
0606     rightInterpolationLinear_[i++] = rightX;
0607     leftInterpolationLinear_[i] = leftY;
0608     rightInterpolationLinear_[i++] = rightY;
0609     leftInterpolationLinear_[i] = leftZ;
0610     rightInterpolationLinear_[i++] = rightZ;
0611     leftInterpolationLinear_[i] = leftT;
0612     rightInterpolationLinear_[i++] = rightT;
0613     leftInterpolationLinear_[i] = leftV;
0614     rightInterpolationLinear_[i++] = rightV;
0615     assert(i == dim_);
0616 
0617     allConstInterpolated_ = allConstInterpolated();
0618   }
0619 
0620   template <typename Numeric, class Axis>
0621   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(const Axis& xAxis,
0622                                                                 bool leftX,
0623                                                                 bool rightX,
0624                                                                 const Axis& yAxis,
0625                                                                 bool leftY,
0626                                                                 bool rightY,
0627                                                                 const Axis& zAxis,
0628                                                                 bool leftZ,
0629                                                                 bool rightZ,
0630                                                                 const Axis& tAxis,
0631                                                                 bool leftT,
0632                                                                 bool rightT,
0633                                                                 const char* label)
0634       : data_(Private::makeTableShape(xAxis, yAxis, zAxis, tAxis)), functionLabel_(label ? label : ""), dim_(4U) {
0635     axes_.reserve(dim_);
0636     axes_.push_back(xAxis);
0637     axes_.push_back(yAxis);
0638     axes_.push_back(zAxis);
0639     axes_.push_back(tAxis);
0640 
0641     unsigned i = 0;
0642     leftInterpolationLinear_[i] = leftX;
0643     rightInterpolationLinear_[i++] = rightX;
0644     leftInterpolationLinear_[i] = leftY;
0645     rightInterpolationLinear_[i++] = rightY;
0646     leftInterpolationLinear_[i] = leftZ;
0647     rightInterpolationLinear_[i++] = rightZ;
0648     leftInterpolationLinear_[i] = leftT;
0649     rightInterpolationLinear_[i++] = rightT;
0650     assert(i == dim_);
0651 
0652     allConstInterpolated_ = allConstInterpolated();
0653   }
0654 
0655   template <typename Numeric, class Axis>
0656   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(const Axis& xAxis,
0657                                                                 bool leftX,
0658                                                                 bool rightX,
0659                                                                 const Axis& yAxis,
0660                                                                 bool leftY,
0661                                                                 bool rightY,
0662                                                                 const Axis& zAxis,
0663                                                                 bool leftZ,
0664                                                                 bool rightZ,
0665                                                                 const char* label)
0666       : data_(Private::makeTableShape(xAxis, yAxis, zAxis)), functionLabel_(label ? label : ""), dim_(3U) {
0667     axes_.reserve(dim_);
0668     axes_.push_back(xAxis);
0669     axes_.push_back(yAxis);
0670     axes_.push_back(zAxis);
0671 
0672     unsigned i = 0;
0673     leftInterpolationLinear_[i] = leftX;
0674     rightInterpolationLinear_[i++] = rightX;
0675     leftInterpolationLinear_[i] = leftY;
0676     rightInterpolationLinear_[i++] = rightY;
0677     leftInterpolationLinear_[i] = leftZ;
0678     rightInterpolationLinear_[i++] = rightZ;
0679     assert(i == dim_);
0680 
0681     allConstInterpolated_ = allConstInterpolated();
0682   }
0683 
0684   template <typename Numeric, class Axis>
0685   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(
0686       const Axis& xAxis, bool leftX, bool rightX, const Axis& yAxis, bool leftY, bool rightY, const char* label)
0687       : data_(Private::makeTableShape(xAxis, yAxis)), functionLabel_(label ? label : ""), dim_(2U) {
0688     axes_.reserve(dim_);
0689     axes_.push_back(xAxis);
0690     axes_.push_back(yAxis);
0691 
0692     unsigned i = 0;
0693     leftInterpolationLinear_[i] = leftX;
0694     rightInterpolationLinear_[i++] = rightX;
0695     leftInterpolationLinear_[i] = leftY;
0696     rightInterpolationLinear_[i++] = rightY;
0697     assert(i == dim_);
0698 
0699     allConstInterpolated_ = allConstInterpolated();
0700   }
0701 
0702   template <typename Numeric, class Axis>
0703   LinInterpolatedTableND<Numeric, Axis>::LinInterpolatedTableND(const Axis& xAxis,
0704                                                                 bool leftX,
0705                                                                 bool rightX,
0706                                                                 const char* label)
0707       : data_(Private::makeTableShape(xAxis)), functionLabel_(label ? label : ""), dim_(1U) {
0708     axes_.reserve(dim_);
0709     axes_.push_back(xAxis);
0710 
0711     leftInterpolationLinear_[0] = leftX;
0712     rightInterpolationLinear_[0] = rightX;
0713 
0714     allConstInterpolated_ = allConstInterpolated();
0715   }
0716 
0717   template <typename Numeric, class Axis>
0718   template <typename ConvertibleToUnsigned>
0719   std::unique_ptr<LinInterpolatedTableND<Numeric, Axis> > LinInterpolatedTableND<Numeric, Axis>::invertWRTAxis(
0720       const ConvertibleToUnsigned axisNumC,
0721       const Axis& replacementAxis,
0722       const bool leftLinear,
0723       const bool rightLinear,
0724       const char* functionLabel) const {
0725     const unsigned axisNumber = static_cast<unsigned>(axisNumC);
0726 
0727     if (axisNumber >= dim_)
0728       throw npstat::NpstatOutOfRange(
0729           "In npstat::LinInterpolatedTableND::invertAxis: "
0730           "axis number is out of range");
0731 
0732     // Generate the new set of axes
0733     std::vector<Axis> newAxes(axes_);
0734     newAxes[axisNumber] = replacementAxis;
0735 
0736     std::vector<std::pair<bool, bool> > iType(interpolationType());
0737     iType[axisNumber] = std::pair<bool, bool>(leftLinear, rightLinear);
0738 
0739     // Create the new table
0740     std::unique_ptr<LinInterpolatedTableND> pTable(new LinInterpolatedTableND(newAxes, iType, functionLabel));
0741 
0742     if (dim_ > 1U) {
0743       // Prepare array slices
0744       unsigned sliceIndex[CHAR_BIT * sizeof(unsigned long)];
0745       unsigned fixedIndices[CHAR_BIT * sizeof(unsigned long)];
0746       unsigned count = 0;
0747       for (unsigned i = 0; i < dim_; ++i)
0748         if (i != axisNumber) {
0749           sliceIndex[count] = data_.span(i);
0750           fixedIndices[count++] = i;
0751         }
0752       ArrayND<Numeric> parentSlice(data_, fixedIndices, count);
0753       ArrayND<Numeric> dauSlice(pTable->data_, fixedIndices, count);
0754 
0755       // Cycle over the slices
0756       for (ArrayNDScanner scan(sliceIndex, count); scan.isValid(); ++scan) {
0757         scan.getIndex(sliceIndex, count);
0758         data_.exportSlice(&parentSlice, fixedIndices, sliceIndex, count);
0759         Private::lind_invert1DSlice(
0760             parentSlice, axes_[axisNumber], replacementAxis, leftLinear, rightLinear, &dauSlice);
0761         pTable->data_.importSlice(dauSlice, fixedIndices, sliceIndex, count);
0762       }
0763     } else
0764       Private::lind_invert1DSlice(data_, axes_[0], replacementAxis, leftLinear, rightLinear, &pTable->data_);
0765     return pTable;
0766   }
0767 
0768   template <typename Numeric, class Axis>
0769   template <class Functor1, class Functor2>
0770   std::unique_ptr<LinInterpolatedTableND<Numeric, Axis> > LinInterpolatedTableND<Numeric, Axis>::invertRatioResponse(
0771       const unsigned axisNumber,
0772       const Axis& replacementAxis,
0773       const bool leftLinear,
0774       const bool rightLinear,
0775       Functor1 invg,
0776       Functor2 invh,
0777       const char* functionLabel) const {
0778     if (axisNumber >= dim_)
0779       throw npstat::NpstatOutOfRange(
0780           "In npstat::LinInterpolatedTableND::invertRatioResponse: "
0781           "axis number is out of range");
0782 
0783     // Generate the new set of axes
0784     std::vector<Axis> newAxes(axes_);
0785     newAxes[axisNumber] = replacementAxis;
0786 
0787     std::vector<std::pair<bool, bool> > iType(interpolationType());
0788     iType[axisNumber] = std::pair<bool, bool>(leftLinear, rightLinear);
0789 
0790     // Transform the original axis to the raw x values
0791     const Axis& oldAxis(axes_[axisNumber]);
0792     std::vector<double> rawx;
0793     const unsigned nCoords = oldAxis.nCoords();
0794     rawx.reserve(nCoords);
0795     for (unsigned i = 0; i < nCoords; ++i) {
0796       const double x = invg(oldAxis.coordinate(i));
0797       if (x < 0.0)
0798         throw npstat::NpstatInvalidArgument(
0799             "In npstat::LinInterpolatedTableND::invertRatioResponse: "
0800             "invalid original axis definition (negative transformed "
0801             "coordinate)");
0802       rawx.push_back(x);
0803     }
0804 
0805     // Transform the new axis to the raw f(x) values
0806     std::vector<double> rawf;
0807     const unsigned nFuncs = replacementAxis.nCoords();
0808     rawf.reserve(nFuncs);
0809     for (unsigned i = 0; i < nFuncs; ++i) {
0810       const double f = invh(replacementAxis.coordinate(i));
0811       if (f < 0.0)
0812         throw npstat::NpstatInvalidArgument(
0813             "In npstat::LinInterpolatedTableND::invertRatioResponse: "
0814             "invalid new axis definition (negative transformed "
0815             "coordinate)");
0816       rawf.push_back(f);
0817     }
0818 
0819     // Workspace needed for the inversion code
0820     std::vector<double> workspace(nCoords);
0821 
0822     // Create the new table
0823     std::unique_ptr<LinInterpolatedTableND> pTable(new LinInterpolatedTableND(newAxes, iType, functionLabel));
0824 
0825     if (dim_ > 1U) {
0826       // Prepare array slices
0827       unsigned sliceIndex[CHAR_BIT * sizeof(unsigned long)];
0828       unsigned fixedIndices[CHAR_BIT * sizeof(unsigned long)];
0829       unsigned count = 0;
0830       for (unsigned i = 0; i < dim_; ++i)
0831         if (i != axisNumber) {
0832           sliceIndex[count] = data_.span(i);
0833           fixedIndices[count++] = i;
0834         }
0835       ArrayND<Numeric> parentSlice(data_, fixedIndices, count);
0836       ArrayND<Numeric> dauSlice(pTable->data_, fixedIndices, count);
0837 
0838       // Cycle over the slices
0839       for (ArrayNDScanner scan(sliceIndex, count); scan.isValid(); ++scan) {
0840         scan.getIndex(sliceIndex, count);
0841         data_.exportSlice(&parentSlice, fixedIndices, sliceIndex, count);
0842         invert1DResponse(parentSlice,
0843                          oldAxis,
0844                          replacementAxis,
0845                          leftLinear,
0846                          rightLinear,
0847                          invg,
0848                          &rawx[0],
0849                          &rawf[0],
0850                          &workspace[0],
0851                          &dauSlice);
0852         pTable->data_.importSlice(dauSlice, fixedIndices, sliceIndex, count);
0853       }
0854     } else
0855       invert1DResponse(data_,
0856                        oldAxis,
0857                        replacementAxis,
0858                        leftLinear,
0859                        rightLinear,
0860                        invg,
0861                        &rawx[0],
0862                        &rawf[0],
0863                        &workspace[0],
0864                        &pTable->data_);
0865     return pTable;
0866   }
0867 
0868   template <typename Numeric, class Axis>
0869   void LinInterpolatedTableND<Numeric, Axis>::getCoords(const unsigned long linearIndex,
0870                                                         double* coords,
0871                                                         const unsigned coordsBufferSize) const {
0872     if (coordsBufferSize < dim_)
0873       throw npstat::NpstatInvalidArgument(
0874           "In LinInterpolatedTableND::getCoords: "
0875           "insufficient buffer size");
0876     assert(coords);
0877     unsigned index[CHAR_BIT * sizeof(unsigned long)];
0878     data_.convertLinearIndex(linearIndex, index, dim_);
0879     for (unsigned i = 0; i < dim_; ++i)
0880       coords[i] = axes_[i].coordinate(index[i]);
0881   }
0882 
0883   template <typename Numeric, class Axis>
0884   bool LinInterpolatedTableND<Numeric, Axis>::isWithinLimits(const double* point, const unsigned len) const {
0885     if (len != dim_)
0886       throw npstat::NpstatInvalidArgument(
0887           "In npstat::LinInterpolatedTableND::isWithinLimits: "
0888           "incompatible point dimensionality");
0889     assert(point);
0890 
0891     for (unsigned i = 0; i < dim_; ++i)
0892       if (point[i] < axes_[i].min() || point[i] > axes_[i].max())
0893         return false;
0894     return true;
0895   }
0896 
0897   template <typename Numeric, class Axis>
0898   Numeric LinInterpolatedTableND<Numeric, Axis>::operator()(const double* point, const unsigned len) const {
0899     typedef typename ProperDblFromCmpl<Numeric>::type proper_double;
0900 
0901     if (len != dim_)
0902       throw npstat::NpstatInvalidArgument(
0903           "In npstat::LinInterpolatedTableND::operator(): "
0904           "incompatible point dimensionality");
0905     assert(point);
0906 
0907     bool interpolateArray = true;
0908     if (!allConstInterpolated_)
0909       for (unsigned i = 0; i < dim_; ++i)
0910         if ((leftInterpolationLinear_[i] && point[i] < axes_[i].min()) ||
0911             (rightInterpolationLinear_[i] && point[i] > axes_[i].max())) {
0912           interpolateArray = false;
0913           break;
0914         }
0915 
0916     if (interpolateArray) {
0917       // Translate coordinates into the array system and
0918       // simply use the ArrayND interpolation facilities
0919       double buf[CHAR_BIT * sizeof(unsigned long)];
0920       for (unsigned i = 0; i < dim_; ++i) {
0921         const std::pair<unsigned, double>& pair = axes_[i].getInterval(point[i]);
0922         buf[i] = pair.first + 1U - pair.second;
0923       }
0924       return data_.interpolate1(buf, dim_);
0925     } else {
0926       unsigned ix[CHAR_BIT * sizeof(unsigned long)];
0927       double weight[CHAR_BIT * sizeof(unsigned long)];
0928       for (unsigned i = 0; i < dim_; ++i) {
0929         const bool linear = (leftInterpolationLinear_[i] && point[i] < axes_[i].min()) ||
0930                             (rightInterpolationLinear_[i] && point[i] > axes_[i].max());
0931         const std::pair<unsigned, double>& pair =
0932             linear ? axes_[i].linearInterval(point[i]) : axes_[i].getInterval(point[i]);
0933         ix[i] = pair.first;
0934         weight[i] = pair.second;
0935       }
0936 
0937       Numeric sum = Numeric();
0938       const unsigned long maxcycle = 1UL << dim_;
0939       const unsigned long* strides = data_.strides();
0940       const Numeric* dat = data_.data();
0941       for (unsigned long icycle = 0UL; icycle < maxcycle; ++icycle) {
0942         double w = 1.0;
0943         unsigned long icell = 0UL;
0944         for (unsigned i = 0; i < dim_; ++i) {
0945           if (icycle & (1UL << i)) {
0946             w *= (1.0 - weight[i]);
0947             icell += strides[i] * (ix[i] + 1U);
0948           } else {
0949             w *= weight[i];
0950             icell += strides[i] * ix[i];
0951           }
0952         }
0953         sum += dat[icell] * static_cast<proper_double>(w);
0954       }
0955       return sum;
0956     }
0957   }
0958 
0959   template <typename Numeric, class Axis>
0960   Numeric LinInterpolatedTableND<Numeric, Axis>::operator()(const double& x0) const {
0961     const unsigned nArgs = 1U;
0962     if (dim_ != nArgs)
0963       throw npstat::NpstatInvalidArgument(
0964           "In npstat::LinInterpolatedTableND::operator(): number of "
0965           "arguments, 1, is incompatible with the interpolator dimensionality");
0966     double tmp[nArgs];
0967     tmp[0] = x0;
0968     return operator()(tmp, nArgs);
0969   }
0970 
0971   template <typename Numeric, class Axis>
0972   Numeric LinInterpolatedTableND<Numeric, Axis>::operator()(const double& x0, const double& x1) const {
0973     const unsigned nArgs = 2U;
0974     if (dim_ != nArgs)
0975       throw npstat::NpstatInvalidArgument(
0976           "In npstat::LinInterpolatedTableND::operator(): number of "
0977           "arguments, 2, is incompatible with the interpolator dimensionality");
0978     double tmp[nArgs];
0979     tmp[0] = x0;
0980     tmp[1] = x1;
0981     return operator()(tmp, nArgs);
0982   }
0983 
0984   template <typename Numeric, class Axis>
0985   Numeric LinInterpolatedTableND<Numeric, Axis>::operator()(const double& x0,
0986                                                             const double& x1,
0987                                                             const double& x2) const {
0988     const unsigned nArgs = 3U;
0989     if (dim_ != nArgs)
0990       throw npstat::NpstatInvalidArgument(
0991           "In npstat::LinInterpolatedTableND::operator(): number of "
0992           "arguments, 3, is incompatible with the interpolator dimensionality");
0993     double tmp[nArgs];
0994     tmp[0] = x0;
0995     tmp[1] = x1;
0996     tmp[2] = x2;
0997     return operator()(tmp, nArgs);
0998   }
0999 
1000   template <typename Numeric, class Axis>
1001   Numeric LinInterpolatedTableND<Numeric, Axis>::operator()(const double& x0,
1002                                                             const double& x1,
1003                                                             const double& x2,
1004                                                             const double& x3) const {
1005     const unsigned nArgs = 4U;
1006     if (dim_ != nArgs)
1007       throw npstat::NpstatInvalidArgument(
1008           "In npstat::LinInterpolatedTableND::operator(): number of "
1009           "arguments, 4, is incompatible with the interpolator dimensionality");
1010     double tmp[nArgs];
1011     tmp[0] = x0;
1012     tmp[1] = x1;
1013     tmp[2] = x2;
1014     tmp[3] = x3;
1015     return operator()(tmp, nArgs);
1016   }
1017 
1018   template <typename Numeric, class Axis>
1019   Numeric LinInterpolatedTableND<Numeric, Axis>::operator()(
1020       const double& x0, const double& x1, const double& x2, const double& x3, const double& x4) const {
1021     const unsigned nArgs = 5U;
1022     if (dim_ != nArgs)
1023       throw npstat::NpstatInvalidArgument(
1024           "In npstat::LinInterpolatedTableND::operator(): number of "
1025           "arguments, 5, is incompatible with the interpolator dimensionality");
1026     double tmp[nArgs];
1027     tmp[0] = x0;
1028     tmp[1] = x1;
1029     tmp[2] = x2;
1030     tmp[3] = x3;
1031     tmp[4] = x4;
1032     return operator()(tmp, nArgs);
1033   }
1034 
1035   template <typename Numeric, class Axis>
1036   template <class Functor1>
1037   double LinInterpolatedTableND<Numeric, Axis>::solveForRatioArg(
1038       const double xmin, const double xmax, const double rmin, const double rmax, const double fval, Functor1 invg) {
1039     // Find two values of x so that f(x0) <= fval <= f(x1)
1040     double x0 = xmin;
1041     double x1 = xmax;
1042     double fmin = invg(xmin) * rmin;
1043     double fmax = invg(xmax) * rmax;
1044     const double step = xmax - xmin;
1045     assert(fmin < fmax);
1046     assert(step > 0.0);
1047 
1048     unsigned stepcount = 0;
1049     const unsigned maxSteps = 1000U;
1050     for (double stepfactor = 1.0; (fval < fmin || fval > fmax) && stepcount < maxSteps; stepfactor *= 2.0, ++stepcount)
1051       if (fval < fmin) {
1052         x1 = x0;
1053         fmax = fmin;
1054         x0 -= stepfactor * step;
1055         fmin = invg(x0) * Private::lind_interpolateSimple(xmin, xmax, rmin, rmax, x0);
1056       } else {
1057         x0 = x1;
1058         fmin = fmax;
1059         x1 += stepfactor * step;
1060         fmax = invg(x1) * Private::lind_interpolateSimple(xmin, xmax, rmin, rmax, x1);
1061       }
1062     if (stepcount == maxSteps)
1063       throw npstat::NpstatRuntimeError(
1064           "In LinInterpolatedTableND::solveForRatioArg: "
1065           "faled to bracket the root");
1066 
1067     assert(x1 >= x0);
1068     while ((x1 - x0) / (std::abs(x1) + std::abs(x0) + DBL_EPSILON) > 4.0 * DBL_EPSILON) {
1069       const double xhalf = (x1 + x0) / 2.0;
1070       const double fhalf = invg(xhalf) * Private::lind_interpolateSimple(xmin, xmax, rmin, rmax, xhalf);
1071       if (fval < fhalf) {
1072         x1 = xhalf;
1073         fmax = fhalf;
1074       } else {
1075         x0 = xhalf;
1076         fmin = fhalf;
1077       }
1078     }
1079     return (x1 + x0) / 2.0;
1080   }
1081 
1082   template <typename Numeric, class Axis>
1083   template <class Functor1>
1084   void LinInterpolatedTableND<Numeric, Axis>::invert1DResponse(const ArrayND<Numeric>& fromSlice,
1085                                                                const Axis& fromAxis,
1086                                                                const Axis& toAxis,
1087                                                                const bool newLeftLinear,
1088                                                                const bool newRightLinear,
1089                                                                Functor1 invg,
1090                                                                const double* rawx,
1091                                                                const double* rawf,
1092                                                                double* workspace,
1093                                                                ArrayND<Numeric>* toSlice) {
1094     assert(toSlice);
1095     assert(fromSlice.rank() == 1U);
1096     assert(toSlice->rank() == 1U);
1097 
1098     const Numeric zero = Numeric();
1099     const Numeric* fromData = fromSlice.data();
1100     const unsigned fromLen = fromSlice.length();
1101     assert(fromLen > 1U);
1102     assert(fromLen == fromAxis.nCoords());
1103     Numeric* toD = const_cast<Numeric*>(toSlice->data());
1104     const unsigned nAxisPoints = toAxis.nCoords();
1105     assert(toSlice->length() == nAxisPoints);
1106 
1107     for (unsigned i = 0; i < fromLen; ++i) {
1108       if (fromData[i] <= zero)
1109         throw npstat::NpstatDomainError(
1110             "In LinInterpolatedTableND::invert1DResponse: "
1111             "non-positive response found. This ratio "
1112             "response table is not invertible.");
1113       workspace[i] = rawx[i] * fromData[i];
1114     }
1115 
1116     const double yfirst = workspace[0];
1117     const double ylast = workspace[fromLen - 1U];
1118 
1119     bool adjustZero = false;
1120     unsigned nBelow = 0;
1121     for (unsigned ipt = 0; ipt < nAxisPoints; ++ipt) {
1122       const double y = rawf[ipt];
1123       unsigned i0 = 0;
1124       bool solve = false;
1125       if (y == 0.0) {
1126         assert(ipt == 0U);
1127         if (newLeftLinear)
1128           adjustZero = true;
1129       } else if (y <= yfirst) {
1130         ++nBelow;
1131         solve = newLeftLinear;
1132       } else if (y >= ylast) {
1133         solve = newRightLinear;
1134         i0 = solve ? fromLen - 2 : fromLen - 1;
1135       } else {
1136         solve = true;
1137         i0 = static_cast<unsigned>(std::lower_bound(workspace, workspace + fromLen, y) - workspace) - 1U;
1138       }
1139       if (solve) {
1140         const double x = solveForRatioArg(
1141             fromAxis.coordinate(i0), fromAxis.coordinate(i0 + 1), fromData[i0], fromData[i0 + 1], y, invg);
1142         toD[ipt] = invg(x) / y;
1143       } else
1144         toD[ipt] = 1.0 / fromData[i0];
1145     }
1146     if (adjustZero && nBelow)
1147       toD[0] = toD[1];
1148   }
1149 }  // namespace npstat
1150 
1151 #endif  // NPSTAT_LININTERPOLATEDTABLEND_HH_
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!