npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
LinInterpolatedTableND.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_LININTERPOLATEDTABLEND_HH_
2 #define NPSTAT_LININTERPOLATEDTABLEND_HH_
3 
4 /**
5 // \file LinInterpolatedTableND.hh
6 //
7 // \brief Multilinear interpolation/extrapolation on rectangular grids
8 //
9 // Author: I. Volobouev
10 //
11 // June 2012
12 */
13 
14 #include <climits>
15 #include <vector>
16 #include <utility>
17 
18 #include "geners/CPP11_auto_ptr.hh"
19 
20 #include "npstat/nm/ArrayND.hh"
21 #include "npstat/nm/UniformAxis.hh"
22 
23 namespace npstat {
24  /**
25  // Template for multilinear interpolation/extrapolation of values provided
26  // on a rectangular coordinate grid. "Numeric" is the type stored in
27  // the table. "Axis" should be one of GridAxis, UniformAxis, or DualAxis
28  // classes or a user-provided class with a similar set of methods.
29  */
30  template <class Numeric, class Axis=UniformAxis>
32  {
33  template <typename Num2, typename Axis2>
34  friend class LinInterpolatedTableND;
35 
36  public:
37  typedef Numeric value_type;
38  typedef Axis axis_type;
39 
40  /**
41  // Main constructor for arbitrary-dimensional interpolators.
42  //
43  // "axes" are the axes of the rectangular coordinate grid.
44  //
45  // In each pair provided by the "extrapolationType" argument,
46  // the first element of the pair specifies whether the extrapolation
47  // to negative infinity should be linear (if "true") or constant
48  // (if "false"). The second element of the pair specifies whether
49  // to extrapolate linearly to positive infinity.
50  //
51  // "functionLabel" is an arbitrary string which can potentially
52  // be used by plotting programs and such.
53  */
55  const std::vector<Axis>& axes,
56  const std::vector<std::pair<bool,bool> >& extrapolationType,
57  const char* functionLabel=0);
58 
59  /** Convenience constructor for 1-d interpolators */
60  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
61  const char* functionLabel=0);
62 
63  /** Convenience constructor for 2-d interpolators */
64  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
65  const Axis& yAxis, bool leftY, bool rightY,
66  const char* functionLabel=0);
67 
68  /** Convenience constructor for 3-d interpolators */
69  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
70  const Axis& yAxis, bool leftY, bool rightY,
71  const Axis& zAxis, bool leftZ, bool rightZ,
72  const char* functionLabel=0);
73 
74  /** Convenience constructor for 4-d interpolators */
75  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
76  const Axis& yAxis, bool leftY, bool rightY,
77  const Axis& zAxis, bool leftZ, bool rightZ,
78  const Axis& tAxis, bool leftT, bool rightT,
79  const char* functionLabel=0);
80 
81  /** Convenience constructor for 5-d interpolators */
82  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
83  const Axis& yAxis, bool leftY, bool rightY,
84  const Axis& zAxis, bool leftZ, bool rightZ,
85  const Axis& tAxis, bool leftT, bool rightT,
86  const Axis& vAxis, bool leftV, bool rightV,
87  const char* functionLabel=0);
88 
89  /**
90  // Converting copy constructor from interpolator
91  // with another storage type
92  */
93  template <class Num2>
95 
96  /**
97  // Basic interpolation result. Argument point dimensionality must be
98  // compatible with the interpolator dimensionality.
99  */
100  Numeric operator()(const double* point, unsigned dim) const;
101 
102  //@{
103  /** Convenience function for low-dimensional interpolators */
104  Numeric operator()(const double& x0) const;
105  Numeric operator()(const double& x0, const double& x1) const;
106  Numeric operator()(const double& x0, const double& x1,
107  const double& x2) const;
108  Numeric operator()(const double& x0, const double& x1,
109  const double& x2, const double& x3) const;
110  Numeric operator()(const double& x0, const double& x1,
111  const double& x2, const double& x3,
112  const double& x4) const;
113  //@}
114 
115  //@{
116  /** Examine interpolator contents */
117  inline unsigned dim() const {return dim_;}
118  inline const std::vector<Axis>& axes() const {return axes_;}
119  inline const Axis& axis(const unsigned i) const
120  {return axes_.at(i);}
121  inline unsigned long length() const {return data_.length();}
122  bool leftInterpolationLinear(unsigned i) const;
123  bool rightInterpolationLinear(unsigned i) const;
124  std::vector<std::pair<bool,bool> > interpolationType() const;
125  inline const std::string& functionLabel() const
126  {return functionLabel_;}
127  //@}
128 
129  //@{
130  /** Access the interpolator table data */
131  inline const ArrayND<Numeric>& table() const {return data_;}
132  inline ArrayND<Numeric>& table() {return data_;}
133  //@}
134 
135  /** Convenience function for getting coordinates of the grid points */
136  void getCoords(unsigned long linearIndex,
137  double* coords, unsigned coordsBufferSize) const;
138 
139  /**
140  // This method returns "true" if the method isUniform()
141  // of each interpolator axis returns "true"
142  */
143  bool isUniformlyBinned() const;
144 
145  /**
146  // This method will return "true" if the point
147  // is inside the grid limits or on the boundary
148  */
149  bool isWithinLimits(const double* point, unsigned dim) const;
150 
151  /** Modifier for the function label */
152  inline void setFunctionLabel(const char* newlabel)
153  {functionLabel_ = newlabel ? newlabel : "";}
154 
155  /**
156  // Invert the function w.r.t. the variable represented by one of
157  // the axes. The function values must change monotonously along
158  // the chosen axis. Note that this operation is meaningful only
159  // in case "Numeric" type is either float or double.
160  */
161  template <typename ConvertibleToUnsigned>
162  CPP11_auto_ptr<LinInterpolatedTableND> invertWRTAxis(
163  ConvertibleToUnsigned axisNumber, const Axis& replacementAxis,
164  bool newAxisLeftLinear, bool newAxisRightLinear,
165  const char* functionLabel=0) const;
166 
167  /**
168  // This method inverts the ratio response.
169  // That is, we assume that the table encodes r(x) for
170  // some function f(x) = x*r(x). We also assume that the
171  // original axis does not represent x directly -- instead,
172  // axis coordinates are given by g(x) (in practice, g is
173  // often the natural log). We will also assume that the new
174  // axis is not going to represent f(x) directly -- it
175  // will be h(f(x)) instead. The functors "invg" and "invh"
176  // below must do the inverse: taking the axes coordinates
177  // to the actual values of x and f(x). Both "invg" and "invh"
178  // must be monotonously increasing functions. The code assumes
179  // that x*r(x) -> 0 when x->0 (that is, r(x) is bounded at 0)
180  // and it also assumes (but does not check explicitly)
181  // that x*r(x) is monotonously increasing with x.
182  //
183  // The returned interpolation table encodes the values
184  // of x/f(x). Of course, they are just 1/r(x), but the trick
185  // is to be able to look them up quickly as a function of
186  // h(f(x)). Naturally, this whole operation is meaningful
187  // only in case "Numeric" type is either float or double.
188  */
189  template <class Functor1, class Functor2>
190  CPP11_auto_ptr<LinInterpolatedTableND> invertRatioResponse(
191  unsigned axisNumber, const Axis& replacementAxis,
192  bool newAxisLeftLinear, bool newAxisRightLinear,
193  Functor1 invg, Functor2 invh,
194  const char* functionLabel=0) const;
195 
196  /** Comparison for equality */
197  bool operator==(const LinInterpolatedTableND&) const;
198 
199  /** Logical negation of operator== */
200  inline bool operator!=(const LinInterpolatedTableND& r) const
201  {return !(*this == r);}
202 
203  //@{
204  // Method related to "geners" I/O
205  inline gs::ClassId classId() const {return gs::ClassId(*this);}
206  bool write(std::ostream& of) const;
207  //@}
208 
209  static const char* classname();
210  static inline unsigned version() {return 1;}
211  static LinInterpolatedTableND* read(
212  const gs::ClassId& id, std::istream& in);
213 
214  private:
215  LinInterpolatedTableND();
216 
217  LinInterpolatedTableND(
218  const ArrayND<Numeric>& data,
219  const std::vector<Axis>& axes,
220  const char* leftInterpolation,
221  const char* rightInterpolation,
222  const std::string& label);
223 
224  bool allConstInterpolated() const;
225 
226  ArrayND<Numeric> data_;
227  std::vector<Axis> axes_;
228  std::string functionLabel_;
229  char leftInterpolationLinear_[CHAR_BIT*sizeof(unsigned long)];
230  char rightInterpolationLinear_[CHAR_BIT*sizeof(unsigned long)];
231  unsigned dim_;
232  bool allConstInterpolated_;
233 
234  template <class Functor1>
235  static double solveForRatioArg(double xmin, double xmax,
236  double rmin, double rmax,
237  double fval, Functor1 invg);
238 
239  template <class Functor1>
240  static void invert1DResponse(const ArrayND<Numeric>& fromSlice,
241  const Axis& fromAxis, const Axis& toAxis,
242  bool newLeftLinear, bool newRightLinear,
243  Functor1 invg,
244  const double* rawx, const double* rawf,
245  double* workspace,
246  ArrayND<Numeric>* toSlice);
247  };
248 }
249 
250 #include "npstat/nm/LinInterpolatedTableND.icc"
251 
252 #endif // NPSTAT_LININTERPOLATEDTABLEND_HH_
Arbitrary-dimensional array template.
Uniformly spaced coordinate sets for use in constructing rectangular grids.
unsigned long length() const
Definition: ArrayND.hh:320
Definition: LinInterpolatedTableND.hh:32
unsigned dim() const
Definition: LinInterpolatedTableND.hh:117
LinInterpolatedTableND(const Axis &xAxis, bool leftX, bool rightX, const Axis &yAxis, bool leftY, bool rightY, const char *functionLabel=0)
void setFunctionLabel(const char *newlabel)
Definition: LinInterpolatedTableND.hh:152
bool operator!=(const LinInterpolatedTableND &r) const
Definition: LinInterpolatedTableND.hh:200
LinInterpolatedTableND(const Axis &xAxis, bool leftX, bool rightX, const Axis &yAxis, bool leftY, bool rightY, const Axis &zAxis, bool leftZ, bool rightZ, const char *functionLabel=0)
LinInterpolatedTableND(const Axis &xAxis, bool leftX, bool rightX, const Axis &yAxis, bool leftY, bool rightY, const Axis &zAxis, bool leftZ, bool rightZ, const Axis &tAxis, bool leftT, bool rightT, const Axis &vAxis, bool leftV, bool rightV, const char *functionLabel=0)
LinInterpolatedTableND(const Axis &xAxis, bool leftX, bool rightX, const char *functionLabel=0)
void getCoords(unsigned long linearIndex, double *coords, unsigned coordsBufferSize) const
const ArrayND< Numeric > & table() const
Definition: LinInterpolatedTableND.hh:131
CPP11_auto_ptr< LinInterpolatedTableND > invertRatioResponse(unsigned axisNumber, const Axis &replacementAxis, bool newAxisLeftLinear, bool newAxisRightLinear, Functor1 invg, Functor2 invh, const char *functionLabel=0) const
CPP11_auto_ptr< LinInterpolatedTableND > invertWRTAxis(ConvertibleToUnsigned axisNumber, const Axis &replacementAxis, bool newAxisLeftLinear, bool newAxisRightLinear, const char *functionLabel=0) const
bool operator==(const LinInterpolatedTableND &) const
bool isWithinLimits(const double *point, unsigned dim) const
Numeric operator()(const double *point, unsigned dim) const
LinInterpolatedTableND(const std::vector< Axis > &axes, const std::vector< std::pair< bool, bool > > &extrapolationType, const char *functionLabel=0)
LinInterpolatedTableND(const LinInterpolatedTableND< Num2, Axis > &)
Numeric operator()(const double &x0) const
LinInterpolatedTableND(const Axis &xAxis, bool leftX, bool rightX, const Axis &yAxis, bool leftY, bool rightY, const Axis &zAxis, bool leftZ, bool rightZ, const Axis &tAxis, bool leftT, bool rightT, const char *functionLabel=0)
Definition: AbsArrayProjector.hh:14
Definition: SimpleFunctors.hh:58
Definition: SimpleFunctors.hh:89