1 #ifndef NPSTAT_LININTERPOLATEDTABLEND_HH_
2 #define NPSTAT_LININTERPOLATEDTABLEND_HH_
18 #include "Alignment/Geners/interface/CPP11_auto_ptr.hh"
30 template <
class Numeric,
class Axis=UniformAxis>
33 template <
typename Num2,
typename Axis2>
55 const std::vector<Axis>&
axes,
56 const std::vector<std::pair<bool,bool> >& extrapolationType,
65 const Axis& yAxis,
bool leftY,
bool rightY,
70 const Axis& yAxis,
bool leftY,
bool rightY,
71 const Axis& zAxis,
bool leftZ,
bool rightZ,
76 const Axis& yAxis,
bool leftY,
bool rightY,
77 const Axis& zAxis,
bool leftZ,
bool rightZ,
78 const Axis& tAxis,
bool leftT,
bool rightT,
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,
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;
118 inline const std::vector<Axis>&
axes()
const {
return axes_;}
119 inline const Axis&
axis(
const unsigned i)
const
120 {
return axes_.at(i);}
136 void getCoords(
unsigned long linearIndex,
137 double* coords,
unsigned coordsBufferSize)
const;
161 template <
typename ConvertibleToUn
signed>
163 ConvertibleToUnsigned axisNumber,
const Axis& replacementAxis,
164 bool newAxisLeftLinear,
bool newAxisRightLinear,
189 template <
class Functor1,
class Functor2>
191 unsigned axisNumber,
const Axis& replacementAxis,
192 bool newAxisLeftLinear,
bool newAxisRightLinear,
201 {
return !(*
this ==
r);}
205 inline gs::ClassId
classId()
const {
return gs::ClassId(*
this);}
206 bool write(std::ostream& of)
const;
212 const gs::ClassId&
id, std::istream&
in);
219 const std::vector<Axis>&
axes,
220 const char* leftInterpolation,
221 const char* rightInterpolation,
234 template <
class Functor1>
236 double rmin,
double rmax,
239 template <
class Functor1>
241 const Axis& fromAxis,
const Axis& toAxis,
242 bool newLeftLinear,
bool newRightLinear,
244 const double* rawx,
const double* rawf,
254 #include <functional>
256 #include "Alignment/Geners/interface/binaryIO.hh"
263 template <
class Axis>
266 const unsigned n = axes.size();
269 for (
unsigned i=0;
i<
n; ++
i)
270 result.push_back(axes[
i].nCoords());
274 template <
class Axis>
279 result.push_back(xAxis.nCoords());
283 template <
class Axis>
288 result.push_back(xAxis.nCoords());
289 result.push_back(yAxis.nCoords());
293 template <
class Axis>
300 result.push_back(xAxis.nCoords());
301 result.push_back(yAxis.nCoords());
302 result.push_back(zAxis.nCoords());
306 template <
class Axis>
308 const Axis& zAxis,
const Axis& tAxis)
312 result.push_back(xAxis.nCoords());
313 result.push_back(yAxis.nCoords());
314 result.push_back(zAxis.nCoords());
315 result.push_back(tAxis.nCoords());
319 template <
class Axis>
321 const Axis& zAxis,
const Axis& tAxis,
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());
335 const double x0,
const double x1,
336 const double y0,
const double y1,
339 return y0 + (y1 - y0)*((x - x0)/(x1 - x0));
342 template <
typename Numeric,
class Axis>
345 const Axis& fromAxis,
const Axis& toAxis,
346 const bool leftLinear,
const bool rightLinear,
350 assert(fromSlice.
rank() == 1U);
351 assert(toSlice->
rank() == 1U);
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;
360 "In npstat::Private::lind_invert1DSlice: "
361 "slice data is not monotonous and can not be inverted");
363 const Numeric yfirst = fromData[0];
364 const Numeric ylast = fromData[fromLen - 1U];
365 const bool increasing = yfirst < ylast;
367 Numeric* toD =
const_cast<Numeric*
>(toSlice->
data());
368 const unsigned nAxisPoints = toAxis.nCoords();
369 assert(toSlice->
length() == nAxisPoints);
371 for (
unsigned ipt=0; ipt<nAxisPoints; ++ipt)
373 const Numeric
y =
static_cast<Numeric
>(toAxis.coordinate(ipt));
380 yfirst, fromData[1], fromAxis.coordinate(0),
381 fromAxis.coordinate(1),
y);
383 toD[ipt] = fromAxis.coordinate(0);
389 ylast, fromData[fromLen - 2U],
390 fromAxis.coordinate(fromLen - 1U),
391 fromAxis.coordinate(fromLen - 2U),
y);
393 toD[ipt] = fromAxis.coordinate(fromLen - 1U);
397 const unsigned i = std::lower_bound(fromData,fromDataEnd,y)-
400 fromData[i-1U], fromData[i],
401 fromAxis.coordinate(i-1U),
402 fromAxis.coordinate(i),
y);
413 ylast, fromData[fromLen - 2U],
414 fromAxis.coordinate(fromLen - 1U),
415 fromAxis.coordinate(fromLen - 2U),
y);
417 toD[ipt] = fromAxis.coordinate(fromLen - 1U);
419 else if (y >= yfirst)
424 fromAxis.coordinate(0),
425 fromAxis.coordinate(1),
y);
427 toD[ipt] = fromAxis.coordinate(0);
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);
443 template <
class Numeric,
class Axis>
446 for (
unsigned i=0;
i<dim_; ++
i)
447 if (leftInterpolationLinear_[
i] || rightInterpolationLinear_[
i])
452 template <
class Numeric,
class Axis>
458 for (
unsigned i=0;
i<dim_; ++
i)
465 return data_ == r.
data_ &&
470 template <
typename Numeric,
class Axis>
473 static const std::string myClass(gs::template_class_name<Numeric,Axis>(
474 "npstat::LinInterpolatedTableND"));
475 return myClass.c_str();
478 template<
typename Numeric,
class Axis>
481 const bool status = data_.classId().write(of) &&
483 gs::write_obj_vector(of, axes_);
486 gs::write_pod_array(of, leftInterpolationLinear_, dim_);
487 gs::write_pod_array(of, rightInterpolationLinear_, dim_);
488 gs::write_pod(of, functionLabel_);
490 return status && !of.fail();
493 template<
typename Numeric,
class Axis>
496 const gs::ClassId&
id, std::istream&
in)
498 static const gs::ClassId
current(
500 current.ensureSameId(
id);
502 gs::ClassId ida(in, 1);
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");
521 data, axes, leftInterpolation, rightInterpolation, label);
524 template<
typename Numeric,
class Axis>
526 const unsigned i)
const
529 "In npstat::LinInterpolatedTableND::leftInterpolationLinear: "
530 "index out of range");
531 return leftInterpolationLinear_[
i];
534 template<
typename Numeric,
class Axis>
536 const unsigned i)
const
539 "In npstat::LinInterpolatedTableND::rightInterpolationLinear: "
540 "index out of range");
541 return rightInterpolationLinear_[
i];
544 template<
typename Numeric,
class Axis>
547 for (
unsigned i=0;
i<dim_; ++
i)
548 if (!axes_[
i].isUniform())
553 template<
typename Numeric,
class Axis>
554 std::vector<std::pair<bool,bool> >
557 std::vector<std::pair<bool,bool> > vec;
559 for (
unsigned i=0;
i<dim_; ++
i)
560 vec.push_back(std::pair<bool, bool>(leftInterpolationLinear_[
i],
561 rightInterpolationLinear_[i]));
565 template<
typename Numeric,
class Axis>
567 const std::vector<Axis>& axes,
568 const std::vector<std::pair<bool,bool> >& interpolationType,
572 functionLabel_(label ? label :
""),
575 if (
dim_ == 0 ||
dim_ >= CHAR_BIT*
sizeof(
unsigned long))
577 "In npstat::LinInterpolatedTableND constructor: requested "
578 "table dimensionality is not supported");
581 "In npstat::LinInterpolatedTableND constructor: "
582 "incompatible number of interpolation specifications");
583 for (
unsigned i=0;
i<
dim_; ++
i)
593 template<
typename Numeric,
class Axis>
594 template <
class Num2>
599 functionLabel_(r.functionLabel_),
601 allConstInterpolated_(r.allConstInterpolated_)
603 for (
unsigned i=0;
i<
dim_; ++
i)
610 template<
typename Numeric,
class Axis>
613 const std::vector<Axis>& axes,
614 const char* leftInterpolation,
615 const char* rightInterpolation,
619 functionLabel_(label),
622 for (
unsigned i=0;
i<
dim_; ++
i)
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,
638 : data_(Private::
makeTableShape(xAxis, yAxis, zAxis, tAxis, vAxis)),
639 functionLabel_(label ? label :
""),
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);
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,
673 functionLabel_(label ? label :
""),
677 axes_.push_back(xAxis);
678 axes_.push_back(yAxis);
679 axes_.push_back(zAxis);
680 axes_.push_back(tAxis);
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,
703 functionLabel_(label ? label :
""),
707 axes_.push_back(xAxis);
708 axes_.push_back(yAxis);
709 axes_.push_back(zAxis);
723 template<
typename Numeric,
class Axis>
725 const Axis& xAxis,
bool leftX,
bool rightX,
726 const Axis& yAxis,
bool leftY,
bool rightY,
729 functionLabel_(label ? label :
""),
733 axes_.push_back(xAxis);
734 axes_.push_back(yAxis);
746 template<
typename Numeric,
class Axis>
748 const Axis& xAxis,
bool leftX,
bool rightX,
751 functionLabel_(label ? label :
""),
755 axes_.push_back(xAxis);
763 template<
typename Numeric,
class Axis>
764 template <
typename ConvertibleToUn
signed>
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
771 const unsigned axisNumber =
static_cast<unsigned>(axisNumC);
773 if (axisNumber >= dim_)
775 "In npstat::LinInterpolatedTableND::invertAxis: "
776 "axis number is out of range");
779 std::vector<Axis> newAxes(axes_);
780 newAxes[axisNumber] = replacementAxis;
782 std::vector<std::pair<bool,bool> > iType(interpolationType());
783 iType[axisNumber] = std::pair<bool,bool>(leftLinear, rightLinear);
786 CPP11_auto_ptr<LinInterpolatedTableND> pTable(
792 unsigned sliceIndex[CHAR_BIT*
sizeof(
unsigned long)];
793 unsigned fixedIndices[CHAR_BIT*
sizeof(
unsigned long)];
795 for (
unsigned i=0;
i<dim_; ++
i)
798 sliceIndex[
count] = data_.span(
i);
799 fixedIndices[count++] =
i;
807 scan.getIndex(sliceIndex, count);
808 data_.exportSlice(&parentSlice, fixedIndices,
811 parentSlice, axes_[axisNumber], replacementAxis,
812 leftLinear, rightLinear, &dauSlice);
813 pTable->data_.importSlice(dauSlice, fixedIndices,
819 data_, axes_[0], replacementAxis,
820 leftLinear, rightLinear, &pTable->data_);
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,
831 const char* functionLabel)
const
833 if (axisNumber >= dim_)
835 "In npstat::LinInterpolatedTableND::invertRatioResponse: "
836 "axis number is out of range");
839 std::vector<Axis> newAxes(axes_);
840 newAxes[axisNumber] = replacementAxis;
842 std::vector<std::pair<bool,bool> > iType(interpolationType());
843 iType[axisNumber] = std::pair<bool,bool>(leftLinear, rightLinear);
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)
852 const double x = invg(oldAxis.coordinate(
i));
855 "In npstat::LinInterpolatedTableND::invertRatioResponse: "
856 "invalid original axis definition (negative transformed "
862 std::vector<double> rawf;
863 const unsigned nFuncs = replacementAxis.nCoords();
864 rawf.reserve(nFuncs);
865 for (
unsigned i=0;
i<nFuncs; ++
i)
867 const double f = invh(replacementAxis.coordinate(
i));
870 "In npstat::LinInterpolatedTableND::invertRatioResponse: "
871 "invalid new axis definition (negative transformed "
877 std::vector<double> workspace(nCoords);
880 CPP11_auto_ptr<LinInterpolatedTableND> pTable(
886 unsigned sliceIndex[CHAR_BIT*
sizeof(
unsigned long)];
887 unsigned fixedIndices[CHAR_BIT*
sizeof(
unsigned long)];
889 for (
unsigned i=0;
i<dim_; ++
i)
892 sliceIndex[
count] = data_.span(
i);
893 fixedIndices[count++] =
i;
901 scan.getIndex(sliceIndex, count);
902 data_.exportSlice(&parentSlice, fixedIndices,
904 invert1DResponse(parentSlice, oldAxis,
905 replacementAxis, leftLinear, rightLinear,
906 invg, &rawx[0], &rawf[0], &workspace[0],
908 pTable->data_.importSlice(dauSlice, fixedIndices,
913 invert1DResponse(data_, oldAxis, replacementAxis, leftLinear,
914 rightLinear, invg, &rawx[0], &rawf[0],
915 &workspace[0], &pTable->data_);
919 template<
typename Numeric,
class Axis>
921 const unsigned long linearIndex,
922 double* coords,
const unsigned coordsBufferSize)
const
925 "In LinInterpolatedTableND::getCoords: "
926 "insufficient buffer size");
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]);
934 template<
typename Numeric,
class Axis>
936 const double*
point,
const unsigned len)
const
940 "In npstat::LinInterpolatedTableND::isWithinLimits: "
941 "incompatible point dimensionality");
944 for (
unsigned i=0;
i<dim_; ++
i)
945 if (point[
i] < axes_[
i].
min() || point[
i] > axes_[
i].
max())
950 template<
typename Numeric,
class Axis>
952 const double*
point,
const unsigned len)
const
958 "In npstat::LinInterpolatedTableND::operator(): "
959 "incompatible point dimensionality");
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()))
968 interpolateArray =
false;
972 if (interpolateArray)
976 double buf[CHAR_BIT*
sizeof(
unsigned long)];
977 for (
unsigned i=0;
i<dim_; ++
i)
979 const std::pair<unsigned,double>& pair =
980 axes_[
i].getInterval(point[
i]);
981 buf[
i] = pair.first + 1U - pair.second;
983 return data_.interpolate1(buf, dim_);
987 unsigned ix[CHAR_BIT*
sizeof(
unsigned long)];
988 double weight[CHAR_BIT*
sizeof(
unsigned long)];
989 for (
unsigned i=0;
i<dim_; ++
i)
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]);
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)
1009 unsigned long icell = 0UL;
1010 for (
unsigned i=0;
i<dim_; ++
i)
1012 if (icycle & (1UL <<
i))
1015 icell += strides[
i]*(ix[
i] + 1U);
1020 icell += strides[
i]*ix[
i];
1023 sum += dat[icell]*
static_cast<proper_double
>(
w);
1029 template<
typename Numeric,
class Axis>
1031 const double& x0)
const
1033 const unsigned nArgs = 1U;
1035 "In npstat::LinInterpolatedTableND::operator(): number of "
1036 "arguments, 1, is incompatible with the interpolator dimensionality");
1039 return operator()(tmp, nArgs);
1042 template<
typename Numeric,
class Axis>
1044 const double& x0,
const double& x1)
const
1046 const unsigned nArgs = 2U;
1048 "In npstat::LinInterpolatedTableND::operator(): number of "
1049 "arguments, 2, is incompatible with the interpolator dimensionality");
1053 return operator()(tmp, nArgs);
1056 template<
typename Numeric,
class Axis>
1058 const double& x0,
const double& x1,
const double& x2)
const
1060 const unsigned nArgs = 3U;
1062 "In npstat::LinInterpolatedTableND::operator(): number of "
1063 "arguments, 3, is incompatible with the interpolator dimensionality");
1068 return operator()(tmp, nArgs);
1071 template<
typename Numeric,
class Axis>
1073 const double& x0,
const double& x1,
1074 const double& x2,
const double& x3)
const
1076 const unsigned nArgs = 4U;
1078 "In npstat::LinInterpolatedTableND::operator(): number of "
1079 "arguments, 4, is incompatible with the interpolator dimensionality");
1085 return operator()(tmp, nArgs);
1088 template<
typename Numeric,
class Axis>
1090 const double& x0,
const double& x1,
const double& x2,
1091 const double& x3,
const double& x4)
const
1093 const unsigned nArgs = 5U;
1095 "In npstat::LinInterpolatedTableND::operator(): number of "
1096 "arguments, 5, is incompatible with the interpolator dimensionality");
1103 return operator()(tmp, nArgs);
1106 template<
typename Numeric,
class Axis>
1107 template <
class Functor1>
1109 const double xmin,
const double xmax,
1110 const double rmin,
const double rmax,
1116 double fmin = invg(xmin)*rmin;
1117 double fmax = invg(xmax)*rmax;
1119 assert(fmin < fmax);
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)
1130 x0 -= stepfactor*
step;
1132 xmin, xmax, rmin, rmax, x0);
1138 x1 += stepfactor*
step;
1140 xmin, xmax, rmin, rmax, x1);
1143 "In LinInterpolatedTableND::solveForRatioArg: "
1144 "faled to bracket the root");
1147 while ((x1 - x0)/(
std::abs(x1) +
std::abs(x0) + DBL_EPSILON) > 4.0*DBL_EPSILON)
1149 const double xhalf = (x1 + x0)/2.0;
1151 xmin, xmax, rmin, rmax, xhalf);
1163 return (x1 + x0)/2.0;
1166 template<
typename Numeric,
class Axis>
1167 template <
class Functor1>
1170 const Axis& fromAxis,
const Axis& toAxis,
1171 const bool newLeftLinear,
const bool newRightLinear,
1172 Functor1 invg,
const double* rawx,
const double* rawf,
1177 assert(fromSlice.
rank() == 1U);
1178 assert(toSlice->
rank() == 1U);
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);
1189 for (
unsigned i=0;
i<fromLen; ++
i)
1192 "In LinInterpolatedTableND::invert1DResponse: "
1193 "non-positive response found. This ratio "
1194 "response table is not invertible.");
1195 workspace[
i] = rawx[
i]*fromData[
i];
1198 const double yfirst = workspace[0];
1199 const double ylast = workspace[fromLen - 1U];
1201 bool adjustZero =
false;
1202 unsigned nBelow = 0;
1203 for (
unsigned ipt=0; ipt<nAxisPoints; ++ipt)
1205 const double y = rawf[ipt];
1214 else if (y <= yfirst)
1217 solve = newLeftLinear;
1219 else if (y >= ylast)
1221 solve = newRightLinear;
1222 i0 = solve ? fromLen-2 : fromLen-1;
1227 i0 =
static_cast<unsigned>(std::lower_bound(
1228 workspace,workspace+fromLen,y) - workspace) - 1U;
1232 const double x = solveForRatioArg(fromAxis.coordinate(i0),
1233 fromAxis.coordinate(i0+1),
1234 fromData[i0], fromData[i0+1],
1236 toD[ipt] = invg(x)/
y;
1239 toD[ipt] = 1.0/fromData[i0];
1241 if (adjustZero && nBelow)
1247 #endif // NPSTAT_LININTERPOLATEDTABLEND_HH_
bool rightInterpolationLinear(unsigned i) const
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)
CPP11_auto_ptr< LinInterpolatedTableND > invertRatioResponse(unsigned axisNumber, const Axis &replacementAxis, bool newAxisLeftLinear, bool newAxisRightLinear, Functor1 invg, Functor2 invh, const char *functionLabel=0) const
bool allConstInterpolated_
void getCoords(unsigned long linearIndex, double *coords, unsigned coordsBufferSize) const
std::vector< unsigned > ArrayShape
unsigned long length() const
std::string functionLabel_
const ArrayND< Numeric > & table() const
static double solveForRatioArg(double xmin, double xmax, double rmin, double rmax, double fval, Functor1 invg)
Numeric operator()(const double *point, unsigned dim) const
unsigned long length() const
bool write(std::ostream &of) const
A few simple template functions for checking monotonicity of container values.
const T & max(const T &a, const T &b)
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
std::vector< Axis > axes_
Abs< T >::type abs(const T &t)
const Numeric * data() const
void setFunctionLabel(const char *newlabel)
char rightInterpolationLinear_[CHAR_BIT *sizeof(unsigned long)]
static unsigned version()
bool isUniformlyBinned() 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)
ArrayND< Numeric > & table()
bool isStrictlyMonotonous(Iter const begin, Iter const end)
static void restore(const gs::ClassId &id, std::istream &in, ArrayND *array)
bool allConstInterpolated() const
CPP11_auto_ptr< LinInterpolatedTableND > invertWRTAxis(ConvertibleToUnsigned axisNumber, const Axis &replacementAxis, bool newAxisLeftLinear, bool newAxisRightLinear, const char *functionLabel=0) const
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
bool isWithinLimits(const double *point, unsigned dim) const
char data[epos_bytes_allocation]
const Axis & axis(const unsigned i) const
const std::string & functionLabel() const
Iteration over indices of a multidimensional array.
static const char * classname()
tuple size
Write out results.
gs::ClassId classId() const
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