CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
LinInterpolatedTableND.h
Go to the documentation of this file.
1 #ifndef NPSTAT_LININTERPOLATEDTABLEND_HH_
2 #define NPSTAT_LININTERPOLATEDTABLEND_HH_
3 
14 #include <climits>
15 #include <vector>
16 #include <utility>
17 
18 #include "Alignment/Geners/interface/CPP11_auto_ptr.hh"
19 
22 
23 namespace npstat {
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 
55  const std::vector<Axis>& axes,
56  const std::vector<std::pair<bool,bool> >& extrapolationType,
57  const char* functionLabel=0);
58 
60  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
61  const char* functionLabel=0);
62 
64  LinInterpolatedTableND(const Axis& xAxis, bool leftX, bool rightX,
65  const Axis& yAxis, bool leftY, bool rightY,
66  const char* functionLabel=0);
67 
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 
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 
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 
93  template <class Num2>
95 
100  Numeric operator()(const double* point, unsigned dim) const;
101 
103 
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;
114 
116 
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_;}
128 
130 
131  inline const ArrayND<Numeric>& table() const {return data_;}
132  inline ArrayND<Numeric>& table() {return data_;}
134 
136  void getCoords(unsigned long linearIndex,
137  double* coords, unsigned coordsBufferSize) const;
138 
143  bool isUniformlyBinned() const;
144 
149  bool isWithinLimits(const double* point, unsigned dim) const;
150 
152  inline void setFunctionLabel(const char* newlabel)
153  {functionLabel_ = newlabel ? newlabel : "";}
154 
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 
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 
197  bool operator==(const LinInterpolatedTableND&) const;
198 
200  inline bool operator!=(const LinInterpolatedTableND& r) const
201  {return !(*this == r);}
202 
204  // Method related to "geners" I/O
205  inline gs::ClassId classId() const {return gs::ClassId(*this);}
206  bool write(std::ostream& of) const;
208 
209  static const char* classname();
210  static inline unsigned version() {return 1;}
212  const gs::ClassId& id, std::istream& in);
213 
214  private:
216 
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 
227  std::vector<Axis> axes_;
229  char leftInterpolationLinear_[CHAR_BIT*sizeof(unsigned long)];
230  char rightInterpolationLinear_[CHAR_BIT*sizeof(unsigned long)];
231  unsigned dim_;
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 <cmath>
251 #include <cfloat>
252 #include <cassert>
253 #include <algorithm>
254 #include <functional>
255 
256 #include "Alignment/Geners/interface/binaryIO.hh"
257 
260 
261 namespace npstat {
262  namespace Private {
263  template <class Axis>
264  ArrayShape makeTableShape(const std::vector<Axis>& axes)
265  {
266  const unsigned n = axes.size();
268  result.reserve(n);
269  for (unsigned i=0; i<n; ++i)
270  result.push_back(axes[i].nCoords());
271  return result;
272  }
273 
274  template <class Axis>
275  ArrayShape makeTableShape(const Axis& xAxis)
276  {
278  result.reserve(1U);
279  result.push_back(xAxis.nCoords());
280  return result;
281  }
282 
283  template <class Axis>
284  ArrayShape makeTableShape(const Axis& xAxis, const Axis& yAxis)
285  {
287  result.reserve(2U);
288  result.push_back(xAxis.nCoords());
289  result.push_back(yAxis.nCoords());
290  return result;
291  }
292 
293  template <class Axis>
294  ArrayShape makeTableShape(const Axis& xAxis,
295  const Axis& yAxis,
296  const Axis& zAxis)
297  {
299  result.reserve(3U);
300  result.push_back(xAxis.nCoords());
301  result.push_back(yAxis.nCoords());
302  result.push_back(zAxis.nCoords());
303  return result;
304  }
305 
306  template <class Axis>
307  ArrayShape makeTableShape(const Axis& xAxis, const Axis& yAxis,
308  const Axis& zAxis, const Axis& tAxis)
309  {
311  result.reserve(4U);
312  result.push_back(xAxis.nCoords());
313  result.push_back(yAxis.nCoords());
314  result.push_back(zAxis.nCoords());
315  result.push_back(tAxis.nCoords());
316  return result;
317  }
318 
319  template <class Axis>
320  ArrayShape makeTableShape(const Axis& xAxis, const Axis& yAxis,
321  const Axis& zAxis, const Axis& tAxis,
322  const Axis& vAxis)
323  {
325  result.reserve(5U);
326  result.push_back(xAxis.nCoords());
327  result.push_back(yAxis.nCoords());
328  result.push_back(zAxis.nCoords());
329  result.push_back(tAxis.nCoords());
330  result.push_back(vAxis.nCoords());
331  return result;
332  }
333 
334  inline double lind_interpolateSimple(
335  const double x0, const double x1,
336  const double y0, const double y1,
337  const double x)
338  {
339  return y0 + (y1 - y0)*((x - x0)/(x1 - x0));
340  }
341 
342  template <typename Numeric, class Axis>
344  const ArrayND<Numeric>& fromSlice,
345  const Axis& fromAxis, const Axis& toAxis,
346  const bool leftLinear, const bool rightLinear,
347  ArrayND<Numeric>* toSlice)
348  {
349  assert(toSlice);
350  assert(fromSlice.rank() == 1U);
351  assert(toSlice->rank() == 1U);
352 
353  const Numeric* fromData = fromSlice.data();
354  const unsigned fromLen = fromSlice.length();
355  assert(fromLen > 1U);
356  assert(fromLen == fromAxis.nCoords());
357  const Numeric* fromDataEnd = fromData + fromLen;
358  if (!isStrictlyMonotonous(fromData, fromDataEnd))
360  "In npstat::Private::lind_invert1DSlice: "
361  "slice data is not monotonous and can not be inverted");
362 
363  const Numeric yfirst = fromData[0];
364  const Numeric ylast = fromData[fromLen - 1U];
365  const bool increasing = yfirst < ylast;
366 
367  Numeric* toD = const_cast<Numeric*>(toSlice->data());
368  const unsigned nAxisPoints = toAxis.nCoords();
369  assert(toSlice->length() == nAxisPoints);
370 
371  for (unsigned ipt=0; ipt<nAxisPoints; ++ipt)
372  {
373  const Numeric y = static_cast<Numeric>(toAxis.coordinate(ipt));
374  if (increasing)
375  {
376  if (y <= yfirst)
377  {
378  if (leftLinear)
380  yfirst, fromData[1], fromAxis.coordinate(0),
381  fromAxis.coordinate(1), y);
382  else
383  toD[ipt] = fromAxis.coordinate(0);
384  }
385  else if (y >= ylast)
386  {
387  if (rightLinear)
389  ylast, fromData[fromLen - 2U],
390  fromAxis.coordinate(fromLen - 1U),
391  fromAxis.coordinate(fromLen - 2U), y);
392  else
393  toD[ipt] = fromAxis.coordinate(fromLen - 1U);
394  }
395  else
396  {
397  const unsigned i = std::lower_bound(fromData,fromDataEnd,y)-
398  fromData;
400  fromData[i-1U], fromData[i],
401  fromAxis.coordinate(i-1U),
402  fromAxis.coordinate(i), y);
403  }
404  }
405  else
406  {
407  // The role of left and right are exchanged
408  // with respect to first and last point
409  if (y <= ylast)
410  {
411  if (leftLinear)
413  ylast, fromData[fromLen - 2U],
414  fromAxis.coordinate(fromLen - 1U),
415  fromAxis.coordinate(fromLen - 2U), y);
416  else
417  toD[ipt] = fromAxis.coordinate(fromLen - 1U);
418  }
419  else if (y >= yfirst)
420  {
421  if (rightLinear)
423  yfirst, fromData[1],
424  fromAxis.coordinate(0),
425  fromAxis.coordinate(1), y);
426  else
427  toD[ipt] = fromAxis.coordinate(0);
428  }
429  else
430  {
431  const unsigned i = std::lower_bound(fromData,fromDataEnd,
432  y,std::greater<Numeric>())-fromData;
434  fromData[i-1U], fromData[i],
435  fromAxis.coordinate(i-1U),
436  fromAxis.coordinate(i), y);
437  }
438  }
439  }
440  }
441  }
442 
443  template <class Numeric, class Axis>
445  {
446  for (unsigned i=0; i<dim_; ++i)
447  if (leftInterpolationLinear_[i] || rightInterpolationLinear_[i])
448  return false;
449  return true;
450  }
451 
452  template <class Numeric, class Axis>
454  const LinInterpolatedTableND& r) const
455  {
456  if (dim_ != r.dim_)
457  return false;
458  for (unsigned i=0; i<dim_; ++i)
459  {
460  if (leftInterpolationLinear_[i] != r.leftInterpolationLinear_[i])
461  return false;
462  if (rightInterpolationLinear_[i] != r.rightInterpolationLinear_[i])
463  return false;
464  }
465  return data_ == r.data_ &&
466  axes_ == r.axes_ &&
467  functionLabel_ == r.functionLabel_;
468  }
469 
470  template <typename Numeric, class Axis>
472  {
473  static const std::string myClass(gs::template_class_name<Numeric,Axis>(
474  "npstat::LinInterpolatedTableND"));
475  return myClass.c_str();
476  }
477 
478  template<typename Numeric, class Axis>
479  bool LinInterpolatedTableND<Numeric,Axis>::write(std::ostream& of) const
480  {
481  const bool status = data_.classId().write(of) &&
482  data_.write(of) &&
483  gs::write_obj_vector(of, axes_);
484  if (status)
485  {
486  gs::write_pod_array(of, leftInterpolationLinear_, dim_);
487  gs::write_pod_array(of, rightInterpolationLinear_, dim_);
488  gs::write_pod(of, functionLabel_);
489  }
490  return status && !of.fail();
491  }
492 
493  template<typename Numeric, class Axis>
496  const gs::ClassId& id, std::istream& in)
497  {
498  static const gs::ClassId current(
499  gs::ClassId::makeId<LinInterpolatedTableND<Numeric,Axis> >());
500  current.ensureSameId(id);
501 
502  gs::ClassId ida(in, 1);
504  ArrayND<Numeric>::restore(ida, in, &data);
505  std::vector<Axis> axes;
506  gs::read_heap_obj_vector_as_placed(in, &axes);
507  const unsigned dim = axes.size();
508  if (dim > CHAR_BIT*sizeof(unsigned long) || data.rank() != dim)
509  throw gs::IOInvalidData(
510  "In npstat::LinInterpolatedTableND::read: "
511  "read back invalid dimensionality");
512  char leftInterpolation[CHAR_BIT*sizeof(unsigned long)];
513  gs::read_pod_array(in, leftInterpolation, dim);
514  char rightInterpolation[CHAR_BIT*sizeof(unsigned long)];
515  gs::read_pod_array(in, rightInterpolation, dim);
517  gs::read_pod(in, &label);
518  if (in.fail()) throw gs::IOReadFailure(
519  "In npstat::LinInterpolatedTableND::read: input stream failure");
520  return new LinInterpolatedTableND(
521  data, axes, leftInterpolation, rightInterpolation, label);
522  }
523 
524  template<typename Numeric, class Axis>
526  const unsigned i) const
527  {
528  if (i >= dim_) throw npstat::NpstatOutOfRange(
529  "In npstat::LinInterpolatedTableND::leftInterpolationLinear: "
530  "index out of range");
531  return leftInterpolationLinear_[i];
532  }
533 
534  template<typename Numeric, class Axis>
536  const unsigned i) const
537  {
538  if (i >= dim_) throw npstat::NpstatOutOfRange(
539  "In npstat::LinInterpolatedTableND::rightInterpolationLinear: "
540  "index out of range");
541  return rightInterpolationLinear_[i];
542  }
543 
544  template<typename Numeric, class Axis>
546  {
547  for (unsigned i=0; i<dim_; ++i)
548  if (!axes_[i].isUniform())
549  return false;
550  return true;
551  }
552 
553  template<typename Numeric, class Axis>
554  std::vector<std::pair<bool,bool> >
556  {
557  std::vector<std::pair<bool,bool> > vec;
558  vec.reserve(dim_);
559  for (unsigned i=0; i<dim_; ++i)
560  vec.push_back(std::pair<bool, bool>(leftInterpolationLinear_[i],
561  rightInterpolationLinear_[i]));
562  return vec;
563  }
564 
565  template<typename Numeric, class Axis>
567  const std::vector<Axis>& axes,
568  const std::vector<std::pair<bool,bool> >& interpolationType,
569  const char* label)
570  : data_(Private::makeTableShape(axes)),
571  axes_(axes),
572  functionLabel_(label ? label : ""),
573  dim_(axes.size())
574  {
575  if (dim_ == 0 || dim_ >= CHAR_BIT*sizeof(unsigned long))
577  "In npstat::LinInterpolatedTableND constructor: requested "
578  "table dimensionality is not supported");
579  if (dim_ != interpolationType.size())
581  "In npstat::LinInterpolatedTableND constructor: "
582  "incompatible number of interpolation specifications");
583  for (unsigned i=0; i<dim_; ++i)
584  {
585  const std::pair<bool,bool>& pair(interpolationType[i]);
586  leftInterpolationLinear_[i] = pair.first;
587  rightInterpolationLinear_[i] = pair.second;
588  }
589 
591  }
592 
593  template<typename Numeric, class Axis>
594  template <class Num2>
597  : data_(r.data_),
598  axes_(r.axes_),
599  functionLabel_(r.functionLabel_),
600  dim_(r.dim_),
601  allConstInterpolated_(r.allConstInterpolated_)
602  {
603  for (unsigned i=0; i<dim_; ++i)
604  {
607  }
608  }
609 
610  template<typename Numeric, class Axis>
612  const ArrayND<Numeric>& data,
613  const std::vector<Axis>& axes,
614  const char* leftInterpolation,
615  const char* rightInterpolation,
616  const std::string& label)
617  : data_(data),
618  axes_(axes),
619  functionLabel_(label),
620  dim_(data.rank())
621  {
622  for (unsigned i=0; i<dim_; ++i)
623  {
624  leftInterpolationLinear_[i] = leftInterpolation[i];
625  rightInterpolationLinear_[i] = rightInterpolation[i];
626  }
628  }
629 
630  template<typename Numeric, class Axis>
632  const Axis& xAxis, bool leftX, bool rightX,
633  const Axis& yAxis, bool leftY, bool rightY,
634  const Axis& zAxis, bool leftZ, bool rightZ,
635  const Axis& tAxis, bool leftT, bool rightT,
636  const Axis& vAxis, bool leftV, bool rightV,
637  const char* label)
638  : data_(Private::makeTableShape(xAxis, yAxis, zAxis, tAxis, vAxis)),
639  functionLabel_(label ? label : ""),
640  dim_(5U)
641  {
642  axes_.reserve(dim_);
643  axes_.push_back(xAxis);
644  axes_.push_back(yAxis);
645  axes_.push_back(zAxis);
646  axes_.push_back(tAxis);
647  axes_.push_back(vAxis);
648 
649  unsigned i = 0;
650  leftInterpolationLinear_[i] = leftX;
651  rightInterpolationLinear_[i++] = rightX;
652  leftInterpolationLinear_[i] = leftY;
653  rightInterpolationLinear_[i++] = rightY;
654  leftInterpolationLinear_[i] = leftZ;
655  rightInterpolationLinear_[i++] = rightZ;
656  leftInterpolationLinear_[i] = leftT;
657  rightInterpolationLinear_[i++] = rightT;
658  leftInterpolationLinear_[i] = leftV;
659  rightInterpolationLinear_[i++] = rightV;
660  assert(i == dim_);
661 
663  }
664 
665  template<typename Numeric, class Axis>
667  const Axis& xAxis, bool leftX, bool rightX,
668  const Axis& yAxis, bool leftY, bool rightY,
669  const Axis& zAxis, bool leftZ, bool rightZ,
670  const Axis& tAxis, bool leftT, bool rightT,
671  const char* label)
672  : data_(Private::makeTableShape(xAxis, yAxis, zAxis, tAxis)),
673  functionLabel_(label ? label : ""),
674  dim_(4U)
675  {
676  axes_.reserve(dim_);
677  axes_.push_back(xAxis);
678  axes_.push_back(yAxis);
679  axes_.push_back(zAxis);
680  axes_.push_back(tAxis);
681 
682  unsigned i = 0;
683  leftInterpolationLinear_[i] = leftX;
684  rightInterpolationLinear_[i++] = rightX;
685  leftInterpolationLinear_[i] = leftY;
686  rightInterpolationLinear_[i++] = rightY;
687  leftInterpolationLinear_[i] = leftZ;
688  rightInterpolationLinear_[i++] = rightZ;
689  leftInterpolationLinear_[i] = leftT;
690  rightInterpolationLinear_[i++] = rightT;
691  assert(i == dim_);
692 
694  }
695 
696  template<typename Numeric, class Axis>
698  const Axis& xAxis, bool leftX, bool rightX,
699  const Axis& yAxis, bool leftY, bool rightY,
700  const Axis& zAxis, bool leftZ, bool rightZ,
701  const char* label)
702  : data_(Private::makeTableShape(xAxis, yAxis, zAxis)),
703  functionLabel_(label ? label : ""),
704  dim_(3U)
705  {
706  axes_.reserve(dim_);
707  axes_.push_back(xAxis);
708  axes_.push_back(yAxis);
709  axes_.push_back(zAxis);
710 
711  unsigned i = 0;
712  leftInterpolationLinear_[i] = leftX;
713  rightInterpolationLinear_[i++] = rightX;
714  leftInterpolationLinear_[i] = leftY;
715  rightInterpolationLinear_[i++] = rightY;
716  leftInterpolationLinear_[i] = leftZ;
717  rightInterpolationLinear_[i++] = rightZ;
718  assert(i == dim_);
719 
721  }
722 
723  template<typename Numeric, class Axis>
725  const Axis& xAxis, bool leftX, bool rightX,
726  const Axis& yAxis, bool leftY, bool rightY,
727  const char* label)
728  : data_(Private::makeTableShape(xAxis, yAxis)),
729  functionLabel_(label ? label : ""),
730  dim_(2U)
731  {
732  axes_.reserve(dim_);
733  axes_.push_back(xAxis);
734  axes_.push_back(yAxis);
735 
736  unsigned i = 0;
737  leftInterpolationLinear_[i] = leftX;
738  rightInterpolationLinear_[i++] = rightX;
739  leftInterpolationLinear_[i] = leftY;
740  rightInterpolationLinear_[i++] = rightY;
741  assert(i == dim_);
742 
744  }
745 
746  template<typename Numeric, class Axis>
748  const Axis& xAxis, bool leftX, bool rightX,
749  const char* label)
750  : data_(Private::makeTableShape(xAxis)),
751  functionLabel_(label ? label : ""),
752  dim_(1U)
753  {
754  axes_.reserve(dim_);
755  axes_.push_back(xAxis);
756 
757  leftInterpolationLinear_[0] = leftX;
758  rightInterpolationLinear_[0] = rightX;
759 
761  }
762 
763  template<typename Numeric, class Axis>
764  template <typename ConvertibleToUnsigned>
765  CPP11_auto_ptr<LinInterpolatedTableND<Numeric,Axis> >
767  const ConvertibleToUnsigned axisNumC, const Axis& replacementAxis,
768  const bool leftLinear, const bool rightLinear,
769  const char* functionLabel) const
770  {
771  const unsigned axisNumber = static_cast<unsigned>(axisNumC);
772 
773  if (axisNumber >= dim_)
775  "In npstat::LinInterpolatedTableND::invertAxis: "
776  "axis number is out of range");
777 
778  // Generate the new set of axes
779  std::vector<Axis> newAxes(axes_);
780  newAxes[axisNumber] = replacementAxis;
781 
782  std::vector<std::pair<bool,bool> > iType(interpolationType());
783  iType[axisNumber] = std::pair<bool,bool>(leftLinear, rightLinear);
784 
785  // Create the new table
786  CPP11_auto_ptr<LinInterpolatedTableND> pTable(
787  new LinInterpolatedTableND(newAxes, iType, functionLabel));
788 
789  if (dim_ > 1U)
790  {
791  // Prepare array slices
792  unsigned sliceIndex[CHAR_BIT*sizeof(unsigned long)];
793  unsigned fixedIndices[CHAR_BIT*sizeof(unsigned long)];
794  unsigned count = 0;
795  for (unsigned i=0; i<dim_; ++i)
796  if (i != axisNumber)
797  {
798  sliceIndex[count] = data_.span(i);
799  fixedIndices[count++] = i;
800  }
801  ArrayND<Numeric> parentSlice(data_, fixedIndices, count);
802  ArrayND<Numeric> dauSlice(pTable->data_, fixedIndices, count);
803 
804  // Cycle over the slices
805  for (ArrayNDScanner scan(sliceIndex,count); scan.isValid(); ++scan)
806  {
807  scan.getIndex(sliceIndex, count);
808  data_.exportSlice(&parentSlice, fixedIndices,
809  sliceIndex, count);
811  parentSlice, axes_[axisNumber], replacementAxis,
812  leftLinear, rightLinear, &dauSlice);
813  pTable->data_.importSlice(dauSlice, fixedIndices,
814  sliceIndex, count);
815  }
816  }
817  else
819  data_, axes_[0], replacementAxis,
820  leftLinear, rightLinear, &pTable->data_);
821  return pTable;
822  }
823 
824  template<typename Numeric, class Axis>
825  template <class Functor1, class Functor2>
826  CPP11_auto_ptr<LinInterpolatedTableND<Numeric,Axis> >
828  const unsigned axisNumber, const Axis& replacementAxis,
829  const bool leftLinear, const bool rightLinear,
830  Functor1 invg, Functor2 invh,
831  const char* functionLabel) const
832  {
833  if (axisNumber >= dim_)
835  "In npstat::LinInterpolatedTableND::invertRatioResponse: "
836  "axis number is out of range");
837 
838  // Generate the new set of axes
839  std::vector<Axis> newAxes(axes_);
840  newAxes[axisNumber] = replacementAxis;
841 
842  std::vector<std::pair<bool,bool> > iType(interpolationType());
843  iType[axisNumber] = std::pair<bool,bool>(leftLinear, rightLinear);
844 
845  // Transform the original axis to the raw x values
846  const Axis& oldAxis(axes_[axisNumber]);
847  std::vector<double> rawx;
848  const unsigned nCoords = oldAxis.nCoords();
849  rawx.reserve(nCoords);
850  for (unsigned i=0; i<nCoords; ++i)
851  {
852  const double x = invg(oldAxis.coordinate(i));
853  if (x < 0.0)
855  "In npstat::LinInterpolatedTableND::invertRatioResponse: "
856  "invalid original axis definition (negative transformed "
857  "coordinate)");
858  rawx.push_back(x);
859  }
860 
861  // Transform the new axis to the raw f(x) values
862  std::vector<double> rawf;
863  const unsigned nFuncs = replacementAxis.nCoords();
864  rawf.reserve(nFuncs);
865  for (unsigned i=0; i<nFuncs; ++i)
866  {
867  const double f = invh(replacementAxis.coordinate(i));
868  if (f < 0.0)
870  "In npstat::LinInterpolatedTableND::invertRatioResponse: "
871  "invalid new axis definition (negative transformed "
872  "coordinate)");
873  rawf.push_back(f);
874  }
875 
876  // Workspace needed for the inversion code
877  std::vector<double> workspace(nCoords);
878 
879  // Create the new table
880  CPP11_auto_ptr<LinInterpolatedTableND> pTable(
881  new LinInterpolatedTableND(newAxes, iType, functionLabel));
882 
883  if (dim_ > 1U)
884  {
885  // Prepare array slices
886  unsigned sliceIndex[CHAR_BIT*sizeof(unsigned long)];
887  unsigned fixedIndices[CHAR_BIT*sizeof(unsigned long)];
888  unsigned count = 0;
889  for (unsigned i=0; i<dim_; ++i)
890  if (i != axisNumber)
891  {
892  sliceIndex[count] = data_.span(i);
893  fixedIndices[count++] = i;
894  }
895  ArrayND<Numeric> parentSlice(data_, fixedIndices, count);
896  ArrayND<Numeric> dauSlice(pTable->data_, fixedIndices, count);
897 
898  // Cycle over the slices
899  for (ArrayNDScanner scan(sliceIndex,count); scan.isValid(); ++scan)
900  {
901  scan.getIndex(sliceIndex, count);
902  data_.exportSlice(&parentSlice, fixedIndices,
903  sliceIndex, count);
904  invert1DResponse(parentSlice, oldAxis,
905  replacementAxis, leftLinear, rightLinear,
906  invg, &rawx[0], &rawf[0], &workspace[0],
907  &dauSlice);
908  pTable->data_.importSlice(dauSlice, fixedIndices,
909  sliceIndex, count);
910  }
911  }
912  else
913  invert1DResponse(data_, oldAxis, replacementAxis, leftLinear,
914  rightLinear, invg, &rawx[0], &rawf[0],
915  &workspace[0], &pTable->data_);
916  return pTable;
917  }
918 
919  template<typename Numeric, class Axis>
921  const unsigned long linearIndex,
922  double* coords, const unsigned coordsBufferSize) const
923  {
924  if (coordsBufferSize < dim_) throw npstat::NpstatInvalidArgument(
925  "In LinInterpolatedTableND::getCoords: "
926  "insufficient buffer size");
927  assert(coords);
928  unsigned index[CHAR_BIT*sizeof(unsigned long)];
929  data_.convertLinearIndex(linearIndex, index, dim_);
930  for (unsigned i=0; i<dim_; ++i)
931  coords[i] = axes_[i].coordinate(index[i]);
932  }
933 
934  template<typename Numeric, class Axis>
936  const double* point, const unsigned len) const
937  {
938  if (len != dim_)
940  "In npstat::LinInterpolatedTableND::isWithinLimits: "
941  "incompatible point dimensionality");
942  assert(point);
943 
944  for (unsigned i=0; i<dim_; ++i)
945  if (point[i] < axes_[i].min() || point[i] > axes_[i].max())
946  return false;
947  return true;
948  }
949 
950  template<typename Numeric, class Axis>
952  const double* point, const unsigned len) const
953  {
954  typedef typename ProperDblFromCmpl<Numeric>::type proper_double;
955 
956  if (len != dim_)
958  "In npstat::LinInterpolatedTableND::operator(): "
959  "incompatible point dimensionality");
960  assert(point);
961 
962  bool interpolateArray = true;
963  if (!allConstInterpolated_)
964  for (unsigned i=0; i<dim_; ++i)
965  if ((leftInterpolationLinear_[i] && point[i] < axes_[i].min()) ||
966  (rightInterpolationLinear_[i] && point[i] > axes_[i].max()))
967  {
968  interpolateArray = false;
969  break;
970  }
971 
972  if (interpolateArray)
973  {
974  // Translate coordinates into the array system and
975  // simply use the ArrayND interpolation facilities
976  double buf[CHAR_BIT*sizeof(unsigned long)];
977  for (unsigned i=0; i<dim_; ++i)
978  {
979  const std::pair<unsigned,double>& pair =
980  axes_[i].getInterval(point[i]);
981  buf[i] = pair.first + 1U - pair.second;
982  }
983  return data_.interpolate1(buf, dim_);
984  }
985  else
986  {
987  unsigned ix[CHAR_BIT*sizeof(unsigned long)];
988  double weight[CHAR_BIT*sizeof(unsigned long)];
989  for (unsigned i=0; i<dim_; ++i)
990  {
991  const bool linear = (leftInterpolationLinear_[i] &&
992  point[i] < axes_[i].min()) ||
993  (rightInterpolationLinear_[i] &&
994  point[i] > axes_[i].max());
995  const std::pair<unsigned,double>& pair = linear ?
996  axes_[i].linearInterval(point[i]) :
997  axes_[i].getInterval(point[i]);
998  ix[i] = pair.first;
999  weight[i] = pair.second;
1000  }
1001 
1002  Numeric sum = Numeric();
1003  const unsigned long maxcycle = 1UL << dim_;
1004  const unsigned long* strides = data_.strides();
1005  const Numeric* dat = data_.data();
1006  for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
1007  {
1008  double w = 1.0;
1009  unsigned long icell = 0UL;
1010  for (unsigned i=0; i<dim_; ++i)
1011  {
1012  if (icycle & (1UL << i))
1013  {
1014  w *= (1.0 - weight[i]);
1015  icell += strides[i]*(ix[i] + 1U);
1016  }
1017  else
1018  {
1019  w *= weight[i];
1020  icell += strides[i]*ix[i];
1021  }
1022  }
1023  sum += dat[icell]*static_cast<proper_double>(w);
1024  }
1025  return sum;
1026  }
1027  }
1028 
1029  template<typename Numeric, class Axis>
1031  const double& x0) const
1032  {
1033  const unsigned nArgs = 1U;
1034  if (dim_ != nArgs) throw npstat::NpstatInvalidArgument(
1035  "In npstat::LinInterpolatedTableND::operator(): number of "
1036  "arguments, 1, is incompatible with the interpolator dimensionality");
1037  double tmp[nArgs];
1038  tmp[0] = x0;
1039  return operator()(tmp, nArgs);
1040  }
1041 
1042  template<typename Numeric, class Axis>
1044  const double& x0, const double& x1) const
1045  {
1046  const unsigned nArgs = 2U;
1047  if (dim_ != nArgs) throw npstat::NpstatInvalidArgument(
1048  "In npstat::LinInterpolatedTableND::operator(): number of "
1049  "arguments, 2, is incompatible with the interpolator dimensionality");
1050  double tmp[nArgs];
1051  tmp[0] = x0;
1052  tmp[1] = x1;
1053  return operator()(tmp, nArgs);
1054  }
1055 
1056  template<typename Numeric, class Axis>
1058  const double& x0, const double& x1, const double& x2) const
1059  {
1060  const unsigned nArgs = 3U;
1061  if (dim_ != nArgs) throw npstat::NpstatInvalidArgument(
1062  "In npstat::LinInterpolatedTableND::operator(): number of "
1063  "arguments, 3, is incompatible with the interpolator dimensionality");
1064  double tmp[nArgs];
1065  tmp[0] = x0;
1066  tmp[1] = x1;
1067  tmp[2] = x2;
1068  return operator()(tmp, nArgs);
1069  }
1070 
1071  template<typename Numeric, class Axis>
1073  const double& x0, const double& x1,
1074  const double& x2, const double& x3) const
1075  {
1076  const unsigned nArgs = 4U;
1077  if (dim_ != nArgs) throw npstat::NpstatInvalidArgument(
1078  "In npstat::LinInterpolatedTableND::operator(): number of "
1079  "arguments, 4, is incompatible with the interpolator dimensionality");
1080  double tmp[nArgs];
1081  tmp[0] = x0;
1082  tmp[1] = x1;
1083  tmp[2] = x2;
1084  tmp[3] = x3;
1085  return operator()(tmp, nArgs);
1086  }
1087 
1088  template<typename Numeric, class Axis>
1090  const double& x0, const double& x1, const double& x2,
1091  const double& x3, const double& x4) const
1092  {
1093  const unsigned nArgs = 5U;
1094  if (dim_ != nArgs) throw npstat::NpstatInvalidArgument(
1095  "In npstat::LinInterpolatedTableND::operator(): number of "
1096  "arguments, 5, is incompatible with the interpolator dimensionality");
1097  double tmp[nArgs];
1098  tmp[0] = x0;
1099  tmp[1] = x1;
1100  tmp[2] = x2;
1101  tmp[3] = x3;
1102  tmp[4] = x4;
1103  return operator()(tmp, nArgs);
1104  }
1105 
1106  template<typename Numeric, class Axis>
1107  template <class Functor1>
1109  const double xmin, const double xmax,
1110  const double rmin, const double rmax,
1111  const double fval, Functor1 invg)
1112  {
1113  // Find two values of x so that f(x0) <= fval <= f(x1)
1114  double x0 = xmin;
1115  double x1 = xmax;
1116  double fmin = invg(xmin)*rmin;
1117  double fmax = invg(xmax)*rmax;
1118  const double step = xmax - xmin;
1119  assert(fmin < fmax);
1120  assert(step > 0.0);
1121 
1122  unsigned stepcount = 0;
1123  const unsigned maxSteps = 1000U;
1124  for (double stepfactor = 1.0; (fval < fmin || fval > fmax) &&
1125  stepcount < maxSteps; stepfactor *= 2.0, ++stepcount)
1126  if (fval < fmin)
1127  {
1128  x1 = x0;
1129  fmax = fmin;
1130  x0 -= stepfactor*step;
1131  fmin = invg(x0)*Private::lind_interpolateSimple(
1132  xmin, xmax, rmin, rmax, x0);
1133  }
1134  else
1135  {
1136  x0 = x1;
1137  fmin = fmax;
1138  x1 += stepfactor*step;
1139  fmax = invg(x1)*Private::lind_interpolateSimple(
1140  xmin, xmax, rmin, rmax, x1);
1141  }
1142  if (stepcount == maxSteps) throw npstat::NpstatRuntimeError(
1143  "In LinInterpolatedTableND::solveForRatioArg: "
1144  "faled to bracket the root");
1145 
1146  assert(x1 >= x0);
1147  while ((x1 - x0)/(std::abs(x1) + std::abs(x0) + DBL_EPSILON) > 4.0*DBL_EPSILON)
1148  {
1149  const double xhalf = (x1 + x0)/2.0;
1150  const double fhalf = invg(xhalf)*Private::lind_interpolateSimple(
1151  xmin, xmax, rmin, rmax, xhalf);
1152  if (fval < fhalf)
1153  {
1154  x1 = xhalf;
1155  fmax = fhalf;
1156  }
1157  else
1158  {
1159  x0 = xhalf;
1160  fmin = fhalf;
1161  }
1162  }
1163  return (x1 + x0)/2.0;
1164  }
1165 
1166  template<typename Numeric, class Axis>
1167  template <class Functor1>
1169  const ArrayND<Numeric>& fromSlice,
1170  const Axis& fromAxis, const Axis& toAxis,
1171  const bool newLeftLinear, const bool newRightLinear,
1172  Functor1 invg, const double* rawx, const double* rawf,
1173  double* workspace,
1174  ArrayND<Numeric>* toSlice)
1175  {
1176  assert(toSlice);
1177  assert(fromSlice.rank() == 1U);
1178  assert(toSlice->rank() == 1U);
1179 
1180  const Numeric zero = Numeric();
1181  const Numeric* fromData = fromSlice.data();
1182  const unsigned fromLen = fromSlice.length();
1183  assert(fromLen > 1U);
1184  assert(fromLen == fromAxis.nCoords());
1185  Numeric* toD = const_cast<Numeric*>(toSlice->data());
1186  const unsigned nAxisPoints = toAxis.nCoords();
1187  assert(toSlice->length() == nAxisPoints);
1188 
1189  for (unsigned i=0; i<fromLen; ++i)
1190  {
1191  if (fromData[i] <= zero) throw npstat::NpstatDomainError(
1192  "In LinInterpolatedTableND::invert1DResponse: "
1193  "non-positive response found. This ratio "
1194  "response table is not invertible.");
1195  workspace[i] = rawx[i]*fromData[i];
1196  }
1197 
1198  const double yfirst = workspace[0];
1199  const double ylast = workspace[fromLen - 1U];
1200 
1201  bool adjustZero = false;
1202  unsigned nBelow = 0;
1203  for (unsigned ipt=0; ipt<nAxisPoints; ++ipt)
1204  {
1205  const double y = rawf[ipt];
1206  unsigned i0 = 0;
1207  bool solve = false;
1208  if (y == 0.0)
1209  {
1210  assert(ipt == 0U);
1211  if (newLeftLinear)
1212  adjustZero = true;
1213  }
1214  else if (y <= yfirst)
1215  {
1216  ++nBelow;
1217  solve = newLeftLinear;
1218  }
1219  else if (y >= ylast)
1220  {
1221  solve = newRightLinear;
1222  i0 = solve ? fromLen-2 : fromLen-1;
1223  }
1224  else
1225  {
1226  solve = true;
1227  i0 = static_cast<unsigned>(std::lower_bound(
1228  workspace,workspace+fromLen,y) - workspace) - 1U;
1229  }
1230  if (solve)
1231  {
1232  const double x = solveForRatioArg(fromAxis.coordinate(i0),
1233  fromAxis.coordinate(i0+1),
1234  fromData[i0], fromData[i0+1],
1235  y, invg);
1236  toD[ipt] = invg(x)/y;
1237  }
1238  else
1239  toD[ipt] = 1.0/fromData[i0];
1240  }
1241  if (adjustZero && nBelow)
1242  toD[0] = toD[1];
1243  }
1244 }
1245 
1246 
1247 #endif // NPSTAT_LININTERPOLATEDTABLEND_HH_
1248 
bool rightInterpolationLinear(unsigned i) const
int i
Definition: DBlmapReader.cc:9
Uniformly spaced coordinate sets for use in constructing rectangular grids.
static void invert1DResponse(const ArrayND< Numeric > &fromSlice, const Axis &fromAxis, const Axis &toAxis, bool newLeftLinear, bool newRightLinear, Functor1 invg, const double *rawx, const double *rawf, double *workspace, ArrayND< Numeric > *toSlice)
const double w
Definition: UKUtility.cc:23
unsigned rank() const
Definition: ArrayND.h:240
assert(m_qm.get())
void getCoords(unsigned long linearIndex, double *coords, unsigned coordsBufferSize) const
std::unique_ptr< LinInterpolatedTableND > invertRatioResponse(unsigned axisNumber, const Axis &replacementAxis, bool newAxisLeftLinear, bool newAxisRightLinear, Functor1 invg, Functor2 invh, const char *functionLabel=0) const
std::vector< unsigned > ArrayShape
Definition: ArrayShape.h:21
unsigned long length() const
Definition: ArrayND.h:231
tuple result
Definition: mps_fire.py:84
const ArrayND< Numeric > & table() const
static double solveForRatioArg(double xmin, double xmax, double rmin, double rmax, double fval, Functor1 invg)
T x() const
Cartesian x coordinate.
Numeric operator()(const double *point, unsigned dim) const
bool write(std::ostream &of) const
A few simple template functions for checking monotonicity of container values.
bool operator==(const LinInterpolatedTableND &) const
std::vector< std::pair< bool, bool > > interpolationType() const
const std::vector< Axis > & axes() const
static LinInterpolatedTableND * read(const gs::ClassId &id, std::istream &in)
bool operator!=(const LinInterpolatedTableND &r) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const Numeric * data() const
Definition: ArrayND.h:234
double f[11][100]
void setFunctionLabel(const char *newlabel)
T min(T a, T b)
Definition: MathUtil.h:58
char rightInterpolationLinear_[CHAR_BIT *sizeof(unsigned long)]
std::unique_ptr< LinInterpolatedTableND > invertWRTAxis(ConvertibleToUnsigned axisNumber, const Axis &replacementAxis, bool newAxisLeftLinear, bool newAxisRightLinear, const char *functionLabel=0) const
char leftInterpolationLinear_[CHAR_BIT *sizeof(unsigned long)]
bool leftInterpolationLinear(unsigned i) const
double lind_interpolateSimple(const double x0, const double x1, const double y0, const double y1, const double x)
bool isStrictlyMonotonous(Iter const begin, Iter const end)
Definition: isMonotonous.h:46
static void restore(const gs::ClassId &id, std::istream &in, ArrayND *array)
Definition: ArrayND.h:6052
ArrayShape makeTableShape(const std::vector< Axis > &axes)
void lind_invert1DSlice(const ArrayND< Numeric > &fromSlice, const Axis &fromAxis, const Axis &toAxis, const bool leftLinear, const bool rightLinear, ArrayND< Numeric > *toSlice)
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
bool isWithinLimits(const double *point, unsigned dim) const
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple zAxis
Definition: MetAnalyzer.py:56
float linear(float x)
const Axis & axis(const unsigned i) const
int weight
Definition: histoStyle.py:50
const std::string & functionLabel() const
Iteration over indices of a multidimensional array.
tuple size
Write out results.
Arbitrary-dimensional array template.
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
tuple status
Definition: mps_update.py:57