NumericalApproximation.cpp

#include "NumericalApproximation.h"

#include <cmath>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>

#include "SpecialPixel.h"
#include "IException.h"
#include "IString.h"

using namespace std;
namespace Isis {
  NumericalApproximation::FunctorList NumericalApproximation::p_interpFunctors;

  NumericalApproximation::NumericalApproximation(const NumericalApproximation::InterpType &itype) {
    // Validate parameter and initialize class variables
    try {
      Init(itype);
      p_dataValidated = false;
    }
    catch(IException &e) {  // catch exception from Init()
      throw IException(e, e.errorType(),
                       "NumericalApproximation() - Unable to construct NumericalApproximation object",
                       _FILEINFO_);
    }
  }

  NumericalApproximation::NumericalApproximation(unsigned int n, double *x, double *y,
      const NumericalApproximation::InterpType &itype) {
    try {
      Init(itype);
      AddData(n, x, y);
      ValidateDataSet();
    }
    catch(IException &e) { // catch exception from Init(), AddData(), ValidateDataSet()
      throw IException(e, e.errorType(),
                       "NumericalApproximation() - Unable to construct object using the given arrays, size and interpolation type",
                       _FILEINFO_);
    }
  }

  NumericalApproximation::NumericalApproximation(const vector <double> &x, const vector <double> &y,
      const NumericalApproximation::InterpType &itype) {
    try {
      Init(itype);
      AddData(x, y); // size of x = size of y validated in this method
      ValidateDataSet();
    }
    catch(IException &e) { // catch exception from Init(), AddData(), ValidateDataSet()
      throw IException(e, e.errorType(),
                       "NumericalApproximation() - Unable to construct an object using the given vectors and interpolation type",
                       _FILEINFO_);
    }
  }

  NumericalApproximation::NumericalApproximation(const NumericalApproximation &oldObject) {
    try {
      // initialize new object and set type to that of old object
      Init(oldObject.p_itype);
      // fill data set of new object
      p_x = oldObject.p_x;
      p_y = oldObject.p_y;
      // if this data was previously validated, no need to do this again
      p_dataValidated = oldObject.p_dataValidated;
      // copy values for interpolation-specific variables
      p_clampedComputed = oldObject.p_clampedComputed;
      p_clampedEndptsSet = oldObject.p_clampedEndptsSet;
      p_clampedSecondDerivs = oldObject.p_clampedSecondDerivs;
      p_clampedDerivFirstPt = oldObject.p_clampedDerivFirstPt;
      p_clampedDerivLastPt = oldObject.p_clampedDerivLastPt;
      p_polyNevError = oldObject.p_polyNevError;
      p_fprimeOfx = oldObject.p_fprimeOfx;
    }
    catch(IException &e) { // catch exception from Init()
      throw IException(e,
                       e.errorType(),
                       "NumericalApproximation() - Unable to copy the given object",
                       _FILEINFO_);
    }
  }

  NumericalApproximation &NumericalApproximation::operator=(const NumericalApproximation &oldObject) {
    try {
      if(&oldObject != this) {
        if(GslInterpType(p_itype)) GslDeallocation();
        SetInterpType(oldObject.p_itype);
        // set class variables
        p_x = oldObject.p_x;
        p_y = oldObject.p_y;
        p_dataValidated = oldObject.p_dataValidated;
        p_clampedSecondDerivs = oldObject.p_clampedSecondDerivs;
        p_clampedDerivFirstPt = oldObject.p_clampedDerivFirstPt;
        p_clampedDerivLastPt = oldObject.p_clampedDerivLastPt;
        p_polyNevError = oldObject.p_polyNevError;
        p_fprimeOfx = oldObject.p_fprimeOfx;
      }
      return (*this);
    }
    catch(IException &e) { // catch exception from SetInterpType()
      throw IException(e,
                       e.errorType(),
                       "operator=() - Unable to copy the given object",
                       _FILEINFO_);
    }

  }

  NumericalApproximation::~NumericalApproximation() {
    if(GslInterpType(p_itype)) GslDeallocation();
  }

  string NumericalApproximation::Name() const {
    if(p_itype == NumericalApproximation::CubicNeighborhood) {
      return "cspline-neighborhood";
    }
    else if(p_itype == NumericalApproximation::CubicClamped) {
      return "cspline-clamped";
    }
    else if(p_itype == NumericalApproximation::PolynomialNeville) {
      return "polynomial-Neville's";
    }
    else if(p_itype == NumericalApproximation::CubicHermite) {
      return "cspline-Hermite";
    }
    try {
      string name = (string(GslFunctor(p_itype)->name));
      if(p_itype == NumericalApproximation::CubicNatural) {
        return name + "-natural";
      }
      else return name;
    }
    catch(IException &e) { // catch exception from GslFunctor()
      throw IException(e,
                       e.errorType(),
                       "Name() - GSL interpolation type not found",
                       _FILEINFO_);
    }
  }

  int NumericalApproximation::MinPoints() {
    if(p_itype == NumericalApproximation::CubicNeighborhood) {
      return 4;
    }
    else if(p_itype == NumericalApproximation::CubicClamped) {
      return 3;
    }
    else if(p_itype == NumericalApproximation::PolynomialNeville) {
      return 3;
    }
    else if(p_itype == NumericalApproximation::CubicHermite) {
      return 2;
    }
    try {
      return (GslFunctor(p_itype)->min_size);
    }
    catch(IException &e) { // catch exception from GslFunctor()
      throw IException(e,
                       e.errorType(),
                       "MinPoints() - GSL interpolation not found",
                       _FILEINFO_);
    }
  }

  int NumericalApproximation::MinPoints(NumericalApproximation::InterpType itype) {
    //Validates the parameter
    if(GslInterpType(itype)) {
      try {
        //GslFunctor() Validates parameter for Gsl interp types
        NumericalApproximation nam(itype);
        return (nam.GslFunctor(itype)->min_size);
      }
      catch(IException &e) { // catch exception from GslFunctor()
        throw IException(e,
                         e.errorType(),
                         "MinPoints() - GSL interpolation type not found",
                         _FILEINFO_);
      }
    }
    else if(itype == NumericalApproximation::CubicNeighborhood) {
      return 4;
    }
    else if(itype == NumericalApproximation::CubicClamped) {
      return 3;
    }
    else if(itype == NumericalApproximation::PolynomialNeville) {
      return 3;
    }
    else if(itype == NumericalApproximation::CubicHermite) {
      return 2;
    }
    throw IException(IException::Programmer,
                     "MinPoints() - Invalid argument. Unknown interpolation type: "
                     + IString(NumericalApproximation::InterpType(itype)),
                     _FILEINFO_);
  }

  void NumericalApproximation::AddData(const double x, const double y) {
    p_x.push_back(x);
    p_y.push_back(y);
    p_clampedComputed = false;
    p_clampedEndptsSet = false;
    p_dataValidated = false;
    p_clampedSecondDerivs.clear();
    p_clampedDerivFirstPt = 0;
    p_clampedDerivLastPt = 0;
    p_polyNevError.clear();
    p_interp = 0;
    p_acc = 0;
    return;
  }

  void NumericalApproximation::AddData(unsigned int n, double *x, double *y) {
    for(unsigned int i = 0; i < n; i++) {
      p_x.push_back(x[i]);
      p_y.push_back(y[i]);
    }
    p_clampedComputed = false;
    p_clampedEndptsSet = false;
    p_dataValidated = false;
    p_clampedSecondDerivs.clear();
    p_clampedDerivFirstPt = 0;
    p_clampedDerivLastPt = 0;
    p_polyNevError.clear();
    p_interp = 0;
    p_acc = 0;
    return;
  }

  void NumericalApproximation::AddData(const vector <double> &x,
                                       const vector <double> &y)
                                        {
    int n = x.size();
    int m = y.size();
    if(n != m) {
      ReportException(IException::Programmer, "AddData()",
                      "Invalid arguments. The sizes of the input vectors do "
                      "not match",
                      _FILEINFO_);
    }

    // Avoid push_back if at all possible. These calls were consuming 10% of
    //   cam2map's run time on a line scan camera.
    if(!p_x.empty() || !p_y.empty()) {
      for(int i = 0; i < n; i++) {
        p_x.push_back(x[i]);
        p_y.push_back(y[i]);
      }
    }
    else {
      p_x = x;
      p_y = y;
    }

    p_clampedComputed = false;
    p_clampedEndptsSet = false;
    p_dataValidated = false;
    p_clampedSecondDerivs.clear();
    p_clampedDerivFirstPt = 0;
    p_clampedDerivLastPt = 0;
    p_polyNevError.clear();
    p_interp = 0;
    p_acc = 0;
    return;
  }

  void NumericalApproximation::SetCubicClampedEndptDeriv(double yp1, double ypn) {
    if(p_itype != NumericalApproximation::CubicClamped) {
      ReportException(IException::Programmer, "SetCubicClampedEndptDeriv()",
                      "This method is only valid for cspline-clamped interpolation, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    p_clampedDerivFirstPt = yp1;
    p_clampedDerivLastPt = ypn;
    p_clampedEndptsSet = true;
    return;
  }

  void NumericalApproximation::AddCubicHermiteDeriv(unsigned int n, double *fprimeOfx) {
    if(p_itype != NumericalApproximation::CubicHermite) {
      ReportException(IException::Programmer, "SetCubicHermiteDeriv()",
                      "This method is only valid for cspline-Hermite interpolation, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    for(unsigned int i = 0; i < n; i++) {
      p_fprimeOfx.push_back(fprimeOfx[i]);
    }
    return;
  }
  void NumericalApproximation::AddCubicHermiteDeriv(
      const vector <double> &fprimeOfx) {
    if(p_itype != NumericalApproximation::CubicHermite) {
      ReportException(IException::Programmer, "SetCubicHermiteDeriv()",
                      "This method is only valid for cspline-Hermite interpolation, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }

    // Avoid push_back if at all possible.
    if(!p_fprimeOfx.empty()) {
      for(unsigned int i = 0; i < fprimeOfx.size(); i++) {
        p_fprimeOfx.push_back(fprimeOfx[i]);
      }
    }
    else {
      p_fprimeOfx = fprimeOfx;
    }

    return;
  }
  void NumericalApproximation::AddCubicHermiteDeriv(
      double fprimeOfx) {
    if(p_itype != NumericalApproximation::CubicHermite) {
      ReportException(IException::Programmer, "SetCubicHermiteDeriv()",
                      "This method is only valid for cspline-Hermite interpolation, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    p_fprimeOfx.push_back(fprimeOfx);
    return;
  }

  vector <double> NumericalApproximation::CubicClampedSecondDerivatives() {
    try {
      if(p_itype != NumericalApproximation::CubicClamped)
        ReportException(IException::Programmer, "CubicClampedSecondDerivatives()",
                        "This method is only valid for cspline-clamped interpolation type may not be used for "
                        + Name() + " interpolation",
                        _FILEINFO_);
      if(!p_clampedComputed) ComputeCubicClamped();
      return p_clampedSecondDerivs;
    }
    catch(IException &e) {  // catch exception from ComputeCubicClamped()
      throw IException(e,
                       e.errorType(),
                       "CubicClampedSecondDerivatives() - Unable to compute clamped cubic spline interpolation",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::DomainMinimum() {
    try {
      if(GslInterpType(p_itype)) {
        // if GSL interpolation type, we need to compute, if not already done
        if(!GslComputed()) ComputeGsl();
        return (p_interp->interp->xmin);
      }
      else {
        if(!p_dataValidated) ValidateDataSet();
        return *min_element(p_x.begin(), p_x.end());
      }
    }
    catch(IException &e) { // catch exception from ComputeGsl() or ValidateDataSet()
      throw IException(e,
                       e.errorType(),
                       "DomainMinimum() - Unable to calculate the domain minimum for the data set",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::DomainMaximum() {
    try {
      if(GslInterpType(p_itype)) {
        // if GSL interpolation type, we need to compute, if not already done
        if(!GslComputed()) ComputeGsl();
        return (p_interp->interp->xmax);
      }
      else {
        if(!p_dataValidated) ValidateDataSet();
        return *max_element(p_x.begin(), p_x.end());
      }
    }
    catch(IException &e) { // catch exception from ComputeGsl() or ValidateDataSet()
      throw IException(e,
                       e.errorType(),
                       "DomainMaximum() - Unable to calculate the domain maximum for the data set",
                       _FILEINFO_);
    }
  }

  bool NumericalApproximation::Contains(double x) {
    return binary_search(p_x.begin(), p_x.end(), x);
  }

  double NumericalApproximation::Evaluate(const double a, const ExtrapType &etype) {
    try {
      // a is const, so we must set up temporary variable in case value needs to be changed
      double a0;
      if(InsideDomain(a))  // this will validate data set, if not already done
        a0 = a;
      else a0 = ValueToExtrapolate(a, etype);
      // perform interpolation/extrapoltion
      if(p_itype == NumericalApproximation::CubicNeighborhood) {
        return EvaluateCubicNeighborhood(a0);
      }
      else if(p_itype == NumericalApproximation::PolynomialNeville) {
        p_polyNevError.clear();
        return EvaluatePolynomialNeville(a0);
      }
      else if(p_itype == NumericalApproximation::CubicClamped) {
        // if cubic clamped, we need to compute, if not already done
        if(!p_clampedComputed) ComputeCubicClamped();
        return EvaluateCubicClamped(a0);
      }
      else if(p_itype == NumericalApproximation::CubicHermite) {
        return EvaluateCubicHermite(a0);
      }
      // if GSL interpolation type we need to compute, if not already done
      if(!GslComputed()) ComputeGsl();
      double result;
      GslIntegrityCheck(gsl_spline_eval_e(p_interp, a0, p_acc, &result), _FILEINFO_);
      return (result);
    }
    catch(IException &e) { // catch exception from EvaluateCubicNeighborhood(), EvaluateCubicClamped(), EvaluatePolynomialNeville(), GslIntegrityCheck()
      throw IException(e,
                       e.errorType(),
                       "Evaluate() - Unable to evaluate the function at the point a = "
                       + IString(a),
                       _FILEINFO_);
    }
  }

  vector <double> NumericalApproximation::Evaluate(const vector <double> &a, const ExtrapType &etype) {
    try {
      if(p_itype == NumericalApproximation::CubicNeighborhood) {
        // cubic neighborhood has it's own method that will take entire vector
        // this is faster than looping through values calling Evaluate(double)
        // for each component of the passed in vector
        return EvaluateCubicNeighborhood(a, etype);
      }
      vector <double> result(a.size());
      if(p_itype == NumericalApproximation::PolynomialNeville) {
        // cannot loop through values calling Evaluate(double)
        // because it will clear the p_polyNevError vector each time
        p_polyNevError.clear();
        for(unsigned int i = 0; i < result.size(); i++) {
          double a0;
          if(InsideDomain(a[i]))
            a0 = a[i];
          else a0 = ValueToExtrapolate(a[i], etype);
          result[i] = EvaluatePolynomialNeville(a0);
        }
        return result;
      }
      // cubic-clamped, cubic-Hermite and gsl types can be done by calling Evaluate(double)
      // for each value of the passed in vector
      for(unsigned int i = 0; i < result.size(); i++) {
        result[i] = Evaluate(a[i], etype);
      }
      return result;
    }
    catch(IException &e) { // catch exception from EvaluateCubicNeighborhood(), EvaluateCubicClamped(), EvaluatePolynomialNeville(), GslIntegrityCheck()
      throw IException(e,
                       e.errorType(),
                       "Evaluate() - Unable to evaluate the function at the given vector of points",
                       _FILEINFO_);
    }
  }

  vector <double> NumericalApproximation::PolynomialNevilleErrorEstimate() {
    if(p_itype != NumericalApproximation::PolynomialNeville) {
      ReportException(IException::Programmer, "PolynomialNevilleErrorEstimate()",
                      "This method is only valid for polynomial-Neville's, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    if(p_polyNevError.empty()) {
      ReportException(IException::Programmer, "PolynomialNevilleErrorEstimate()",
                      "Error not calculated. This method only valid after Evaluate() has been called",
                      _FILEINFO_);
    }
    return p_polyNevError;
  }

  double NumericalApproximation::GslFirstDerivative(const double a) {
    try { // we need to compute, if not already done
      if(!GslComputed()) ComputeGsl();
    }
    catch(IException &e) { // catch exception from ComputeGsl()
      throw IException(e,
                       e.errorType(),
                       "GslFirstDerivative() - Unable to compute the first derivative at a = "
                       + IString(a) + " using the GSL interpolation",
                       _FILEINFO_);
    }
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "GslFirstDerivative()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!GslInterpType(p_itype)) {
      ReportException(IException::Programmer, "GslFirstDerivative()",
                      "Method only valid for GSL interpolation types, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    try {
      double value;
      GslIntegrityCheck(gsl_spline_eval_deriv_e(p_interp, a, p_acc, &value), _FILEINFO_);
      return (value);
    }
    catch(IException &e) { // catch exception from GslIntegrityCheck()
      throw IException(e,
                       e.errorType(),
                       "GslFirstDerivative() - Unable to compute the first derivative at a = "
                       + IString(a) + ". GSL integrity check failed",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::EvaluateCubicHermiteFirstDeriv(const double a) {
    if(p_itype != NumericalApproximation::CubicHermite) { //??? is this necessary??? create single derivative method with GSL?
      ReportException(IException::User, "EvaluateCubicHermiteFirstDeriv()",
                      "This method is only valid for cspline-Hermite interpolation, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    if(p_fprimeOfx.size() != Size()) {
      ReportException(IException::User, "EvaluateCubicHermiteFirstDeriv()",
                      "Invalid arguments. The size of the first derivative vector does not match the number of (x,y) data points.",
                      _FILEINFO_);
    }
    // find the interval in which "a" exists
    int lowerIndex = FindIntervalLowerIndex(a);

    // we know that "a" is within the domain since this is verified in
    // Evaluate() before this method is called, thus n <= Size()
    if(a == p_x[lowerIndex]) {
      return p_fprimeOfx[lowerIndex];
    }
    if(a == p_x[lowerIndex+1]) {
      return p_fprimeOfx[lowerIndex+1];
    }

    double x0, x1, y0, y1, m0, m1;
    // a is contained within the interval (x0,x1)
    x0 = p_x[lowerIndex];
    x1 = p_x[lowerIndex+1];
    // the corresponding known y-values for x0 and x1
    y0 = p_y[lowerIndex];
    y1 = p_y[lowerIndex+1];
    // the corresponding known tangents (slopes) at (x0,y0) and (x1,y1)
    m0 = p_fprimeOfx[lowerIndex];
    m1 = p_fprimeOfx[lowerIndex+1];

    double h, t;
    h = x1 - x0;
    t = (a - x0) / h;
    if(h != 0.) {
      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;
    }
    else {
      return 0;  // Should never happen
    }
  }

  double NumericalApproximation::BackwardFirstDifference(const double a, const unsigned int n, const double h) {
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "BackwardFirstDifference()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a - (n - 1)*h)) {
      ReportException(IException::Programmer, "BackwardFirstDifference()",
                      "Formula steps outside of domain. For "
                      + IString((int) n) + "-point backward difference, a-(n-1)h = "
                      + IString(a - (n - 1)*h) + " is smaller than domain min = "
                      + IString(DomainMinimum())
                      + ".  Try forward difference or use smaller value for h or n",
                      _FILEINFO_);
    }
    if(!p_dataValidated) ValidateDataSet();
    vector <double> f;
    double xi;
    try {
      for(double i = 0; i < n; i++) {
        xi = a + h * (i - (n - 1));
        f.push_back(Evaluate(xi)); // allow ExtrapType = ThrowError (default)
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "BackwardFirstDifference() - Unable to calculate backward first difference for (a, n, h) = ("
                       + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
                       _FILEINFO_);
    }
    switch(n) {
      case 2:
        return (-f[0] + f[1]) / h;              //2pt backward
      case 3:
        return (3 * f[2] - 4 * f[1] + f[0]) / (2 * h); //3pt backward
      default:
        throw IException(IException::Programmer,
                         "BackwardFirstDifference() - Invalid argument. There is no "
                         + IString((int) n) + "-point backward difference formula in use",
                         _FILEINFO_);
    }
  }


  double NumericalApproximation::ForwardFirstDifference(const double a, const unsigned int n, const double h) {
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "ForwardFirstDifference()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a + (n - 1)*h)) {
      ReportException(IException::Programmer, "ForwardFirstDifference()",
                      "Formula steps outside of domain. For "
                      + IString((int) n) + "-point forward difference, a+(n-1)h = "
                      + IString(a + (n - 1)*h) + " is greater than domain max = "
                      + IString(DomainMaximum())
                      + ".  Try backward difference or use smaller value for h or n",
                      _FILEINFO_);
    }
    if(!p_dataValidated) ValidateDataSet();
    vector <double> f;
    double xi;
    try {
      for(double i = 0; i < n; i++) {
        xi = a + h * i;
        f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "ForwardFirstDifference() - Unable to calculate forward first difference for (a, n, h) = ("
                       + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
                       _FILEINFO_);
    }
    switch(n) {
      case 2:
        return (-f[0] + f[1]) / h;              //2pt forward
      case 3:
        return (-3 * f[0] + 4 * f[1] - f[2]) / (2 * h); //3pt forward
      default:
        throw IException(IException::Programmer,
                         "ForwardFirstDifference() - Invalid argument. There is no "
                         + IString((int) n) + "-point forward difference formula in use",
                         _FILEINFO_);
    }
  }


  double NumericalApproximation::CenterFirstDifference(const double a, const unsigned int n, const double h) {
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "CenterFirstDifference()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a + (n - 1)*h) || !InsideDomain(a - (n - 1)*h)) {
      ReportException(IException::Programmer, "CenterFirstDifference()",
                      "Formula steps outside of domain. For "
                      + IString((int) n) + "-point center difference, a-(n-1)h = "
                      + IString(a - (n - 1)*h) + " or a+(n-1)h = "
                      + IString(a + (n - 1)*h) + " is out of domain = ["
                      + IString(DomainMinimum()) + ", " + IString(DomainMaximum())
                      + "].  Use smaller value for h or n",
                      _FILEINFO_);
    }
    if(!p_dataValidated) ValidateDataSet();
    vector <double> f;
    double xi;
    try {
      for(double i = 0; i < n; i++) {
        xi = a + h * (i - (n - 1) / 2);
        f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "CenterFirstDifference() - Unable to calculate center first difference for (a, n, h) = ("
                       + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
                       _FILEINFO_);
    }
    switch(n) {
      case 3:
        return (-f[0] + f[2]) / (2 * h);              //3pt center
      case 5:
        return (f[0] - 8 * f[1] + 8 * f[3] - f[4]) / (12 * h); //5pt center
      default:
        throw IException(IException::Programmer,
                         "CenterFirstDifference() - Invalid argument. There is no "
                         + IString((int) n) + "-point center difference formula in use",
                         _FILEINFO_);
    }
  }

  double NumericalApproximation::GslSecondDerivative(const double a) {
    try { // we need to compute, if not already done
      if(!GslComputed()) ComputeGsl();
    }
    catch(IException &e) { // catch exception from ComputeGsl()
      throw IException(e,
                       e.errorType(),
                       "GslSecondDerivative() - Unable to compute the second derivative at a = "
                       + IString(a) + " using the GSL interpolation",
                       _FILEINFO_);
    }
    if(!InsideDomain(a))
      ReportException(IException::Programmer, "GslSecondDerivative()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    if(!GslInterpType(p_itype))
      ReportException(IException::Programmer, "GslSecondDerivative()",
                      "Method only valid for GSL interpolation types, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    try {
      // we need to compute, if not already done
      if(!GslComputed()) ComputeGsl();
      double value;
      GslIntegrityCheck(gsl_spline_eval_deriv2_e(p_interp, a, p_acc, &value), _FILEINFO_);
      return (value);
    }
    catch(IException &e) { // catch exception from GslIntegrityCheck()
      throw IException(e,
                       e.errorType(),
                       "GslSecondDerivative() - Unable to compute the second derivative at a = "
                       + IString(a) + ". GSL integrity check failed",
                       _FILEINFO_);
    }
  }


  double NumericalApproximation::EvaluateCubicHermiteSecDeriv(const double a) {
    if(p_itype != NumericalApproximation::CubicHermite) {
      ReportException(IException::User, "EvaluateCubicHermiteSecDeriv()",
                      "This method is only valid for cspline-Hermite interpolation, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    if(p_fprimeOfx.size() != Size()) {
      ReportException(IException::User, "EvaluateCubicHermiteSecDeriv()",
                      "Invalid arguments. The size of the first derivative vector does not match the number of (x,y) data points.",
                      _FILEINFO_);
    }
    // find the interval in which "a" exists
    int lowerIndex = FindIntervalLowerIndex(a);
    double x0, x1, y0, y1, m0, m1;
    // a is contained within the interval (x0,x1)
    x0 = p_x[lowerIndex];
    x1 = p_x[lowerIndex+1];
    // the corresponding known y-values for x0 and x1
    y0 = p_y[lowerIndex];
    y1 = p_y[lowerIndex+1];
    // the corresponding known tangents (slopes) at (x0,y0) and (x1,y1)
    m0 = p_fprimeOfx[lowerIndex];
    m1 = p_fprimeOfx[lowerIndex+1];

    double h, t;
    h = x1 - x0;
    t = (a - x0) / h;
    if(h != 0.) {
      return ((12 * t - 6) * y0 + (6 * t - 4) * h * m0 + (-12 * t + 6) * y1 + (6 * t - 2) * h * m1) / h;
    }
    else {
      return 0; // Should never happen
    }
  }

  double NumericalApproximation::BackwardSecondDifference(const double a, const unsigned int n, const double h) {
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "BackwardSecondDifference()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a - (n - 1)*h)) {
      ReportException(IException::Programmer, "BackwardSecondDifference()",
                      "Formula steps outside of domain. For "
                      + IString((int) n) + "-point backward difference, a-(n-1)h = "
                      + IString(a - (n - 1)*h) + " is smaller than domain min = "
                      + IString(DomainMinimum())
                      + ".  Try forward difference or use smaller value for h or n",
                      _FILEINFO_);
    }
    if(!p_dataValidated) ValidateDataSet();
    vector <double> f;
    double xi;
    try {
      for(double i = 0; i < n; i++) {
        xi = a + h * (i - (n - 1));
        f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "BackwardSecondDifference() - Unable to calculate backward second difference for (a, n, h) = ("
                       + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
                       _FILEINFO_);
    }
    switch(n) {
      case 3:
        return (f[0] - 2 * f[1] + f[2]) / (h * h);                   //3pt backward
      default:
        throw IException(IException::Programmer,
                         "BackwardSecondDifference() - Invalid argument. There is no "
                         + IString((int) n) + "-point backward second difference formula in use",
                         _FILEINFO_);
    }
  }


  double NumericalApproximation::ForwardSecondDifference(const double a, const unsigned int n, const double h) {
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "ForwardSecondDifference()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a + (n - 1)*h)) {
      ReportException(IException::Programmer, "ForwardSecondDifference()",
                      "Formula steps outside of domain. For "
                      + IString((int) n) + "-point forward difference, a+(n-1)h = "
                      + IString(a + (n - 1)*h) + " is greater than domain max = "
                      + IString(DomainMaximum())
                      + ".  Try backward difference or use smaller value for h or n",
                      _FILEINFO_);
    }
    if(!p_dataValidated) ValidateDataSet();
    vector <double> f;
    double xi;
    try {
      for(double i = 0; i < n; i++) {
        xi = a + h * i;
        f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "ForwardSecondDifference() - Unable to calculate forward second difference for (a, n, h) = ("
                       + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
                       _FILEINFO_);
    }
    switch(n) {
      case 3:
        return (f[0] - 2 * f[1] + f[2]) / (h * h);                   //3pt forward
      default:
        throw IException(IException::Programmer,
                         "ForwardSecondDifference() - Invalid argument. There is no "
                         + IString((int) n) + "-point forward second difference formula in use",
                         _FILEINFO_);
    }
  }


  double NumericalApproximation::CenterSecondDifference(const double a, const unsigned int n, const double h) {
    if(!InsideDomain(a)) {
      ReportException(IException::Programmer, "CenterSecondDifference()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a + (n - 1)*h) || !InsideDomain(a - (n - 1)*h)) {
      ReportException(IException::Programmer, "CenterSecondDifference()",
                      "Formula steps outside of domain. For "
                      + IString((int) n) + "-point center difference, a-(n-1)h = "
                      + IString(a - (n - 1)*h) + " or a+(n-1)h = "
                      + IString(a + (n - 1)*h) + " is out of domain = ["
                      + IString(DomainMinimum()) + ", " + IString(DomainMaximum())
                      + "].  Use smaller value for h or n",
                      _FILEINFO_);
    }
    if(!p_dataValidated) ValidateDataSet();
    vector <double> f;
    double xi;
    try {
      for(double i = 0; i < n; i++) {
        xi = a + h * (i - (n - 1) / 2);
        f.push_back(Evaluate(xi));// allow ExtrapType = ThrowError (default)
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "CenterSecondDifference() - Unable to calculate center second difference for (a, n, h) = ("
                       + IString(a) + ", " + IString((int) n) + ", " + IString(h) + ")",
                       _FILEINFO_);
    }
    switch(n) {
      case 3:
        return (f[0] - 2 * f[1] + f[2]) / (h * h);                   //3pt center
      case 5:
        return (-f[0] + 16 * f[1] - 30 * f[2] + 16 * f[3] - f[4]) / (12 * h * h); //5pt center
      default:
        throw IException(IException::Programmer,
                         "CenterSecondDifference() - Invalid argument. There is no "
                         + IString((int) n) + "-point center second difference formula in use",
                         _FILEINFO_);
    }
  }

  double NumericalApproximation::GslIntegral(const double a, const double b) {
    try { // we need to compute, if not already done
      if(!GslComputed()) ComputeGsl();
    }
    catch(IException &e) { // catch exception from ComputeGsl()
      throw IException(e,
                       e.errorType(),
                       "GslIntegral() - Unable to compute the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using the GSL interpolation",
                       _FILEINFO_);
    }
    if(a > b) {
      ReportException(IException::Programmer, "GslIntegral()",
                      "Invalid interval entered: [a,b] = ["
                      + IString(a) + ", " + IString(b) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a) || !InsideDomain(b)) {
      ReportException(IException::Programmer, "GslIntegral()",
                      "Invalid arguments. Interval entered ["
                      + IString(a) + ", " + IString(b)
                      + "] is not contained within domain ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(!GslInterpType(p_itype)) {
      ReportException(IException::Programmer, "GslIntegral()",
                      "Method only valid for GSL interpolation types, may not be used for "
                      + Name() + " interpolation",
                      _FILEINFO_);
    }
    try {
      double value;
      GslIntegrityCheck(gsl_spline_eval_integ_e(p_interp, a, b, p_acc, &value), _FILEINFO_);
      return (value);
    }
    catch(IException &e) { // catch exception from GslIntegrityCheck()
      throw IException(e,
                       e.errorType(),
                       "GslIntegral() - Unable to compute the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + "). GSL integrity check failed",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::TrapezoidalRule(const double a, const double b) {
    try {
      //n is the number of points used in the formula of the integration type chosen
      unsigned int n = 2;
      double result = 0;
      vector <double> f = EvaluateForIntegration(a, b, n);
      double h = f.back();
      f.pop_back();
      //Compute the integral using composite trapezoid formula
      int ii;
      for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
        ii = (i + 1) * (n - 1);
        result += (f[ii-1] + f[ii]) * h / 2;
      }
      return result;
    }
    catch(IException &e) { // catch exception from EvaluateForIntegration()
      throw IException(e,
                       e.errorType(),
                       "TrapezoidalRule() - Unable to calculate the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using the trapeziodal rule",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::Simpsons3PointRule(const double a, const double b) {
    try {
      //n is the number of points used in the formula of the integration type chosen
      unsigned int n = 3;
      double result = 0;
      vector <double> f = EvaluateForIntegration(a, b, n);
      double h = f.back();
      f.pop_back();
      int ii;
      for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
        ii = (i + 1) * (n - 1);
        result += (f[ii-2] + 4 * f[ii-1] + f[ii]) * h / 3;
      }
      return result;
    }
    catch(IException &e) { // catch exception from EvaluateForIntegration()
      throw IException(e,
                       e.errorType(),
                       "Simpsons3PointRule() - Unable to calculate the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using Simpson's 3 point rule",
                       _FILEINFO_);
    }
  }


  double NumericalApproximation::Simpsons4PointRule(const double a, const double b) {
    try {
      //n is the number of points used in the formula of the integration type chosen
      unsigned int n = 4;
      double result = 0;
      vector <double> f = EvaluateForIntegration(a, b, n);
      double h = f.back();
      f.pop_back();
      int ii;
      for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
        ii = (i + 1) * (n - 1);
        result += (f[ii-3] + 3 * f[ii-2] + 3 * f[ii-1] + f[ii]) * h * 3 / 8;
      }
      return result;
    }
    catch(IException &e) { // catch exception from EvaluateForIntegration()
      throw IException(e,
                       e.errorType(),
                       "Simpsons4PointRule() - Unable to calculate the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using Simpson's 4 point rule",
                       _FILEINFO_);
    }
  }


  double NumericalApproximation::BoolesRule(const double a, const double b) {
    try {
      //n is the number of points used in the formula of the integration type chosen
      unsigned int n = 5;
      double result = 0;
      vector <double> f = EvaluateForIntegration(a, b, n);
      double h = f.back();
      f.pop_back();
      int ii;
      for(unsigned int i = 0; i < (f.size() - 1) / (n - 1); i++) {
        ii = (i + 1) * (n - 1);
        result += (7 * f[ii-4] + 32 * f[ii-3] + 12 * f[ii-2] + 32 * f[ii-1] + 7 * f[ii]) * h * 2 / 45;
      }
      return result;
    }
    catch(IException &e) { // catch exception from EvaluateForIntegration()
      throw IException(e,
                       e.errorType(),
                       "BoolesRule() - Unable to calculate the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using Boole's rule",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::RefineExtendedTrap(double a, double b, double s, unsigned int n) {
    // This method was derived from an algorithm in the text
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 4.2 by Flannery, Press, Teukolsky, and Vetterling
    try {
      if(n == 1) {
        double begin, end;
        if(GslInterpType(p_itype) || p_itype == NumericalApproximation::CubicNeighborhood) {
          // if a or b are outside the domain, return y-value of nearest endpoint
          begin = Evaluate(a, NumericalApproximation::NearestEndpoint);
          end = Evaluate(b, NumericalApproximation::NearestEndpoint);
        }
        else {
          // if a or b are outside the domain, return extrapolated y-value
          begin = Evaluate(a, NumericalApproximation::Extrapolate);
          end = Evaluate(b, NumericalApproximation::Extrapolate);
        }
        return (0.5 * (b - a) * (begin + end));
      }
      else {
        int it;
        double delta, tnm, x, sum;
        it = (int)(pow(2.0, (double)(n - 2)));
        tnm = it;
        delta = (b - a) / tnm; // spacing of the points to be added
        x = a + 0.5 * delta;
        sum = 0.0;
        for(int i = 0; i < it; i++) {
          if(GslInterpType(p_itype) || p_itype == NumericalApproximation::CubicNeighborhood) {
            sum = sum + Evaluate(x, NumericalApproximation::NearestEndpoint);
          }
          else {
            sum = sum + Evaluate(x, NumericalApproximation::Extrapolate);
          }
          x = x + delta;
        }
        return (0.5 * (s + (b - a) * sum / tnm));// return refined value of s
      }
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "RefineExtendedTrap() - Unable to calculate the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using the extended trapeziodal rule",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::RombergsMethod(double a, double b) {
    // This method was derived from an algorithm in the text
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 4.3 by Flannery, Press, Teukolsky, and Vetterling
    int maxits = 20;
    double dss = 0; // error estimate for
    double h[maxits+1]; // relative stepsizes for trap
    double trap[maxits+1]; // successive trapeziodal approximations
    double epsilon; // desired fractional accuracy
    double epsilon2;// desired fractional accuracy
    double ss; // result

    epsilon = 1.0e-4;
    epsilon2 = 1.0e-6;

    h[0] = 1.0;
    try {
      NumericalApproximation interp(NumericalApproximation::PolynomialNeville);
      for(int i = 0; i < maxits; i++) {
        // i will determine number of trapezoidal partitions of area
        // under curve for "integration" using refined trapezoidal rule
        trap[i] = RefineExtendedTrap(a, b, trap[i], i + 1); // validates data here
        if(i >= 4) {
          interp.AddData(5, &h[i-4], &trap[i-4]);
          // PolynomialNeville can extrapolate data outside of domain
          ss = interp.Evaluate(0.0, NumericalApproximation::Extrapolate);
          dss = interp.PolynomialNevilleErrorEstimate()[0];
          interp.Reset();
          // we work only until our necessary accuracy is achieved
          if(fabs(dss) <= epsilon * fabs(ss)) return ss;
          if(fabs(dss) <= epsilon2) return ss;
        }
        trap[i+1] = trap[i];
        h[i+1] = 0.25 * h[i];
        // This is a key step:  the factor is 0.25d0 even though
        // the stepsize is decreased by 0.5d0.  This makes the
        // extraplolation a polynomial in h-squared as allowed
        // by the equation from Numerical Recipes 4.2.1 pg.132,
        // not just a polynomial in h.
      }
    }
    catch(IException &e) { // catch error from RefineExtendedTrap, Constructor, Evaluate, PolynomialNevilleErrorEstimate
      throw IException(e,
                       e.errorType(),
                       "RombergsMethod() - Unable to calculate the integral on the interval (a,b) = ("
                       + IString(a) + ", " + IString(b)
                       + ") using Romberg's method",
                       _FILEINFO_);
    }
    throw IException(IException::Programmer,
                     "RombergsMethod() - Unable to calculate the integral using RombergsMethod() - Failed to converge in "
                     + IString(maxits) + " iterations",
                     _FILEINFO_);
  }

  void NumericalApproximation::Reset() {
    if(GslInterpType(p_itype)) GslDeallocation();
    p_clampedComputed = false;
    p_clampedEndptsSet = false;
    p_dataValidated = false;
    p_x.clear();
    p_y.clear();
    p_clampedSecondDerivs.clear();
    p_clampedDerivFirstPt = 0;
    p_clampedDerivLastPt = 0;
    p_polyNevError.clear();
    p_fprimeOfx.clear();
    return;
  }

  void NumericalApproximation::Reset(NumericalApproximation::InterpType itype) {
    try {
      Reset();
      SetInterpType(itype);
      return;
    }
    catch(IException &e) { // catch exception from SetInterpType()
      throw IException(e,
                       e.errorType(),
                       "Reset() - Unable to reset interpolation type",
                       _FILEINFO_);
    }
  }

  void NumericalApproximation::SetInterpType(NumericalApproximation::InterpType itype) {
    //  Validates the parameter
    if(GslInterpType(itype)) {
      try {
        GslFunctor(itype);
      }
      catch(IException &e) { // catch exception from GslFunctor()
        throw IException(e,
                         e.errorType(),
                         "SetInterpType() - Unable to set interpolation type",
                         _FILEINFO_);
      }
    }
    else if(itype > 9) {    // there are currently 9 interpolation types
      ReportException(IException::Programmer, "SetInterpType()",
                      "Invalid argument. Unknown interpolation type: "
                      + IString(NumericalApproximation::InterpType(itype)),
                      _FILEINFO_);
    }
    // p_x, p_y are kept and p_itype is replaced
    p_itype = itype;
    // reset state of class variables that are InterpType dependent  //??? should we keep some of this info?
    p_dataValidated = false;
    p_clampedComputed = false;
    p_clampedEndptsSet = false;
    p_clampedSecondDerivs.clear();
    p_clampedDerivFirstPt = 0;
    p_clampedDerivLastPt = 0;
    p_polyNevError.clear();
    p_fprimeOfx.clear();
  }

  void NumericalApproximation::Init(NumericalApproximation::InterpType itype) {
    if(p_interpFunctors.empty()) {
      p_interpFunctors.insert(make_pair(Linear, gsl_interp_linear));
      p_interpFunctors.insert(make_pair(Polynomial, gsl_interp_polynomial));
      p_interpFunctors.insert(make_pair(CubicNatural, gsl_interp_cspline));
      p_interpFunctors.insert(make_pair(CubicNatPeriodic, gsl_interp_cspline_periodic));
      p_interpFunctors.insert(make_pair(Akima, gsl_interp_akima));
      p_interpFunctors.insert(make_pair(AkimaPeriodic, gsl_interp_akima_periodic));
    }

    p_acc = 0;
    p_interp = 0;
    try {
      SetInterpType(itype);
    }
    catch(IException &e) { // catch exception from SetInterpType()
      throw IException(e,
                       e.errorType(),
                       "Init() - Unable to initialize NumericalApproximation object",
                       _FILEINFO_);
    }
    // Turn all GSL error handling off...repeatedly, every time this routine is
    // called.
    gsl_set_error_handler_off();
  }

  bool NumericalApproximation::GslInterpType(NumericalApproximation::InterpType itype) const {
    if(itype == NumericalApproximation::Linear)        return true;
    if(itype == NumericalApproximation::Polynomial)    return true;
    if(itype == NumericalApproximation::CubicNatural)  return true;
    if(itype == NumericalApproximation::CubicNatPeriodic) return true;
    if(itype == NumericalApproximation::Akima)         return true;
    if(itype == NumericalApproximation::AkimaPeriodic) return true;
    return false;
  }

  void NumericalApproximation::GslAllocation(unsigned int npoints) {
    GslDeallocation();
    //get pointer to accelerator object (iterator for interpolation lookups)
    p_acc = gsl_interp_accel_alloc();
    //get pointer to interpolation object of interp type given for npoints datapoints
    p_interp = gsl_spline_alloc(GslFunctor(p_itype), npoints);
    return;
  }

  void NumericalApproximation::GslDeallocation() {
    if(p_interp) gsl_spline_free(p_interp);
    if(p_acc) gsl_interp_accel_free(p_acc);
    p_acc = 0;
    p_interp = 0;
    return;
  }

  NumericalApproximation::InterpFunctor NumericalApproximation::GslFunctor(NumericalApproximation::InterpType itype)
  const  {
    FunctorConstIter fItr = p_interpFunctors.find(itype);
    if(fItr == p_interpFunctors.end()) {
      ReportException(IException::Programmer, "GslFunctor()",
                      "Invalid argument. Unable to find GSL interpolator with id = "
                      + IString(NumericalApproximation::InterpType(itype)),
                      _FILEINFO_);
    }
    return (fItr->second);
  }

  void NumericalApproximation::GslIntegrityCheck(int gsl_status, const char *src, int line)
   {
    if(gsl_status != GSL_SUCCESS) {
      if(gsl_status != GSL_EDOM) {
        ReportException(IException::Programmer, "GslIntegrityCheck(int,char,int)",
                        "GslIntegrityCheck(): GSL error occured: "
                        + string(gsl_strerror(gsl_status)), src, line);
      }
    }
    return;
  }

  void NumericalApproximation::ValidateDataSet() {
    if((int) Size() < MinPoints()) {
      ReportException(IException::Programmer, "ValidateDataSet()",
                      Name() + " interpolation requires a minimum of "
                      + IString(MinPoints()) + " data points - currently have "
                      + IString((int) Size()),
                      _FILEINFO_);
    }
    for(unsigned int i = 1; i < Size(); i++) {
      // Check for uniqueness -- this applies to all interpolation types
      if(p_x[i-1] == p_x[i]) {
        ReportException(IException::Programmer, "ValidateDataSet()",
                        "Invalid data set, x-values must be unique: \n\t\tp_x["
                        + IString((int) i - 1) + "] = " + IString(p_x[i-1])
                        + " = p_x[" + IString((int) i) + "]",
                        _FILEINFO_);
      }
      if(p_x[i-1] > p_x[i]) {
        // Verify that data set is in ascending order --
        // this does not apply to PolynomialNeville, which appears to get the same results with unsorted data
        if(p_itype != NumericalApproximation::PolynomialNeville) {
          ReportException(IException::Programmer, "ValidateDataSet()",
                          "Invalid data set, x-values must be in ascending order for "
                          + Name() + " interpolation: \n\t\tx["
                          + IString((int) i - 1) + "] = " + IString(p_x[i-1]) + " > x["
                          + IString((int) i) + "] = " + IString(p_x[i]),
                          _FILEINFO_);
        }
      }
    }
    if(p_itype == NumericalApproximation::CubicNatPeriodic) {
      if(p_y[0] != p_y[Size()-1]) {
        ReportException(IException::Programmer, "ValidateDataSet()",
                        "First and last points of the data set must have the same y-value for "
                        + Name() + "interpolation to prevent discontinuity at the boundary",
                        _FILEINFO_);
      }
    }
    p_dataValidated = true;
    return;
  }


  bool NumericalApproximation::InsideDomain(const double a) {
    try {
      if(a + DBL_EPSILON < DomainMinimum()) {
        return false;
      }
      if(a - DBL_EPSILON > DomainMaximum()) {
        return false;
      }
    }
    catch(IException &e) { // catch exception from DomainMinimum(), DomainMaximum()
      throw IException(e,
                       e.errorType(),
                       "InsideDomain() - Unable to compute domain boundaries",
                       _FILEINFO_);
    }
    return true;
  }

  bool NumericalApproximation::GslComputed() const {
    if(GslInterpType(p_itype)) return ((p_interp) && (p_acc));
    else
      throw IException(IException::Programmer,
                       "GslComputed() - Method only valid for GSL interpolation types, may not be used for "
                       + Name() + " interpolation",
                       _FILEINFO_);
  }

  void NumericalApproximation::ComputeGsl() {
    try {
      if(GslComputed()) return;
      if(!p_dataValidated) ValidateDataSet();
      GslAllocation(Size());
      GslIntegrityCheck(gsl_spline_init(p_interp, &p_x[0], &p_y[0], Size()), _FILEINFO_);
      return;
    }
    catch(IException &e) { // catch exception from ValidateDataSet(), GslIntegrityCheck()
      throw IException(e,
                       e.errorType(),
                       "ComputeGsl() - Unable to compute GSL interpolation",
                       _FILEINFO_);
    }
  }

  void NumericalApproximation::ComputeCubicClamped() {
    // This method was derived from an algorithm in the text
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 3.3 by Flannery, Press, Teukolsky, and Vetterling

    if(!p_dataValidated) {
      try {
        ValidateDataSet();
      }
      catch(IException &e) {  // catch exception from ValidateDataSet()
        throw IException(e,
                         e.errorType(),
                         "ComputeCubicClamped() - Unable to compute cubic clamped interpolation",
                         _FILEINFO_);
      }
    }
    if(!p_clampedEndptsSet) {
      ReportException(IException::Programmer, "ComputeCubicClamped()",
                      "Must set endpoint derivative values after adding data in order to compute cubic spline with clamped boundary conditions",
                      _FILEINFO_);
    }
    int n = Size();
    p_clampedSecondDerivs.resize(n);
    double u[n];
    double p, sig, qn, un;

    if(p_clampedDerivFirstPt > 0.99e30) {
      p_clampedSecondDerivs[0] = 0.0;//natural boundary conditions are used if deriv of first value is greater than 10^30
      u[0] = 0.0;
    }
    else {
      p_clampedSecondDerivs[0] = -0.5;// clamped conditions are used
      u[0] = (3.0 / (p_x[1] - p_x[0])) * ((p_y[1] - p_y[0]) / (p_x[1] - p_x[0]) - p_clampedDerivFirstPt);
    }
    for(int i = 1; i < n - 1; i++) { // decomposition loop of the tridiagonal algorithm
      sig = (p_x[i] - p_x[i-1]) / (p_x[i+1] - p_x[i-1]);
      p = sig * p_clampedSecondDerivs[i-1] + 2.0;
      p_clampedSecondDerivs[i] = (sig - 1.0) / p;
      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]) /
                     (p_x[i] - p_x[i-1])) / (p_x[i+1] - p_x[i-1]) - sig * u[i-1]) / p;
    }
    if(p_clampedDerivLastPt > 0.99e30) { // upper boundary is natural
      qn = 0.0;
      un = 0.0;
    }
    else {// upper boundary is clamped
      qn = 0.5;
      un = (3.0 / (p_x[n-1] - p_x[n-2])) * (p_clampedDerivLastPt - (p_y[n-1] - p_y[n-2]) /
                                            (p_x[n-1] - p_x[n-2]));
    }
    p_clampedSecondDerivs[n-1] = (un - qn * u[n-2]) / (qn * p_clampedSecondDerivs[n-2] + 1.0);
    for(int i = n - 2; i >= 0; i--) { // backsubstitution loop of the tridiagonal algorithm
      p_clampedSecondDerivs[i] = p_clampedSecondDerivs[i] * p_clampedSecondDerivs[i+1] + u[i];
    }
    p_clampedComputed = true;
    return;
  }

  double NumericalApproximation::ValueToExtrapolate(const double a, const ExtrapType &etype) {
    if(etype == NumericalApproximation::ThrowError) {
      ReportException(IException::Programmer, "Evaluate()",
                      "Invalid argument. Value entered, a = "
                      + IString(a) + ", is outside of domain = ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    if(etype == NumericalApproximation::NearestEndpoint) {
      if(a + DBL_EPSILON < DomainMinimum()) {
        return DomainMinimum();
      }
      else {
        return DomainMaximum(); // (a > DomainMaximum())
      }
    }
    else {  // gsl interpolations and CubicNeighborhood cannot extrapolate
      if(GslInterpType(p_itype)
          || p_itype == NumericalApproximation::CubicNeighborhood) {
        ReportException(IException::Programmer, "Evaluate()",
                        "Invalid argument. Cannot extrapolate for type "
                        + Name() + ", must choose to throw error or return nearest neighbor",
                        _FILEINFO_);
      }
      return a;
    }
  }

  double NumericalApproximation::EvaluateCubicNeighborhood(const double a) {
    try {
      int s = 0, s_old, s0;
      vector <double> x0(4), y0(4);
      NumericalApproximation spline(NumericalApproximation::CubicNatural);
      for(unsigned int n = 0; n < Size(); n++) {
        if(p_x[n] < a) s = n;
        else break;
      }
      if(s < 1) s = 1;
      else if(s > (int) Size() - 3) s = Size() - 3;
      s_old = -1;
      s0 = s - 1;
      if(s_old != s0) {
        for(int n = 0; n < 4; n++) {
          x0[n] = p_x[n+s0];
          y0[n] = p_y[n+s0];
          spline.AddData(x0[n], y0[n]);
        }
        s_old = s0;
      }
      // use nearest endpoint extrapolation method for neighborhood spline
      // since CubicNatural can do no other
      return spline.Evaluate(a, NumericalApproximation::NearestEndpoint);
    }
    catch(IException &e) { // catch exception from Constructor, ComputeGsl(), Evaluate()
      throw IException(e,
                       e.errorType(),
                       "EvaluateCubicNeighborhood() - Unable to evaluate cubic neighborhood interpolation at a = "
                       + IString(a),
                       _FILEINFO_);
    }
  }

  vector <double> NumericalApproximation::EvaluateCubicNeighborhood(const vector <double> &a,
      const NumericalApproximation::ExtrapType &etype) {
    vector <double> result(a.size());
    int s_old, s0;
    vector <double> x0(4), y0(4);
    vector <int> s;
    s.clear();
    s.resize(a.size());
    for(unsigned int i = 0; i < a.size(); i++) {
      for(unsigned int n = 0; n < Size(); n++) {
        if(p_x[n] < a[i]) {
          s[i] = n;
        }
      }
      if(s[i] < 1) {
        s[i] = 1;
      }
      if(s[i] > ((int) Size()) - 3) {
        s[i] = Size() - 3;
      }
    }
    s_old = -1;
    try {
      NumericalApproximation spline(NumericalApproximation::CubicNatural);
      for(unsigned int i = 0; i < a.size(); i++) {
        s0 = s[i] - 1;
        if(s_old != s0) {
          spline.Reset();
          for(int n = 0; n < 4; n++) {
            x0[n] = p_x[n+s0];
            y0[n] = p_y[n+s0];
            spline.AddData(x0[n], y0[n]);
          }
          s_old = s0;
        }
        double a0;
        // checks whether this value is in domain of the main spline
        if(InsideDomain(a[i]))
          a0 = a[i];
        else a0 = ValueToExtrapolate(a[i], etype);
        // since neighborhood spline is CubicNatural, the only extrapolation possible is NearestEndpoint
        result[i] = spline.Evaluate(a0, NumericalApproximation::NearestEndpoint);
      }
      return result;
    }
    catch(IException &e) { // catch exception from Constructor, ComputeGsl(), Evaluate()
      throw IException(e,
                       e.errorType(),
                       "EvaluateCubicNeighborhood() - Unable to evaluate the function at the given vector of points using cubic neighborhood interpolation",
                       _FILEINFO_);
    }
  }

  double NumericalApproximation::EvaluateCubicClamped(double a) {
    // This method was derived from an algorithm in the text
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 3.3 by Flannery, Press, Teukolsky, and Vetterling
    int n = Size();
    double result = 0;

    int k;
    int k_Lo;
    int k_Hi;
    double h;
    double A;
    double B;

    k_Lo = 0;
    k_Hi = n - 1;
    while(k_Hi - k_Lo > 1) {
      k = (k_Hi + k_Lo) / 2;
      if(p_x[k] > a) {
        k_Hi = k;
      }
      else {
        k_Lo = k;
      }
    }

    h = p_x[k_Hi] - p_x[k_Lo];
    A = (p_x[k_Hi] - a) / h;
    B = (a - p_x[k_Lo]) / h;
    result = A * p_y[k_Lo] + B * p_y[k_Hi] + ((pow(A, 3.0) - A) *
             p_clampedSecondDerivs[k_Lo] + (pow(B, 3.0) - B) * p_clampedSecondDerivs[k_Hi]) * pow(h, 2.0) / 6.0;
    return result;
  }

  double NumericalApproximation::EvaluateCubicHermite(const double a) {
    //  algorithm was found at en.wikipedia.org/wiki/Cubic_Hermite_spline
    //  it seems to produce same answers, as the NumericalAnalysis book

    if(p_fprimeOfx.size() != Size()) {
      ReportException(IException::User, "EvaluateCubicHermite()",
                      "Invalid arguments. The size of the first derivative vector does not match the number of (x,y) data points.",
                      _FILEINFO_);
    }
    // find the interval in which "a" exists
    int lowerIndex = FindIntervalLowerIndex(a);

    // we know that "a" is within the domain since this is verified in
    // Evaluate() before this method is called, thus n <= Size()
    if(a == p_x[lowerIndex]) {
      return p_y[lowerIndex];
    }
    if(a == p_x[lowerIndex+1]) {
      return p_y[lowerIndex+1];
    }

    double x0, x1, y0, y1, m0, m1;
    // a is contained within the interval (x0,x1)
    x0 = p_x[lowerIndex];
    x1 = p_x[lowerIndex+1];
    // the corresponding known y-values for x0 and x1
    y0 = p_y[lowerIndex];
    y1 = p_y[lowerIndex+1];
    // the corresponding known tangents (slopes) at (x0,y0) and (x1,y1)
    m0 = p_fprimeOfx[lowerIndex];
    m1 = p_fprimeOfx[lowerIndex+1];


    //  following algorithm found at en.wikipedia.org/wiki/Cubic_Hermite_spline
    //  seems to produce same answers, is it faster?

    double h, t;
    h = x1 - x0;
    t = (a - x0) / h;
    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;
  }

  int NumericalApproximation::FindIntervalLowerIndex(const double a) {
    if(InsideDomain(a)) {
      // find the interval in which "a" exists
      std::vector<double>::iterator pos;
      // find position in vector that is greater than or equal to "a"
      pos = upper_bound(p_x.begin(), p_x.end(), a);
      int upperIndex = 0;
      if(pos != p_x.end()) {
        upperIndex = distance(p_x.begin(), pos);
      }
      else {
        upperIndex = Size() - 1;
      }
      return upperIndex - 1;
    }
    else if((a + DBL_EPSILON) < DomainMinimum()) {
      return 0;
    }
    else {
      return Size() - 2;
    }
  }


  double NumericalApproximation::EvaluatePolynomialNeville(const double a) {
    // This method was derived from an algorithm in the text
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 3.1 by Flannery, Press, Teukolsky, and Vetterling
    int n = Size();
    double y;

    int ns;
    double den, dif, dift, c[n], d[n], ho, hp, w;
    double *err = 0;
    ns = 1;
    dif = fabs(a - p_x[0]);
    for(int i = 0; i < n; i++) { // Get the index ns of the closest table entry
      dift = fabs(a - p_x[i]);
      if(dift < dif) {
        ns = i + 1;
        dif = dift;
      }
      c[i] = p_y[i];  // initialize c and d
      d[i] = p_y[i];
    }
    ns = ns - 1;
    y = p_y[ns]; // initial approximation for y
    for(int m = 1; m < n; m++) {
      for(int i = 1; i <= n - m; i++) { // loop over c and d and update them
        ho = p_x[i-1] - a;
        hp = p_x[i+m-1] - a;
        w = c[i] - d[i-1];
        den = ho - hp;
        den = w / den;
        d[i-1] = hp * den; // update c and d
        c[i-1] = ho * den;
      }
      if(2 * ns < n - m) { // After each column in the tableau is completed, we decide
        err = &c[ns];   // which correction, c or d, we want to add to our accumulating
      }                 // value of y, i.e., which path to take through the tableau|
      else {            // forking up or down. We do this in such a way as to take the
        ns = ns - 1;    // most "straight line" route through the tableau to its apex,
        err = &d[ns];   // updating ns accordingly to keep track of where we are.  This
      }                 // route keeps the partial approximations centered (insofar as possible)
      y = y + *err;     // on the target x. The last err added is thus the error indication.
    }
    p_polyNevError.push_back(*err);
    return y;
  }

  vector <double> NumericalApproximation::EvaluateForIntegration(const double a, const double b, const unsigned int n) {
    if(a > b) {
      ReportException(IException::Programmer, "EvaluatedForIntegration()",
                      "Invalid interval entered: [a,b] = ["
                      + IString(a) + ", " + IString(b) + "]",
                      _FILEINFO_);
    }
    if(!InsideDomain(a) || !InsideDomain(b)) {
      ReportException(IException::Programmer, "EvaluateForIntegration()",
                      "Invalid arguments. Interval entered ["
                      + IString(a) + ", " + IString(b)
                      + "] is not contained within domain ["
                      + IString(DomainMinimum()) + ", "
                      + IString(DomainMaximum()) + "]",
                      _FILEINFO_);
    }
    vector <double> f;
    //need total number of segments to be divisible by n-1
    // This way we can use the formula for each interval: 0 to n, n to 2n, 2n to 3n, etc.
    // Notice interval endpoints will overlap for these composite formulas
    int xSegments = Size() - 1;
    while(xSegments % (n - 1) != 0) {
      xSegments++;
    }
    // x must be sorted and unique
    double xMin = a;
    double xMax = b;
    //Uniform step size.
    double h = (xMax - xMin) / xSegments;
    //Compute the interpolates using spline.
    try {
      for(double i = 0; i < (xSegments + 1); i++) {
        double xi = h * i + xMin;
        f.push_back(Evaluate(xi));  // validate data set here, allow default ThrowError extrap type
      }
      f.push_back(h);
      return f;
    }
    catch(IException &e) { // catch exception from Evaluate()
      throw IException(e,
                       e.errorType(),
                       "EvaluateForIntegration() - Unable to evaluate the data set for integration",
                       _FILEINFO_);
    }
  }// for integration using composite of spline (creates evenly spaced points)

  void NumericalApproximation::ReportException(IException::ErrorType type, const string &methodName, const string &message,
      const char *filesrc, int lineno)
  const  {
    string msg = methodName + " - " + message;
    throw IException(type, msg.c_str(), filesrc, lineno);
    return;
  }

}