Isis 3 Programmer Reference
NumericalApproximation.cpp
1 #include "NumericalApproximation.h"
2 
3 #include <cmath>
4 #include <iostream>
5 #include <sstream>
6 #include <string>
7 #include <vector>
8 #include <algorithm>
9 
10 #include "SpecialPixel.h"
11 #include "IException.h"
12 #include "IString.h"
13 
14 using namespace std;
15 namespace Isis {
16  NumericalApproximation::FunctorList NumericalApproximation::p_interpFunctors;
17 
35  NumericalApproximation::NumericalApproximation(const NumericalApproximation::InterpType &itype) {
36  // Validate parameter and initialize class variables
37  try {
38  Init(itype);
39  p_dataValidated = false;
40  }
41  catch(IException &e) { // catch exception from Init()
42  throw IException(e, e.errorType(),
43  "NumericalApproximation() - Unable to construct NumericalApproximation object",
44  _FILEINFO_);
45  }
46  }
47 
81  NumericalApproximation::NumericalApproximation(unsigned int n, double *x, double *y,
83  try {
84  Init(itype);
85  AddData(n, x, y);
86  ValidateDataSet();
87  }
88  catch(IException &e) { // catch exception from Init(), AddData(), ValidateDataSet()
89  throw IException(e, e.errorType(),
90  "NumericalApproximation() - Unable to construct object using the given arrays, size and interpolation type",
91  _FILEINFO_);
92  }
93  }
94 
125  NumericalApproximation::NumericalApproximation(const vector <double> &x, const vector <double> &y,
126  const NumericalApproximation::InterpType &itype) {
127  try {
128  Init(itype);
129  AddData(x, y); // size of x = size of y validated in this method
130  ValidateDataSet();
131  }
132  catch(IException &e) { // catch exception from Init(), AddData(), ValidateDataSet()
133  throw IException(e, e.errorType(),
134  "NumericalApproximation() - Unable to construct an object using the given vectors and interpolation type",
135  _FILEINFO_);
136  }
137  }
138 
164  NumericalApproximation::NumericalApproximation(const NumericalApproximation &oldObject) {
165  try {
166  // initialize new object and set type to that of old object
167  Init(oldObject.p_itype);
168  // fill data set of new object
169  p_x = oldObject.p_x;
170  p_y = oldObject.p_y;
171  // if this data was previously validated, no need to do this again
172  p_dataValidated = oldObject.p_dataValidated;
173  // copy values for interpolation-specific variables
174  p_clampedComputed = oldObject.p_clampedComputed;
175  p_clampedEndptsSet = oldObject.p_clampedEndptsSet;
176  p_clampedSecondDerivs = oldObject.p_clampedSecondDerivs;
177  p_clampedDerivFirstPt = oldObject.p_clampedDerivFirstPt;
178  p_clampedDerivLastPt = oldObject.p_clampedDerivLastPt;
179  p_polyNevError = oldObject.p_polyNevError;
180  p_fprimeOfx = oldObject.p_fprimeOfx;
181  }
182  catch(IException &e) { // catch exception from Init()
183  throw IException(e,
184  e.errorType(),
185  "NumericalApproximation() - Unable to copy the given object",
186  _FILEINFO_);
187  }
188  }
189 
214  NumericalApproximation &NumericalApproximation::operator=(const NumericalApproximation &oldObject) {
215  try {
216  if(&oldObject != this) {
217  if(GslInterpType(p_itype)) GslDeallocation();
218  SetInterpType(oldObject.p_itype);
219  // set class variables
220  p_x = oldObject.p_x;
221  p_y = oldObject.p_y;
222  p_dataValidated = oldObject.p_dataValidated;
223  p_clampedSecondDerivs = oldObject.p_clampedSecondDerivs;
224  p_clampedDerivFirstPt = oldObject.p_clampedDerivFirstPt;
225  p_clampedDerivLastPt = oldObject.p_clampedDerivLastPt;
226  p_polyNevError = oldObject.p_polyNevError;
227  p_fprimeOfx = oldObject.p_fprimeOfx;
228  }
229  return (*this);
230  }
231  catch(IException &e) { // catch exception from SetInterpType()
232  throw IException(e,
233  e.errorType(),
234  "operator=() - Unable to copy the given object",
235  _FILEINFO_);
236  }
237 
238  }
239 
250  NumericalApproximation::~NumericalApproximation() {
251  if(GslInterpType(p_itype)) GslDeallocation();
252  }
253 
276  string NumericalApproximation::Name() const {
277  if(p_itype == NumericalApproximation::CubicNeighborhood) {
278  return "cspline-neighborhood";
279  }
280  else if(p_itype == NumericalApproximation::CubicClamped) {
281  return "cspline-clamped";
282  }
283  else if(p_itype == NumericalApproximation::PolynomialNeville) {
284  return "polynomial-Neville's";
285  }
286  else if(p_itype == NumericalApproximation::CubicHermite) {
287  return "cspline-Hermite";
288  }
289  try {
290  string name = (string(GslFunctor(p_itype)->name));
291  if(p_itype == NumericalApproximation::CubicNatural) {
292  return name + "-natural";
293  }
294  else return name;
295  }
296  catch(IException &e) { // catch exception from GslFunctor()
297  throw IException(e,
298  e.errorType(),
299  "Name() - GSL interpolation type not found",
300  _FILEINFO_);
301  }
302  }
303 
327  int NumericalApproximation::MinPoints() {
328  if(p_itype == NumericalApproximation::CubicNeighborhood) {
329  return 4;
330  }
331  else if(p_itype == NumericalApproximation::CubicClamped) {
332  return 3;
333  }
334  else if(p_itype == NumericalApproximation::PolynomialNeville) {
335  return 3;
336  }
337  else if(p_itype == NumericalApproximation::CubicHermite) {
338  return 2;
339  }
340  try {
341  return (GslFunctor(p_itype)->min_size);
342  }
343  catch(IException &e) { // catch exception from GslFunctor()
344  throw IException(e,
345  e.errorType(),
346  "MinPoints() - GSL interpolation not found",
347  _FILEINFO_);
348  }
349  }
350 
378  int NumericalApproximation::MinPoints(NumericalApproximation::InterpType itype) {
379  //Validates the parameter
380  if(GslInterpType(itype)) {
381  try {
382  //GslFunctor() Validates parameter for Gsl interp types
383  NumericalApproximation nam(itype);
384  return (nam.GslFunctor(itype)->min_size);
385  }
386  catch(IException &e) { // catch exception from GslFunctor()
387  throw IException(e,
388  e.errorType(),
389  "MinPoints() - GSL interpolation type not found",
390  _FILEINFO_);
391  }
392  }
393  else if(itype == NumericalApproximation::CubicNeighborhood) {
394  return 4;
395  }
396  else if(itype == NumericalApproximation::CubicClamped) {
397  return 3;
398  }
399  else if(itype == NumericalApproximation::PolynomialNeville) {
400  return 3;
401  }
402  else if(itype == NumericalApproximation::CubicHermite) {
403  return 2;
404  }
405  throw IException(IException::Programmer,
406  "MinPoints() - Invalid argument. Unknown interpolation type: "
408  _FILEINFO_);
409  }
410 
434  void NumericalApproximation::AddData(const double x, const double y) {
435  p_x.push_back(x);
436  p_y.push_back(y);
437  p_clampedComputed = false;
438  p_clampedEndptsSet = false;
439  p_dataValidated = false;
440  p_clampedSecondDerivs.clear();
441  p_clampedDerivFirstPt = 0;
442  p_clampedDerivLastPt = 0;
443  p_polyNevError.clear();
444  p_interp = 0;
445  p_acc = 0;
446  return;
447  }
448 
467  void NumericalApproximation::AddData(unsigned int n, double *x, double *y) {
468  for(unsigned int i = 0; i < n; i++) {
469  p_x.push_back(x[i]);
470  p_y.push_back(y[i]);
471  }
472  p_clampedComputed = false;
473  p_clampedEndptsSet = false;
474  p_dataValidated = false;
475  p_clampedSecondDerivs.clear();
476  p_clampedDerivFirstPt = 0;
477  p_clampedDerivLastPt = 0;
478  p_polyNevError.clear();
479  p_interp = 0;
480  p_acc = 0;
481  return;
482  }
483 
506  void NumericalApproximation::AddData(const vector <double> &x,
507  const vector <double> &y)
508  {
509  int n = x.size();
510  int m = y.size();
511  if(n != m) {
512  ReportException(IException::Programmer, "AddData()",
513  "Invalid arguments. The sizes of the input vectors do "
514  "not match",
515  _FILEINFO_);
516  }
517 
518  // Avoid push_back if at all possible. These calls were consuming 10% of
519  // cam2map's run time on a line scan camera.
520  if(!p_x.empty() || !p_y.empty()) {
521  for(int i = 0; i < n; i++) {
522  p_x.push_back(x[i]);
523  p_y.push_back(y[i]);
524  }
525  }
526  else {
527  p_x = x;
528  p_y = y;
529  }
530 
531  p_clampedComputed = false;
532  p_clampedEndptsSet = false;
533  p_dataValidated = false;
534  p_clampedSecondDerivs.clear();
535  p_clampedDerivFirstPt = 0;
536  p_clampedDerivLastPt = 0;
537  p_polyNevError.clear();
538  p_interp = 0;
539  p_acc = 0;
540  return;
541  }
542 
559  void NumericalApproximation::SetCubicClampedEndptDeriv(double yp1, double ypn) {
560  if(p_itype != NumericalApproximation::CubicClamped) {
561  ReportException(IException::Programmer, "SetCubicClampedEndptDeriv()",
562  "This method is only valid for cspline-clamped interpolation, may not be used for "
563  + Name() + " interpolation",
564  _FILEINFO_);
565  }
566  p_clampedDerivFirstPt = yp1;
567  p_clampedDerivLastPt = ypn;
568  p_clampedEndptsSet = true;
569  return;
570  }
571 
586  void NumericalApproximation::AddCubicHermiteDeriv(unsigned int n, double *fprimeOfx) {
587  if(p_itype != NumericalApproximation::CubicHermite) {
588  ReportException(IException::Programmer, "SetCubicHermiteDeriv()",
589  "This method is only valid for cspline-Hermite interpolation, may not be used for "
590  + Name() + " interpolation",
591  _FILEINFO_);
592  }
593  for(unsigned int i = 0; i < n; i++) {
594  p_fprimeOfx.push_back(fprimeOfx[i]);
595  }
596  return;
597  }
611  void NumericalApproximation::AddCubicHermiteDeriv(
612  const vector <double> &fprimeOfx) {
613  if(p_itype != NumericalApproximation::CubicHermite) {
614  ReportException(IException::Programmer, "SetCubicHermiteDeriv()",
615  "This method is only valid for cspline-Hermite interpolation, may not be used for "
616  + Name() + " interpolation",
617  _FILEINFO_);
618  }
619 
620  // Avoid push_back if at all possible.
621  if(!p_fprimeOfx.empty()) {
622  for(unsigned int i = 0; i < fprimeOfx.size(); i++) {
623  p_fprimeOfx.push_back(fprimeOfx[i]);
624  }
625  }
626  else {
627  p_fprimeOfx = fprimeOfx;
628  }
629 
630  return;
631  }
645  void NumericalApproximation::AddCubicHermiteDeriv(
646  double fprimeOfx) {
647  if(p_itype != NumericalApproximation::CubicHermite) {
648  ReportException(IException::Programmer, "SetCubicHermiteDeriv()",
649  "This method is only valid for cspline-Hermite interpolation, may not be used for "
650  + Name() + " interpolation",
651  _FILEINFO_);
652  }
653  p_fprimeOfx.push_back(fprimeOfx);
654  return;
655  }
656 
681  vector <double> NumericalApproximation::CubicClampedSecondDerivatives() {
682  try {
683  if(p_itype != NumericalApproximation::CubicClamped)
684  ReportException(IException::Programmer, "CubicClampedSecondDerivatives()",
685  "This method is only valid for cspline-clamped interpolation type may not be used for "
686  + Name() + " interpolation",
687  _FILEINFO_);
688  if(!p_clampedComputed) ComputeCubicClamped();
689  return p_clampedSecondDerivs;
690  }
691  catch(IException &e) { // catch exception from ComputeCubicClamped()
692  throw IException(e,
693  e.errorType(),
694  "CubicClampedSecondDerivatives() - Unable to compute clamped cubic spline interpolation",
695  _FILEINFO_);
696  }
697  }
698 
717  double NumericalApproximation::DomainMinimum() {
718  try {
719  if(GslInterpType(p_itype)) {
720  // if GSL interpolation type, we need to compute, if not already done
721  if(!GslComputed()) ComputeGsl();
722  return (p_interp->interp->xmin);
723  }
724  else {
725  if(!p_dataValidated) ValidateDataSet();
726  return *min_element(p_x.begin(), p_x.end());
727  }
728  }
729  catch(IException &e) { // catch exception from ComputeGsl() or ValidateDataSet()
730  throw IException(e,
731  e.errorType(),
732  "DomainMinimum() - Unable to calculate the domain minimum for the data set",
733  _FILEINFO_);
734  }
735  }
736 
754  double NumericalApproximation::DomainMaximum() {
755  try {
756  if(GslInterpType(p_itype)) {
757  // if GSL interpolation type, we need to compute, if not already done
758  if(!GslComputed()) ComputeGsl();
759  return (p_interp->interp->xmax);
760  }
761  else {
762  if(!p_dataValidated) ValidateDataSet();
763  return *max_element(p_x.begin(), p_x.end());
764  }
765  }
766  catch(IException &e) { // catch exception from ComputeGsl() or ValidateDataSet()
767  throw IException(e,
768  e.errorType(),
769  "DomainMaximum() - Unable to calculate the domain maximum for the data set",
770  _FILEINFO_);
771  }
772  }
773 
787  bool NumericalApproximation::Contains(double x) {
788  return binary_search(p_x.begin(), p_x.end(), x);
789  }
790 
830  double NumericalApproximation::Evaluate(const double a, const ExtrapType &etype) {
831  try {
832  // a is const, so we must set up temporary variable in case value needs to be changed
833  double a0;
834  if(InsideDomain(a)) // this will validate data set, if not already done
835  a0 = a;
836  else a0 = ValueToExtrapolate(a, etype);
837  // perform interpolation/extrapoltion
838  if(p_itype == NumericalApproximation::CubicNeighborhood) {
839  return EvaluateCubicNeighborhood(a0);
840  }
841  else if(p_itype == NumericalApproximation::PolynomialNeville) {
842  p_polyNevError.clear();
843  return EvaluatePolynomialNeville(a0);
844  }
845  else if(p_itype == NumericalApproximation::CubicClamped) {
846  // if cubic clamped, we need to compute, if not already done
847  if(!p_clampedComputed) ComputeCubicClamped();
848  return EvaluateCubicClamped(a0);
849  }
850  else if(p_itype == NumericalApproximation::CubicHermite) {
851  return EvaluateCubicHermite(a0);
852  }
853  // if GSL interpolation type we need to compute, if not already done
854  if(!GslComputed()) ComputeGsl();
855  double result;
856  GslIntegrityCheck(gsl_spline_eval_e(p_interp, a0, p_acc, &result), _FILEINFO_);
857  return (result);
858  }
859  catch(IException &e) { // catch exception from EvaluateCubicNeighborhood(), EvaluateCubicClamped(), EvaluatePolynomialNeville(), GslIntegrityCheck()
860  throw IException(e,
861  e.errorType(),
862  "Evaluate() - Unable to evaluate the function at the point a = "
863  + IString(a),
864  _FILEINFO_);
865  }
866  }
867 
906  vector <double> NumericalApproximation::Evaluate(const vector <double> &a, const ExtrapType &etype) {
907  try {
908  if(p_itype == NumericalApproximation::CubicNeighborhood) {
909  // cubic neighborhood has it's own method that will take entire vector
910  // this is faster than looping through values calling Evaluate(double)
911  // for each component of the passed in vector
912  return EvaluateCubicNeighborhood(a, etype);
913  }
914  vector <double> result(a.size());
915  if(p_itype == NumericalApproximation::PolynomialNeville) {
916  // cannot loop through values calling Evaluate(double)
917  // because it will clear the p_polyNevError vector each time
918  p_polyNevError.clear();
919  for(unsigned int i = 0; i < result.size(); i++) {
920  double a0;
921  if(InsideDomain(a[i]))
922  a0 = a[i];
923  else a0 = ValueToExtrapolate(a[i], etype);
924  result[i] = EvaluatePolynomialNeville(a0);
925  }
926  return result;
927  }
928  // cubic-clamped, cubic-Hermite and gsl types can be done by calling Evaluate(double)
929  // for each value of the passed in vector
930  for(unsigned int i = 0; i < result.size(); i++) {
931  result[i] = Evaluate(a[i], etype);
932  }
933  return result;
934  }
935  catch(IException &e) { // catch exception from EvaluateCubicNeighborhood(), EvaluateCubicClamped(), EvaluatePolynomialNeville(), GslIntegrityCheck()
936  throw IException(e,
937  e.errorType(),
938  "Evaluate() - Unable to evaluate the function at the given vector of points",
939  _FILEINFO_);
940  }
941  }
942 
964  vector <double> NumericalApproximation::PolynomialNevilleErrorEstimate() {
965  if(p_itype != NumericalApproximation::PolynomialNeville) {
966  ReportException(IException::Programmer, "PolynomialNevilleErrorEstimate()",
967  "This method is only valid for polynomial-Neville's, may not be used for "
968  + Name() + " interpolation",
969  _FILEINFO_);
970  }
971  if(p_polyNevError.empty()) {
972  ReportException(IException::Programmer, "PolynomialNevilleErrorEstimate()",
973  "Error not calculated. This method only valid after Evaluate() has been called",
974  _FILEINFO_);
975  }
976  return p_polyNevError;
977  }
978 
1014  double NumericalApproximation::GslFirstDerivative(const double a) {
1015  try { // we need to compute, if not already done
1016  if(!GslComputed()) ComputeGsl();
1017  }
1018  catch(IException &e) { // catch exception from ComputeGsl()
1019  throw IException(e,
1020  e.errorType(),
1021  "GslFirstDerivative() - Unable to compute the first derivative at a = "
1022  + IString(a) + " using the GSL interpolation",
1023  _FILEINFO_);
1024  }
1025  if(!InsideDomain(a)) {
1026  ReportException(IException::Programmer, "GslFirstDerivative()",
1027  "Invalid argument. Value entered, a = "
1028  + IString(a) + ", is outside of domain = ["
1029  + IString(DomainMinimum()) + ", "
1030  + IString(DomainMaximum()) + "]",
1031  _FILEINFO_);
1032  }
1033  if(!GslInterpType(p_itype)) {
1034  ReportException(IException::Programmer, "GslFirstDerivative()",
1035  "Method only valid for GSL interpolation types, may not be used for "
1036  + Name() + " interpolation",
1037  _FILEINFO_);
1038  }
1039  try {
1040  double value;
1041  GslIntegrityCheck(gsl_spline_eval_deriv_e(p_interp, a, p_acc, &value), _FILEINFO_);
1042  return (value);
1043  }
1044  catch(IException &e) { // catch exception from GslIntegrityCheck()
1045  throw IException(e,
1046  e.errorType(),
1047  "GslFirstDerivative() - Unable to compute the first derivative at a = "
1048  + IString(a) + ". GSL integrity check failed",
1049  _FILEINFO_);
1050  }
1051  }
1052 
1065  double NumericalApproximation::EvaluateCubicHermiteFirstDeriv(const double a) {
1066  if(p_itype != NumericalApproximation::CubicHermite) { //??? is this necessary??? create single derivative method with GSL?
1067  ReportException(IException::User, "EvaluateCubicHermiteFirstDeriv()",
1068  "This method is only valid for cspline-Hermite interpolation, may not be used for "
1069  + Name() + " interpolation",
1070  _FILEINFO_);
1071  }
1072  if(p_fprimeOfx.size() != Size()) {
1073  ReportException(IException::User, "EvaluateCubicHermiteFirstDeriv()",
1074  "Invalid arguments. The size of the first derivative vector does not match the number of (x,y) data points.",
1075  _FILEINFO_);
1076  }
1077  // find the interval in which "a" exists
1078  int lowerIndex = FindIntervalLowerIndex(a);
1079 
1080  // we know that "a" is within the domain since this is verified in
1081  // Evaluate() before this method is called, thus n <= Size()
1082  if(a == p_x[lowerIndex]) {
1083  return p_fprimeOfx[lowerIndex];
1084  }
1085  if(a == p_x[lowerIndex+1]) {
1086  return p_fprimeOfx[lowerIndex+1];
1087  }
1088 
1089  double x0, x1, y0, y1, m0, m1;
1090  // a is contained within the interval (x0,x1)
1091  x0 = p_x[lowerIndex];
1092  x1 = p_x[lowerIndex+1];
1093  // the corresponding known y-values for x0 and x1
1094  y0 = p_y[lowerIndex];
1095  y1 = p_y[lowerIndex+1];
1096  // the corresponding known tangents (slopes) at (x0,y0) and (x1,y1)
1097  m0 = p_fprimeOfx[lowerIndex];
1098  m1 = p_fprimeOfx[lowerIndex+1];
1099 
1100  double h, t;
1101  h = x1 - x0;
1102  t = (a - x0) / h;
1103  if(h != 0.) {
1104  return ((6 * t * t - 6 * t) * y0 + (3 * t * t - 4 * t + 1) * h * m0 + (-6 * t * t + 6 * t) * y1 + (3 * t * t - 2 * t) * h * m1) / h;
1105  }
1106  else {
1107  return 0; // Should never happen
1108  }
1109  }
1110 
1153  double NumericalApproximation::BackwardFirstDifference(const double a, const unsigned int n, const double h) {
1154  if(!InsideDomain(a)) {
1155  ReportException(IException::Programmer, "BackwardFirstDifference()",
1156  "Invalid argument. Value entered, a = "
1157  + IString(a) + ", is outside of domain = ["
1158  + IString(DomainMinimum()) + ", "
1159  + IString(DomainMaximum()) + "]",
1160  _FILEINFO_);
1161  }
1162  if(!InsideDomain(a - (n - 1)*h)) {
1163  ReportException(IException::Programmer, "BackwardFirstDifference()",
1164  "Formula steps outside of domain. For "
1165  + IString((int) n) + "-point backward difference, a-(n-1)h = "
1166  + IString(a - (n - 1)*h) + " is smaller than domain min = "
1167  + IString(DomainMinimum())
1168  + ". Try forward difference or use smaller value for h or n",
1169  _FILEINFO_);
1170  }
1171  if(!p_dataValidated) ValidateDataSet();
1172  vector <double> f;
1173  double xi;
1174  try {
1175  for(double i = 0; i < n; i++) {
1176  xi = a + h * (i - (n - 1));
1177  f.push_back(Evaluate(xi)); // allow ExtrapType = ThrowError (default)
1178  }
1179  }
1180  catch(IException &e) { // catch exception from Evaluate()
1181  throw IException(e,
1182  e.errorType(),
1183  "BackwardFirstDifference() - Unable to calculate backward first difference for (a, n, h) = ("
1184  + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
1185  _FILEINFO_);
1186  }
1187  switch(n) {
1188  case 2:
1189  return (-f[0] + f[1]) / h; //2pt backward
1190  case 3:
1191  return (3 * f[2] - 4 * f[1] + f[0]) / (2 * h); //3pt backward
1192  default:
1193  throw IException(IException::Programmer,
1194  "BackwardFirstDifference() - Invalid argument. There is no "
1195  + IString((int) n) + "-point backward difference formula in use",
1196  _FILEINFO_);
1197  }
1198  }
1199 
1200 
1242  double NumericalApproximation::ForwardFirstDifference(const double a, const unsigned int n, const double h) {
1243  if(!InsideDomain(a)) {
1244  ReportException(IException::Programmer, "ForwardFirstDifference()",
1245  "Invalid argument. Value entered, a = "
1246  + IString(a) + ", is outside of domain = ["
1247  + IString(DomainMinimum()) + ", "
1248  + IString(DomainMaximum()) + "]",
1249  _FILEINFO_);
1250  }
1251  if(!InsideDomain(a + (n - 1)*h)) {
1252  ReportException(IException::Programmer, "ForwardFirstDifference()",
1253  "Formula steps outside of domain. For "
1254  + IString((int) n) + "-point forward difference, a+(n-1)h = "
1255  + IString(a + (n - 1)*h) + " is greater than domain max = "
1256  + IString(DomainMaximum())
1257  + ". Try backward difference or use smaller value for h or n",
1258  _FILEINFO_);
1259  }
1260  if(!p_dataValidated) ValidateDataSet();
1261  vector <double> f;
1262  double xi;
1263  try {
1264  for(double i = 0; i < n; i++) {
1265  xi = a + h * i;
1266  f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
1267  }
1268  }
1269  catch(IException &e) { // catch exception from Evaluate()
1270  throw IException(e,
1271  e.errorType(),
1272  "ForwardFirstDifference() - Unable to calculate forward first difference for (a, n, h) = ("
1273  + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
1274  _FILEINFO_);
1275  }
1276  switch(n) {
1277  case 2:
1278  return (-f[0] + f[1]) / h; //2pt forward
1279  case 3:
1280  return (-3 * f[0] + 4 * f[1] - f[2]) / (2 * h); //3pt forward
1281  default:
1282  throw IException(IException::Programmer,
1283  "ForwardFirstDifference() - Invalid argument. There is no "
1284  + IString((int) n) + "-point forward difference formula in use",
1285  _FILEINFO_);
1286  }
1287  }
1288 
1289 
1332  double NumericalApproximation::CenterFirstDifference(const double a, const unsigned int n, const double h) {
1333  if(!InsideDomain(a)) {
1334  ReportException(IException::Programmer, "CenterFirstDifference()",
1335  "Invalid argument. Value entered, a = "
1336  + IString(a) + ", is outside of domain = ["
1337  + IString(DomainMinimum()) + ", "
1338  + IString(DomainMaximum()) + "]",
1339  _FILEINFO_);
1340  }
1341  if(!InsideDomain(a + (n - 1)*h) || !InsideDomain(a - (n - 1)*h)) {
1342  ReportException(IException::Programmer, "CenterFirstDifference()",
1343  "Formula steps outside of domain. For "
1344  + IString((int) n) + "-point center difference, a-(n-1)h = "
1345  + IString(a - (n - 1)*h) + " or a+(n-1)h = "
1346  + IString(a + (n - 1)*h) + " is out of domain = ["
1347  + IString(DomainMinimum()) + ", " + IString(DomainMaximum())
1348  + "]. Use smaller value for h or n",
1349  _FILEINFO_);
1350  }
1351  if(!p_dataValidated) ValidateDataSet();
1352  vector <double> f;
1353  double xi;
1354  try {
1355  for(double i = 0; i < n; i++) {
1356  xi = a + h * (i - (n - 1) / 2);
1357  f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
1358  }
1359  }
1360  catch(IException &e) { // catch exception from Evaluate()
1361  throw IException(e,
1362  e.errorType(),
1363  "CenterFirstDifference() - Unable to calculate center first difference for (a, n, h) = ("
1364  + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
1365  _FILEINFO_);
1366  }
1367  switch(n) {
1368  case 3:
1369  return (-f[0] + f[2]) / (2 * h); //3pt center
1370  case 5:
1371  return (f[0] - 8 * f[1] + 8 * f[3] - f[4]) / (12 * h); //5pt center
1372  default:
1373  throw IException(IException::Programmer,
1374  "CenterFirstDifference() - Invalid argument. There is no "
1375  + IString((int) n) + "-point center difference formula in use",
1376  _FILEINFO_);
1377  }
1378  }
1379 
1416  double NumericalApproximation::GslSecondDerivative(const double a) {
1417  try { // we need to compute, if not already done
1418  if(!GslComputed()) ComputeGsl();
1419  }
1420  catch(IException &e) { // catch exception from ComputeGsl()
1421  throw IException(e,
1422  e.errorType(),
1423  "GslSecondDerivative() - Unable to compute the second derivative at a = "
1424  + IString(a) + " using the GSL interpolation",
1425  _FILEINFO_);
1426  }
1427  if(!InsideDomain(a))
1428  ReportException(IException::Programmer, "GslSecondDerivative()",
1429  "Invalid argument. Value entered, a = "
1430  + IString(a) + ", is outside of domain = ["
1431  + IString(DomainMinimum()) + ", "
1432  + IString(DomainMaximum()) + "]",
1433  _FILEINFO_);
1434  if(!GslInterpType(p_itype))
1435  ReportException(IException::Programmer, "GslSecondDerivative()",
1436  "Method only valid for GSL interpolation types, may not be used for "
1437  + Name() + " interpolation",
1438  _FILEINFO_);
1439  try {
1440  // we need to compute, if not already done
1441  if(!GslComputed()) ComputeGsl();
1442  double value;
1443  GslIntegrityCheck(gsl_spline_eval_deriv2_e(p_interp, a, p_acc, &value), _FILEINFO_);
1444  return (value);
1445  }
1446  catch(IException &e) { // catch exception from GslIntegrityCheck()
1447  throw IException(e,
1448  e.errorType(),
1449  "GslSecondDerivative() - Unable to compute the second derivative at a = "
1450  + IString(a) + ". GSL integrity check failed",
1451  _FILEINFO_);
1452  }
1453  }
1454 
1455 
1468  double NumericalApproximation::EvaluateCubicHermiteSecDeriv(const double a) {
1469  if(p_itype != NumericalApproximation::CubicHermite) {
1470  ReportException(IException::User, "EvaluateCubicHermiteSecDeriv()",
1471  "This method is only valid for cspline-Hermite interpolation, may not be used for "
1472  + Name() + " interpolation",
1473  _FILEINFO_);
1474  }
1475  if(p_fprimeOfx.size() != Size()) {
1476  ReportException(IException::User, "EvaluateCubicHermiteSecDeriv()",
1477  "Invalid arguments. The size of the first derivative vector does not match the number of (x,y) data points.",
1478  _FILEINFO_);
1479  }
1480  // find the interval in which "a" exists
1481  int lowerIndex = FindIntervalLowerIndex(a);
1482  double x0, x1, y0, y1, m0, m1;
1483  // a is contained within the interval (x0,x1)
1484  x0 = p_x[lowerIndex];
1485  x1 = p_x[lowerIndex+1];
1486  // the corresponding known y-values for x0 and x1
1487  y0 = p_y[lowerIndex];
1488  y1 = p_y[lowerIndex+1];
1489  // the corresponding known tangents (slopes) at (x0,y0) and (x1,y1)
1490  m0 = p_fprimeOfx[lowerIndex];
1491  m1 = p_fprimeOfx[lowerIndex+1];
1492 
1493  double h, t;
1494  h = x1 - x0;
1495  t = (a - x0) / h;
1496  if(h != 0.) {
1497  return ((12 * t - 6) * y0 + (6 * t - 4) * h * m0 + (-12 * t + 6) * y1 + (6 * t - 2) * h * m1) / h;
1498  }
1499  else {
1500  return 0; // Should never happen
1501  }
1502  }
1503 
1540  double NumericalApproximation::BackwardSecondDifference(const double a, const unsigned int n, const double h) {
1541  if(!InsideDomain(a)) {
1542  ReportException(IException::Programmer, "BackwardSecondDifference()",
1543  "Invalid argument. Value entered, a = "
1544  + IString(a) + ", is outside of domain = ["
1545  + IString(DomainMinimum()) + ", "
1546  + IString(DomainMaximum()) + "]",
1547  _FILEINFO_);
1548  }
1549  if(!InsideDomain(a - (n - 1)*h)) {
1550  ReportException(IException::Programmer, "BackwardSecondDifference()",
1551  "Formula steps outside of domain. For "
1552  + IString((int) n) + "-point backward difference, a-(n-1)h = "
1553  + IString(a - (n - 1)*h) + " is smaller than domain min = "
1554  + IString(DomainMinimum())
1555  + ". Try forward difference or use smaller value for h or n",
1556  _FILEINFO_);
1557  }
1558  if(!p_dataValidated) ValidateDataSet();
1559  vector <double> f;
1560  double xi;
1561  try {
1562  for(double i = 0; i < n; i++) {
1563  xi = a + h * (i - (n - 1));
1564  f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
1565  }
1566  }
1567  catch(IException &e) { // catch exception from Evaluate()
1568  throw IException(e,
1569  e.errorType(),
1570  "BackwardSecondDifference() - Unable to calculate backward second difference for (a, n, h) = ("
1571  + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
1572  _FILEINFO_);
1573  }
1574  switch(n) {
1575  case 3:
1576  return (f[0] - 2 * f[1] + f[2]) / (h * h); //3pt backward
1577  default:
1578  throw IException(IException::Programmer,
1579  "BackwardSecondDifference() - Invalid argument. There is no "
1580  + IString((int) n) + "-point backward second difference formula in use",
1581  _FILEINFO_);
1582  }
1583  }
1584 
1585 
1622  double NumericalApproximation::ForwardSecondDifference(const double a, const unsigned int n, const double h) {
1623  if(!InsideDomain(a)) {
1624  ReportException(IException::Programmer, "ForwardSecondDifference()",
1625  "Invalid argument. Value entered, a = "
1626  + IString(a) + ", is outside of domain = ["
1627  + IString(DomainMinimum()) + ", "
1628  + IString(DomainMaximum()) + "]",
1629  _FILEINFO_);
1630  }
1631  if(!InsideDomain(a + (n - 1)*h)) {
1632  ReportException(IException::Programmer, "ForwardSecondDifference()",
1633  "Formula steps outside of domain. For "
1634  + IString((int) n) + "-point forward difference, a+(n-1)h = "
1635  + IString(a + (n - 1)*h) + " is greater than domain max = "
1636  + IString(DomainMaximum())
1637  + ". Try backward difference or use smaller value for h or n",
1638  _FILEINFO_);
1639  }
1640  if(!p_dataValidated) ValidateDataSet();
1641  vector <double> f;
1642  double xi;
1643  try {
1644  for(double i = 0; i < n; i++) {
1645  xi = a + h * i;
1646  f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
1647  }
1648  }
1649  catch(IException &e) { // catch exception from Evaluate()
1650  throw IException(e,
1651  e.errorType(),
1652  "ForwardSecondDifference() - Unable to calculate forward second difference for (a, n, h) = ("
1653  + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
1654  _FILEINFO_);
1655  }
1656  switch(n) {
1657  case 3:
1658  return (f[0] - 2 * f[1] + f[2]) / (h * h); //3pt forward
1659  default:
1660  throw IException(IException::Programmer,
1661  "ForwardSecondDifference() - Invalid argument. There is no "
1662  + IString((int) n) + "-point forward second difference formula in use",
1663  _FILEINFO_);
1664  }
1665  }
1666 
1667 
1710  double NumericalApproximation::CenterSecondDifference(const double a, const unsigned int n, const double h) {
1711  if(!InsideDomain(a)) {
1712  ReportException(IException::Programmer, "CenterSecondDifference()",
1713  "Invalid argument. Value entered, a = "
1714  + IString(a) + ", is outside of domain = ["
1715  + IString(DomainMinimum()) + ", "
1716  + IString(DomainMaximum()) + "]",
1717  _FILEINFO_);
1718  }
1719  if(!InsideDomain(a + (n - 1)*h) || !InsideDomain(a - (n - 1)*h)) {
1720  ReportException(IException::Programmer, "CenterSecondDifference()",
1721  "Formula steps outside of domain. For "
1722  + IString((int) n) + "-point center difference, a-(n-1)h = "
1723  + IString(a - (n - 1)*h) + " or a+(n-1)h = "
1724  + IString(a + (n - 1)*h) + " is out of domain = ["
1725  + IString(DomainMinimum()) + ", " + IString(DomainMaximum())
1726  + "]. Use smaller value for h or n",
1727  _FILEINFO_);
1728  }
1729  if(!p_dataValidated) ValidateDataSet();
1730  vector <double> f;
1731  double xi;
1732  try {
1733  for(double i = 0; i < n; i++) {
1734  xi = a + h * (i - (n - 1) / 2);
1735  f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
1736  }
1737  }
1738  catch(IException &e) { // catch exception from Evaluate()
1739  throw IException(e,
1740  e.errorType(),
1741  "CenterSecondDifference() - Unable to calculate center second difference for (a, n, h) = ("
1742  + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
1743  _FILEINFO_);
1744  }
1745  switch(n) {
1746  case 3:
1747  return (f[0] - 2 * f[1] + f[2]) / (h * h); //3pt center
1748  case 5:
1749  return (-f[0] + 16 * f[1] - 30 * f[2] + 16 * f[3] - f[4]) / (12 * h * h); //5pt center
1750  default:
1751  throw IException(IException::Programmer,
1752  "CenterSecondDifference() - Invalid argument. There is no "
1753  + IString((int) n) + "-point center second difference formula in use",
1754  _FILEINFO_);
1755  }
1756  }
1757 
1791  double NumericalApproximation::GslIntegral(const double a, const double b) {
1792  try { // we need to compute, if not already done
1793  if(!GslComputed()) ComputeGsl();
1794  }
1795  catch(IException &e) { // catch exception from ComputeGsl()
1796  throw IException(e,
1797  e.errorType(),
1798  "GslIntegral() - Unable to compute the integral on the interval (a,b) = ("
1799  + IString(a) + ", " + IString(b)
1800  + ") using the GSL interpolation",
1801  _FILEINFO_);
1802  }
1803  if(a > b) {
1804  ReportException(IException::Programmer, "GslIntegral()",
1805  "Invalid interval entered: [a,b] = ["
1806  + IString(a) + ", " + IString(b) + "]",
1807  _FILEINFO_);
1808  }
1809  if(!InsideDomain(a) || !InsideDomain(b)) {
1810  ReportException(IException::Programmer, "GslIntegral()",
1811  "Invalid arguments. Interval entered ["
1812  + IString(a) + ", " + IString(b)
1813  + "] is not contained within domain ["
1814  + IString(DomainMinimum()) + ", "
1815  + IString(DomainMaximum()) + "]",
1816  _FILEINFO_);
1817  }
1818  if(!GslInterpType(p_itype)) {
1819  ReportException(IException::Programmer, "GslIntegral()",
1820  "Method only valid for GSL interpolation types, may not be used for "
1821  + Name() + " interpolation",
1822  _FILEINFO_);
1823  }
1824  try {
1825  double value;
1826  GslIntegrityCheck(gsl_spline_eval_integ_e(p_interp, a, b, p_acc, &value), _FILEINFO_);
1827  return (value);
1828  }
1829  catch(IException &e) { // catch exception from GslIntegrityCheck()
1830  throw IException(e,
1831  e.errorType(),
1832  "GslIntegral() - Unable to compute the integral on the interval (a,b) = ("
1833  + IString(a) + ", " + IString(b)
1834  + "). GSL integrity check failed",
1835  _FILEINFO_);
1836  }
1837  }
1838 
1862  double NumericalApproximation::TrapezoidalRule(const double a, const double b) {
1863  try {
1864  //n is the number of points used in the formula of the integration type chosen
1865  unsigned int n = 2;
1866  double result = 0;
1867  vector <double> f = EvaluateForIntegration(a, b, n);
1868  double h = f.back();
1869  f.pop_back();
1870  //Compute the integral using composite trapezoid formula
1871  int ii;
1872  for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
1873  ii = (i + 1) * (n - 1);
1874  result += (f[ii-1] + f[ii]) * h / 2;
1875  }
1876  return result;
1877  }
1878  catch(IException &e) { // catch exception from EvaluateForIntegration()
1879  throw IException(e,
1880  e.errorType(),
1881  "TrapezoidalRule() - Unable to calculate the integral on the interval (a,b) = ("
1882  + IString(a) + ", " + IString(b)
1883  + ") using the trapeziodal rule",
1884  _FILEINFO_);
1885  }
1886  }
1887 
1914  double NumericalApproximation::Simpsons3PointRule(const double a, const double b) {
1915  try {
1916  //n is the number of points used in the formula of the integration type chosen
1917  unsigned int n = 3;
1918  double result = 0;
1919  vector <double> f = EvaluateForIntegration(a, b, n);
1920  double h = f.back();
1921  f.pop_back();
1922  int ii;
1923  for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
1924  ii = (i + 1) * (n - 1);
1925  result += (f[ii-2] + 4 * f[ii-1] + f[ii]) * h / 3;
1926  }
1927  return result;
1928  }
1929  catch(IException &e) { // catch exception from EvaluateForIntegration()
1930  throw IException(e,
1931  e.errorType(),
1932  "Simpsons3PointRule() - Unable to calculate the integral on the interval (a,b) = ("
1933  + IString(a) + ", " + IString(b)
1934  + ") using Simpson's 3 point rule",
1935  _FILEINFO_);
1936  }
1937  }
1938 
1939 
1966  double NumericalApproximation::Simpsons4PointRule(const double a, const double b) {
1967  try {
1968  //n is the number of points used in the formula of the integration type chosen
1969  unsigned int n = 4;
1970  double result = 0;
1971  vector <double> f = EvaluateForIntegration(a, b, n);
1972  double h = f.back();
1973  f.pop_back();
1974  int ii;
1975  for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
1976  ii = (i + 1) * (n - 1);
1977  result += (f[ii-3] + 3 * f[ii-2] + 3 * f[ii-1] + f[ii]) * h * 3 / 8;
1978  }
1979  return result;
1980  }
1981  catch(IException &e) { // catch exception from EvaluateForIntegration()
1982  throw IException(e,
1983  e.errorType(),
1984  "Simpsons4PointRule() - Unable to calculate the integral on the interval (a,b) = ("
1985  + IString(a) + ", " + IString(b)
1986  + ") using Simpson's 4 point rule",
1987  _FILEINFO_);
1988  }
1989  }
1990 
1991 
2021  double NumericalApproximation::BoolesRule(const double a, const double b) {
2022  try {
2023  //n is the number of points used in the formula of the integration type chosen
2024  unsigned int n = 5;
2025  double result = 0;
2026  vector <double> f = EvaluateForIntegration(a, b, n);
2027  double h = f.back();
2028  f.pop_back();
2029  int ii;
2030  for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
2031  ii = (i + 1) * (n - 1);
2032  result += (7 * f[ii-4] + 32 * f[ii-3] + 12 * f[ii-2] + 32 * f[ii-1] + 7 * f[ii]) * h * 2 / 45;
2033  }
2034  return result;
2035  }
2036  catch(IException &e) { // catch exception from EvaluateForIntegration()
2037  throw IException(e,
2038  e.errorType(),
2039  "BoolesRule() - Unable to calculate the integral on the interval (a,b) = ("
2040  + IString(a) + ", " + IString(b)
2041  + ") using Boole's rule",
2042  _FILEINFO_);
2043  }
2044  }
2045 
2090  double NumericalApproximation::RefineExtendedTrap(double a, double b, double s, unsigned int n) {
2091  // This method was derived from an algorithm in the text
2092  // Numerical Recipes in C: The Art of Scientific Computing
2093  // Section 4.2 by Flannery, Press, Teukolsky, and Vetterling
2094  try {
2095  if(n == 1) {
2096  double begin, end;
2097  if(GslInterpType(p_itype) || p_itype == NumericalApproximation::CubicNeighborhood) {
2098  // if a or b are outside the domain, return y-value of nearest endpoint
2099  begin = Evaluate(a, NumericalApproximation::NearestEndpoint);
2100  end = Evaluate(b, NumericalApproximation::NearestEndpoint);
2101  }
2102  else {
2103  // if a or b are outside the domain, return extrapolated y-value
2104  begin = Evaluate(a, NumericalApproximation::Extrapolate);
2105  end = Evaluate(b, NumericalApproximation::Extrapolate);
2106  }
2107  return (0.5 * (b - a) * (begin + end));
2108  }
2109  else {
2110  int it;
2111  double delta, tnm, x, sum;
2112  it = (int)(pow(2.0, (double)(n - 2)));
2113  tnm = it;
2114  delta = (b - a) / tnm; // spacing of the points to be added
2115  x = a + 0.5 * delta;
2116  sum = 0.0;
2117  for(int i = 0; i < it; i++) {
2118  if(GslInterpType(p_itype) || p_itype == NumericalApproximation::CubicNeighborhood) {
2119  sum = sum + Evaluate(x, NumericalApproximation::NearestEndpoint);
2120  }
2121  else {
2122  sum = sum + Evaluate(x, NumericalApproximation::Extrapolate);
2123  }
2124  x = x + delta;
2125  }
2126  return (0.5 * (s + (b - a) * sum / tnm));// return refined value of s
2127  }
2128  }
2129  catch(IException &e) { // catch exception from Evaluate()
2130  throw IException(e,
2131  e.errorType(),
2132  "RefineExtendedTrap() - Unable to calculate the integral on the interval (a,b) = ("
2133  + IString(a) + ", " + IString(b)
2134  + ") using the extended trapeziodal rule",
2135  _FILEINFO_);
2136  }
2137  }
2138 
2173  double NumericalApproximation::RombergsMethod(double a, double b) {
2174  // This method was derived from an algorithm in the text
2175  // Numerical Recipes in C: The Art of Scientific Computing
2176  // Section 4.3 by Flannery, Press, Teukolsky, and Vetterling
2177  int maxits = 20;
2178  double dss = 0; // error estimate for
2179  double h[maxits+1]; // relative stepsizes for trap
2180  double trap[maxits+1]; // successive trapeziodal approximations
2181  double epsilon; // desired fractional accuracy
2182  double epsilon2;// desired fractional accuracy
2183  double ss; // result
2184 
2185  epsilon = 1.0e-4;
2186  epsilon2 = 1.0e-6;
2187 
2188  h[0] = 1.0;
2189  try {
2190  NumericalApproximation interp(NumericalApproximation::PolynomialNeville);
2191  for(int i = 0; i < maxits; i++) {
2192  // i will determine number of trapezoidal partitions of area
2193  // under curve for "integration" using refined trapezoidal rule
2194  trap[i] = RefineExtendedTrap(a, b, trap[i], i + 1); // validates data here
2195  if(i >= 4) {
2196  interp.AddData(5, &h[i-4], &trap[i-4]);
2197  // PolynomialNeville can extrapolate data outside of domain
2198  ss = interp.Evaluate(0.0, NumericalApproximation::Extrapolate);
2199  dss = interp.PolynomialNevilleErrorEstimate()[0];
2200  interp.Reset();
2201  // we work only until our necessary accuracy is achieved
2202  if(fabs(dss) <= epsilon * fabs(ss)) return ss;
2203  if(fabs(dss) <= epsilon2) return ss;
2204  }
2205  trap[i+1] = trap[i];
2206  h[i+1] = 0.25 * h[i];
2207  // This is a key step: the factor is 0.25d0 even though
2208  // the stepsize is decreased by 0.5d0. This makes the
2209  // extraplolation a polynomial in h-squared as allowed
2210  // by the equation from Numerical Recipes 4.2.1 pg.132,
2211  // not just a polynomial in h.
2212  }
2213  }
2214  catch(IException &e) { // catch error from RefineExtendedTrap, Constructor, Evaluate, PolynomialNevilleErrorEstimate
2215  throw IException(e,
2216  e.errorType(),
2217  "RombergsMethod() - Unable to calculate the integral on the interval (a,b) = ("
2218  + IString(a) + ", " + IString(b)
2219  + ") using Romberg's method",
2220  _FILEINFO_);
2221  }
2222  throw IException(IException::Programmer,
2223  "RombergsMethod() - Unable to calculate the integral using RombergsMethod() - Failed to converge in "
2224  + IString(maxits) + " iterations",
2225  _FILEINFO_);
2226  }
2227 
2245  void NumericalApproximation::Reset() {
2246  if(GslInterpType(p_itype)) GslDeallocation();
2247  p_clampedComputed = false;
2248  p_clampedEndptsSet = false;
2249  p_dataValidated = false;
2250  p_x.clear();
2251  p_y.clear();
2252  p_clampedSecondDerivs.clear();
2253  p_clampedDerivFirstPt = 0;
2254  p_clampedDerivLastPt = 0;
2255  p_polyNevError.clear();
2256  p_fprimeOfx.clear();
2257  return;
2258  }
2259 
2277  void NumericalApproximation::Reset(NumericalApproximation::InterpType itype) {
2278  try {
2279  Reset();
2280  SetInterpType(itype);
2281  return;
2282  }
2283  catch(IException &e) { // catch exception from SetInterpType()
2284  throw IException(e,
2285  e.errorType(),
2286  "Reset() - Unable to reset interpolation type",
2287  _FILEINFO_);
2288  }
2289  }
2290 
2313  void NumericalApproximation::SetInterpType(NumericalApproximation::InterpType itype) {
2314  // Validates the parameter
2315  if(GslInterpType(itype)) {
2316  try {
2317  GslFunctor(itype);
2318  }
2319  catch(IException &e) { // catch exception from GslFunctor()
2320  throw IException(e,
2321  e.errorType(),
2322  "SetInterpType() - Unable to set interpolation type",
2323  _FILEINFO_);
2324  }
2325  }
2326  else if(itype > 9) { // there are currently 9 interpolation types
2327  ReportException(IException::Programmer, "SetInterpType()",
2328  "Invalid argument. Unknown interpolation type: "
2330  _FILEINFO_);
2331  }
2332  // p_x, p_y are kept and p_itype is replaced
2333  p_itype = itype;
2334  // reset state of class variables that are InterpType dependent //??? should we keep some of this info?
2335  p_dataValidated = false;
2336  p_clampedComputed = false;
2337  p_clampedEndptsSet = false;
2338  p_clampedSecondDerivs.clear();
2339  p_clampedDerivFirstPt = 0;
2340  p_clampedDerivLastPt = 0;
2341  p_polyNevError.clear();
2342  p_fprimeOfx.clear();
2343  }
2344 
2369  void NumericalApproximation::Init(NumericalApproximation::InterpType itype) {
2370  if(p_interpFunctors.empty()) {
2371  p_interpFunctors.insert(make_pair(Linear, gsl_interp_linear));
2372  p_interpFunctors.insert(make_pair(Polynomial, gsl_interp_polynomial));
2373  p_interpFunctors.insert(make_pair(CubicNatural, gsl_interp_cspline));
2374  p_interpFunctors.insert(make_pair(CubicNatPeriodic, gsl_interp_cspline_periodic));
2375  p_interpFunctors.insert(make_pair(Akima, gsl_interp_akima));
2376  p_interpFunctors.insert(make_pair(AkimaPeriodic, gsl_interp_akima_periodic));
2377  }
2378 
2379  p_acc = 0;
2380  p_interp = 0;
2381  try {
2382  SetInterpType(itype);
2383  }
2384  catch(IException &e) { // catch exception from SetInterpType()
2385  throw IException(e,
2386  e.errorType(),
2387  "Init() - Unable to initialize NumericalApproximation object",
2388  _FILEINFO_);
2389  }
2390  // Turn all GSL error handling off...repeatedly, every time this routine is
2391  // called.
2392  gsl_set_error_handler_off();
2393  }
2394 
2412  bool NumericalApproximation::GslInterpType(NumericalApproximation::InterpType itype) const {
2413  if(itype == NumericalApproximation::Linear) return true;
2414  if(itype == NumericalApproximation::Polynomial) return true;
2415  if(itype == NumericalApproximation::CubicNatural) return true;
2416  if(itype == NumericalApproximation::CubicNatPeriodic) return true;
2417  if(itype == NumericalApproximation::Akima) return true;
2418  if(itype == NumericalApproximation::AkimaPeriodic) return true;
2419  return false;
2420  }
2421 
2438  void NumericalApproximation::GslAllocation(unsigned int npoints) {
2439  GslDeallocation();
2440  //get pointer to accelerator object (iterator for interpolation lookups)
2441  p_acc = gsl_interp_accel_alloc();
2442  //get pointer to interpolation object of interp type given for npoints datapoints
2443  p_interp = gsl_spline_alloc(GslFunctor(p_itype), npoints);
2444  return;
2445  }
2446 
2460  void NumericalApproximation::GslDeallocation() {
2461  if(p_interp) gsl_spline_free(p_interp);
2462  if(p_acc) gsl_interp_accel_free(p_acc);
2463  p_acc = 0;
2464  p_interp = 0;
2465  return;
2466  }
2467 
2487  const {
2488  FunctorConstIter fItr = p_interpFunctors.find(itype);
2489  if(fItr == p_interpFunctors.end()) {
2490  ReportException(IException::Programmer, "GslFunctor()",
2491  "Invalid argument. Unable to find GSL interpolator with id = "
2493  _FILEINFO_);
2494  }
2495  return (fItr->second);
2496  }
2497 
2517  void NumericalApproximation::GslIntegrityCheck(int gsl_status, const char *src, int line)
2518  {
2519  if(gsl_status != GSL_SUCCESS) {
2520  if(gsl_status != GSL_EDOM) {
2521  ReportException(IException::Programmer, "GslIntegrityCheck(int,char,int)",
2522  "GslIntegrityCheck(): GSL error occured: "
2523  + string(gsl_strerror(gsl_status)), src, line);
2524  }
2525  }
2526  return;
2527  }
2528 
2557  void NumericalApproximation::ValidateDataSet() {
2558  if((int) Size() < MinPoints()) {
2559  ReportException(IException::Programmer, "ValidateDataSet()",
2560  Name() + " interpolation requires a minimum of "
2561  + IString(MinPoints()) + " data points - currently have "
2562  + IString((int) Size()),
2563  _FILEINFO_);
2564  }
2565  for(unsigned int i = 1; i < Size(); i++) {
2566  // Check for uniqueness -- this applies to all interpolation types
2567  if(p_x[i-1] == p_x[i]) {
2568  ReportException(IException::Programmer, "ValidateDataSet()",
2569  "Invalid data set, x-values must be unique: \n\t\tp_x["
2570  + IString((int) i - 1) + "] = " + IString(p_x[i-1])
2571  + " = p_x[" + IString((int) i) + "]",
2572  _FILEINFO_);
2573  }
2574  if(p_x[i-1] > p_x[i]) {
2575  // Verify that data set is in ascending order --
2576  // this does not apply to PolynomialNeville, which appears to get the same results with unsorted data
2577  if(p_itype != NumericalApproximation::PolynomialNeville) {
2578  ReportException(IException::Programmer, "ValidateDataSet()",
2579  "Invalid data set, x-values must be in ascending order for "
2580  + Name() + " interpolation: \n\t\tx["
2581  + IString((int) i - 1) + "] = " + IString(p_x[i-1]) + " > x["
2582  + IString((int) i) + "] = " + IString(p_x[i]),
2583  _FILEINFO_);
2584  }
2585  }
2586  }
2587  if(p_itype == NumericalApproximation::CubicNatPeriodic) {
2588  if(p_y[0] != p_y[Size()-1]) {
2589  ReportException(IException::Programmer, "ValidateDataSet()",
2590  "First and last points of the data set must have the same y-value for "
2591  + Name() + "interpolation to prevent discontinuity at the boundary",
2592  _FILEINFO_);
2593  }
2594  }
2595  p_dataValidated = true;
2596  return;
2597  }
2598 
2599 
2613  bool NumericalApproximation::InsideDomain(const double a) {
2614  try {
2615  if(a + DBL_EPSILON < DomainMinimum()) {
2616  return false;
2617  }
2618  if(a - DBL_EPSILON > DomainMaximum()) {
2619  return false;
2620  }
2621  }
2622  catch(IException &e) { // catch exception from DomainMinimum(), DomainMaximum()
2623  throw IException(e,
2624  e.errorType(),
2625  "InsideDomain() - Unable to compute domain boundaries",
2626  _FILEINFO_);
2627  }
2628  return true;
2629  }
2630 
2647  bool NumericalApproximation::GslComputed() const {
2648  if(GslInterpType(p_itype)) return ((p_interp) && (p_acc));
2649  else
2650  throw IException(IException::Programmer,
2651  "GslComputed() - Method only valid for GSL interpolation types, may not be used for "
2652  + Name() + " interpolation",
2653  _FILEINFO_);
2654  }
2655 
2677  void NumericalApproximation::ComputeGsl() {
2678  try {
2679  if(GslComputed()) return;
2680  if(!p_dataValidated) ValidateDataSet();
2681  GslAllocation(Size());
2682  GslIntegrityCheck(gsl_spline_init(p_interp, &p_x[0], &p_y[0], Size()), _FILEINFO_);
2683  return;
2684  }
2685  catch(IException &e) { // catch exception from ValidateDataSet(), GslIntegrityCheck()
2686  throw IException(e,
2687  e.errorType(),
2688  "ComputeGsl() - Unable to compute GSL interpolation",
2689  _FILEINFO_);
2690  }
2691  }
2692 
2729  void NumericalApproximation::ComputeCubicClamped() {
2730  // This method was derived from an algorithm in the text
2731  // Numerical Recipes in C: The Art of Scientific Computing
2732  // Section 3.3 by Flannery, Press, Teukolsky, and Vetterling
2733 
2734  if(!p_dataValidated) {
2735  try {
2736  ValidateDataSet();
2737  }
2738  catch(IException &e) { // catch exception from ValidateDataSet()
2739  throw IException(e,
2740  e.errorType(),
2741  "ComputeCubicClamped() - Unable to compute cubic clamped interpolation",
2742  _FILEINFO_);
2743  }
2744  }
2745  if(!p_clampedEndptsSet) {
2746  ReportException(IException::Programmer, "ComputeCubicClamped()",
2747  "Must set endpoint derivative values after adding data in order to compute cubic spline with clamped boundary conditions",
2748  _FILEINFO_);
2749  }
2750  int n = Size();
2751  p_clampedSecondDerivs.resize(n);
2752  double u[n];
2753  double p, sig, qn, un;
2754 
2755  if(p_clampedDerivFirstPt > 0.99e30) {
2756  p_clampedSecondDerivs[0] = 0.0;//natural boundary conditions are used if deriv of first value is greater than 10^30
2757  u[0] = 0.0;
2758  }
2759  else {
2760  p_clampedSecondDerivs[0] = -0.5;// clamped conditions are used
2761  u[0] = (3.0 / (p_x[1] - p_x[0])) * ((p_y[1] - p_y[0]) / (p_x[1] - p_x[0]) - p_clampedDerivFirstPt);
2762  }
2763  for(int i = 1; i < n - 1; i++) { // decomposition loop of the tridiagonal algorithm
2764  sig = (p_x[i] - p_x[i-1]) / (p_x[i+1] - p_x[i-1]);
2765  p = sig * p_clampedSecondDerivs[i-1] + 2.0;
2766  p_clampedSecondDerivs[i] = (sig - 1.0) / p;
2767  u[i] = (6.0 * ((p_y[i+1] - p_y[i]) / (p_x[i+1] - p_x[i]) - (p_y[i] - p_y[i-1]) /
2768  (p_x[i] - p_x[i-1])) / (p_x[i+1] - p_x[i-1]) - sig * u[i-1]) / p;
2769  }
2770  if(p_clampedDerivLastPt > 0.99e30) { // upper boundary is natural
2771  qn = 0.0;
2772  un = 0.0;
2773  }
2774  else {// upper boundary is clamped
2775  qn = 0.5;
2776  un = (3.0 / (p_x[n-1] - p_x[n-2])) * (p_clampedDerivLastPt - (p_y[n-1] - p_y[n-2]) /
2777  (p_x[n-1] - p_x[n-2]));
2778  }
2779  p_clampedSecondDerivs[n-1] = (un - qn * u[n-2]) / (qn * p_clampedSecondDerivs[n-2] + 1.0);
2780  for(int i = n - 2; i >= 0; i--) { // backsubstitution loop of the tridiagonal algorithm
2781  p_clampedSecondDerivs[i] = p_clampedSecondDerivs[i] * p_clampedSecondDerivs[i+1] + u[i];
2782  }
2783  p_clampedComputed = true;
2784  return;
2785  }
2786 
2824  double NumericalApproximation::ValueToExtrapolate(const double a, const ExtrapType &etype) {
2825  if(etype == NumericalApproximation::ThrowError) {
2826  ReportException(IException::Programmer, "Evaluate()",
2827  "Invalid argument. Value entered, a = "
2828  + IString(a) + ", is outside of domain = ["
2829  + IString(DomainMinimum()) + ", "
2830  + IString(DomainMaximum()) + "]",
2831  _FILEINFO_);
2832  }
2833  if(etype == NumericalApproximation::NearestEndpoint) {
2834  if(a + DBL_EPSILON < DomainMinimum()) {
2835  return DomainMinimum();
2836  }
2837  else {
2838  return DomainMaximum(); // (a > DomainMaximum())
2839  }
2840  }
2841  else { // gsl interpolations and CubicNeighborhood cannot extrapolate
2842  if(GslInterpType(p_itype)
2843  || p_itype == NumericalApproximation::CubicNeighborhood) {
2844  ReportException(IException::Programmer, "Evaluate()",
2845  "Invalid argument. Cannot extrapolate for type "
2846  + Name() + ", must choose to throw error or return nearest neighbor",
2847  _FILEINFO_);
2848  }
2849  return a;
2850  }
2851  }
2852 
2887  double NumericalApproximation::EvaluateCubicNeighborhood(const double a) {
2888  try {
2889  int s = 0, s_old, s0;
2890  vector <double> x0(4), y0(4);
2891  NumericalApproximation spline(NumericalApproximation::CubicNatural);
2892  for(unsigned int n = 0; n < Size(); n++) {
2893  if(p_x[n] < a) s = n;
2894  else break;
2895  }
2896  if(s < 1) s = 1;
2897  else if(s > (int) Size() - 3) s = Size() - 3;
2898  s_old = -1;
2899  s0 = s - 1;
2900  if(s_old != s0) {
2901  for(int n = 0; n < 4; n++) {
2902  x0[n] = p_x[n+s0];
2903  y0[n] = p_y[n+s0];
2904  spline.AddData(x0[n], y0[n]);
2905  }
2906  s_old = s0;
2907  }
2908  // use nearest endpoint extrapolation method for neighborhood spline
2909  // since CubicNatural can do no other
2910  return spline.Evaluate(a, NumericalApproximation::NearestEndpoint);
2911  }
2912  catch(IException &e) { // catch exception from Constructor, ComputeGsl(), Evaluate()
2913  throw IException(e,
2914  e.errorType(),
2915  "EvaluateCubicNeighborhood() - Unable to evaluate cubic neighborhood interpolation at a = "
2916  + IString(a),
2917  _FILEINFO_);
2918  }
2919  }
2920 
2967  vector <double> NumericalApproximation::EvaluateCubicNeighborhood(const vector <double> &a,
2968  const NumericalApproximation::ExtrapType &etype) {
2969  vector <double> result(a.size());
2970  int s_old, s0;
2971  vector <double> x0(4), y0(4);
2972  vector <int> s;
2973  s.clear();
2974  s.resize(a.size());
2975  for(unsigned int i = 0; i < a.size(); i++) {
2976  for(unsigned int n = 0; n < Size(); n++) {
2977  if(p_x[n] < a[i]) {
2978  s[i] = n;
2979  }
2980  }
2981  if(s[i] < 1) {
2982  s[i] = 1;
2983  }
2984  if(s[i] > ((int) Size()) - 3) {
2985  s[i] = Size() - 3;
2986  }
2987  }
2988  s_old = -1;
2989  try {
2990  NumericalApproximation spline(NumericalApproximation::CubicNatural);
2991  for(unsigned int i = 0; i < a.size(); i++) {
2992  s0 = s[i] - 1;
2993  if(s_old != s0) {
2994  spline.Reset();
2995  for(int n = 0; n < 4; n++) {
2996  x0[n] = p_x[n+s0];
2997  y0[n] = p_y[n+s0];
2998  spline.AddData(x0[n], y0[n]);
2999  }
3000  s_old = s0;
3001  }
3002  double a0;
3003  // checks whether this value is in domain of the main spline
3004  if(InsideDomain(a[i]))
3005  a0 = a[i];
3006  else a0 = ValueToExtrapolate(a[i], etype);
3007  // since neighborhood spline is CubicNatural, the only extrapolation possible is NearestEndpoint
3008  result[i] = spline.Evaluate(a0, NumericalApproximation::NearestEndpoint);
3009  }
3010  return result;
3011  }
3012  catch(IException &e) { // catch exception from Constructor, ComputeGsl(), Evaluate()
3013  throw IException(e,
3014  e.errorType(),
3015  "EvaluateCubicNeighborhood() - Unable to evaluate the function at the given vector of points using cubic neighborhood interpolation",
3016  _FILEINFO_);
3017  }
3018  }
3019 
3054  double NumericalApproximation::EvaluateCubicClamped(double a) {
3055  // This method was derived from an algorithm in the text
3056  // Numerical Recipes in C: The Art of Scientific Computing
3057  // Section 3.3 by Flannery, Press, Teukolsky, and Vetterling
3058  int n = Size();
3059  double result = 0;
3060 
3061  int k;
3062  int k_Lo;
3063  int k_Hi;
3064  double h;
3065  double A;
3066  double B;
3067 
3068  k_Lo = 0;
3069  k_Hi = n - 1;
3070  while(k_Hi - k_Lo > 1) {
3071  k = (k_Hi + k_Lo) / 2;
3072  if(p_x[k] > a) {
3073  k_Hi = k;
3074  }
3075  else {
3076  k_Lo = k;
3077  }
3078  }
3079 
3080  h = p_x[k_Hi] - p_x[k_Lo];
3081  A = (p_x[k_Hi] - a) / h;
3082  B = (a - p_x[k_Lo]) / h;
3083  result = A * p_y[k_Lo] + B * p_y[k_Hi] + ((pow(A, 3.0) - A) *
3084  p_clampedSecondDerivs[k_Lo] + (pow(B, 3.0) - B) * p_clampedSecondDerivs[k_Hi]) * pow(h, 2.0) / 6.0;
3085  return result;
3086  }
3087 
3115  double NumericalApproximation::EvaluateCubicHermite(const double a) {
3116  // algorithm was found at en.wikipedia.org/wiki/Cubic_Hermite_spline
3117  // it seems to produce same answers, as the NumericalAnalysis book
3118 
3119  if(p_fprimeOfx.size() != Size()) {
3120  ReportException(IException::User, "EvaluateCubicHermite()",
3121  "Invalid arguments. The size of the first derivative vector does not match the number of (x,y) data points.",
3122  _FILEINFO_);
3123  }
3124  // find the interval in which "a" exists
3125  int lowerIndex = FindIntervalLowerIndex(a);
3126 
3127  // we know that "a" is within the domain since this is verified in
3128  // Evaluate() before this method is called, thus n <= Size()
3129  if(a == p_x[lowerIndex]) {
3130  return p_y[lowerIndex];
3131  }
3132  if(a == p_x[lowerIndex+1]) {
3133  return p_y[lowerIndex+1];
3134  }
3135 
3136  double x0, x1, y0, y1, m0, m1;
3137  // a is contained within the interval (x0,x1)
3138  x0 = p_x[lowerIndex];
3139  x1 = p_x[lowerIndex+1];
3140  // the corresponding known y-values for x0 and x1
3141  y0 = p_y[lowerIndex];
3142  y1 = p_y[lowerIndex+1];
3143  // the corresponding known tangents (slopes) at (x0,y0) and (x1,y1)
3144  m0 = p_fprimeOfx[lowerIndex];
3145  m1 = p_fprimeOfx[lowerIndex+1];
3146 
3147 
3148  // following algorithm found at en.wikipedia.org/wiki/Cubic_Hermite_spline
3149  // seems to produce same answers, is it faster?
3150 
3151  double h, t;
3152  h = x1 - x0;
3153  t = (a - x0) / h;
3154  return (2 * t * t * t - 3 * t * t + 1) * y0 + (t * t * t - 2 * t * t + t) * h * m0 + (-2 * t * t * t + 3 * t * t) * y1 + (t * t * t - t * t) * h * m1;
3155  }
3156 
3177  int NumericalApproximation::FindIntervalLowerIndex(const double a) {
3178  if(InsideDomain(a)) {
3179  // find the interval in which "a" exists
3180  std::vector<double>::iterator pos;
3181  // find position in vector that is greater than or equal to "a"
3182  pos = upper_bound(p_x.begin(), p_x.end(), a);
3183  int upperIndex = 0;
3184  if(pos != p_x.end()) {
3185  upperIndex = distance(p_x.begin(), pos);
3186  }
3187  else {
3188  upperIndex = Size() - 1;
3189  }
3190  return upperIndex - 1;
3191  }
3192  else if((a + DBL_EPSILON) < DomainMinimum()) {
3193  return 0;
3194  }
3195  else {
3196  return Size() - 2;
3197  }
3198  }
3199 
3200 
3236  double NumericalApproximation::EvaluatePolynomialNeville(const double a) {
3237  // This method was derived from an algorithm in the text
3238  // Numerical Recipes in C: The Art of Scientific Computing
3239  // Section 3.1 by Flannery, Press, Teukolsky, and Vetterling
3240  int n = Size();
3241  double y;
3242 
3243  int ns;
3244  double den, dif, dift, c[n], d[n], ho, hp, w;
3245  double *err = 0;
3246  ns = 1;
3247  dif = fabs(a - p_x[0]);
3248  for(int i = 0; i < n; i++) { // Get the index ns of the closest table entry
3249  dift = fabs(a - p_x[i]);
3250  if(dift < dif) {
3251  ns = i + 1;
3252  dif = dift;
3253  }
3254  c[i] = p_y[i]; // initialize c and d
3255  d[i] = p_y[i];
3256  }
3257  ns = ns - 1;
3258  y = p_y[ns]; // initial approximation for y
3259  for(int m = 1; m < n; m++) {
3260  for(int i = 1; i <= n - m; i++) { // loop over c and d and update them
3261  ho = p_x[i-1] - a;
3262  hp = p_x[i+m-1] - a;
3263  w = c[i] - d[i-1];
3264  den = ho - hp;
3265  den = w / den;
3266  d[i-1] = hp * den; // update c and d
3267  c[i-1] = ho * den;
3268  }
3269  if(2 * ns < n - m) { // After each column in the tableau is completed, we decide
3270  err = &c[ns]; // which correction, c or d, we want to add to our accumulating
3271  } // value of y, i.e., which path to take through the tableau|
3272  else { // forking up or down. We do this in such a way as to take the
3273  ns = ns - 1; // most "straight line" route through the tableau to its apex,
3274  err = &d[ns]; // updating ns accordingly to keep track of where we are. This
3275  } // route keeps the partial approximations centered (insofar as possible)
3276  y = y + *err; // on the target x. The last err added is thus the error indication.
3277  }
3278  p_polyNevError.push_back(*err);
3279  return y;
3280  }
3281 
3300  vector <double> NumericalApproximation::EvaluateForIntegration(const double a, const double b, const unsigned int n) {
3301  if(a > b) {
3302  ReportException(IException::Programmer, "EvaluatedForIntegration()",
3303  "Invalid interval entered: [a,b] = ["
3304  + IString(a) + ", " + IString(b) + "]",
3305  _FILEINFO_);
3306  }
3307  if(!InsideDomain(a) || !InsideDomain(b)) {
3308  ReportException(IException::Programmer, "EvaluateForIntegration()",
3309  "Invalid arguments. Interval entered ["
3310  + IString(a) + ", " + IString(b)
3311  + "] is not contained within domain ["
3312  + IString(DomainMinimum()) + ", "
3313  + IString(DomainMaximum()) + "]",
3314  _FILEINFO_);
3315  }
3316  vector <double> f;
3317  //need total number of segments to be divisible by n-1
3318  // This way we can use the formula for each interval: 0 to n, n to 2n, 2n to 3n, etc.
3319  // Notice interval endpoints will overlap for these composite formulas
3320  int xSegments = Size() - 1;
3321  while(xSegments % (n - 1) != 0) {
3322  xSegments++;
3323  }
3324  // x must be sorted and unique
3325  double xMin = a;
3326  double xMax = b;
3327  //Uniform step size.
3328  double h = (xMax - xMin) / xSegments;
3329  //Compute the interpolates using spline.
3330  try {
3331  for(double i = 0; i < (xSegments + 1); i++) {
3332  double xi = h * i + xMin;
3333  f.push_back(Evaluate(xi)); // validate data set here, allow default ThrowError extrap type
3334  }
3335  f.push_back(h);
3336  return f;
3337  }
3338  catch(IException &e) { // catch exception from Evaluate()
3339  throw IException(e,
3340  e.errorType(),
3341  "EvaluateForIntegration() - Unable to evaluate the data set for integration",
3342  _FILEINFO_);
3343  }
3344  }// for integration using composite of spline (creates evenly spaced points)
3345 
3366  void NumericalApproximation::ReportException(IException::ErrorType type, const string &methodName, const string &message,
3367  const char *filesrc, int lineno)
3368  const {
3369  string msg = methodName + " - " + message;
3370  throw IException(type, msg.c_str(), filesrc, lineno);
3371  return;
3372  }
3373 
3374 }
bool p_clampedEndptsSet
Flag variable to determine whether SetCubicClampedEndptDeriv() has been called after all data was add...
vector< double > p_x
List of X values.
InterpType
This enum defines the types of interpolation supported in this class.
vector< double > p_fprimeOfx
List of first derivatives corresponding to each x value in the data set (i.e. each value in p_x) ...
double p_clampedDerivFirstPt
First derivative of first x-value, p_x[0]. This is only used for the CubicClamped interpolation type...
Namespace for the standard library.
NumericalApproximation provides various numerical analysis methods of interpolation, extrapolation and approximation of a tabulated set of x, y data.
vector< double > p_clampedSecondDerivs
List of second derivatives evaluated at p_x values. This is only used for the CubicClamped interpolat...
double Evaluate(const double a, const ExtrapType &etype=ThrowError)
Calculates interpolated or extrapolated value of tabulated data set for given domain value...
void Reset()
Resets the state of the object.
vector< double > PolynomialNevilleErrorEstimate()
Retrieves the error estimate for the Neville&#39;s polynomial interpolation type.
ExtrapType
This enum defines the manner in which a value outside of the domain should be handled if passed to th...
#define _FILEINFO_
Macro for the filename and line number.
Definition: IException.h:40
double p_clampedDerivLastPt
First derivative of last x-value, p_x[n-1]. This is only used for the CubicClamped interpolation type...
const gsl_interp_type * InterpFunctor
GSL Interpolation specs.
InterpFunctor GslFunctor(NumericalApproximation::InterpType itype) const
Search for a GSL interpolation function.
ErrorType errorType() const
Returns the source of the error for this exception.
Definition: IException.cpp:446
ErrorType
Contains a set of exception error types.
Definition: IException.h:127
void AddData(const double x, const double y)
Add a datapoint to the set.
bool p_dataValidated
Flag variable to determine whether ValidateDataSet() has been called.
Isis exception class.
Definition: IException.h:107
Adds specific functionality to C++ strings.
Definition: IString.h:181
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
vector< double > p_y
List of Y values.
vector< double > p_polyNevError
Estimate of error for interpolation evaluated at x. This is only used for the PolynomialNeville inter...
InterpType p_itype
Interpolation type.
FunctorList::const_iterator FunctorConstIter
GSL Iterator.
bool p_clampedComputed
Flag variable to determine whether ComputeCubicClamped() has been called.