SplineInterpolatorEinspline.h

Go to the documentation of this file.

#ifndef _utl_SplineInterpolatorEinspline_h_
#define _utl_SplineInterpolatorEinspline_h_

#include <iostream>

#include <boost/array.hpp>
#include <boost/multi_array.hpp>

#include <utl/AugerException.h>

/* fix annoying compiler warnings
 * because of useless macro constants
 * set by both einspline and CLHEP */
#undef PACKAGE
#undef PACKAGE_BUGREPORT
#undef PACKAGE_STRING
#undef PACKAGE_NAME
#undef PACKAGE_TARNAME
#undef PACKAGE_VERSION
#undef VERSION
#include <einspline/bspline_base.h>
#include <einspline/bspline_structs.h>
#include <einspline/bspline_create.h>
#include <einspline/bspline_eval_std_d.h>

#include <einspline/multi_bspline_structs.h>
#include <einspline/multi_bspline_create.h>
#include <einspline/multi_bspline_eval_d.h>
#undef PACKAGE
#undef PACKAGE_BUGREPORT
#undef PACKAGE_STRING
#undef PACKAGE_NAME
#undef PACKAGE_TARNAME
#undef PACKAGE_VERSION
#undef VERSION

#include <vector>


namespace utl {

  namespace Spline {

    enum Type {
      eNatural  = NATURAL,         /* minimize total curvature          */
      ePeriodic = PERIODIC,        /* for a polar coordinate            */
      eFirstDerivative  = DERIV1,  /* set 1st derivatives at endpoints  */
      eSecondDerivative = DERIV2,  /* set 2nd derivatives at endpoints  */
    };


    class BoundaryCondition : public BCtype_d {
    public:
      BoundaryCondition(const Type type = eNatural,
                        const double vLeft = 0,
                        const double vRight = 0)
      {
        lCode = static_cast<const bc_code>(type);
        lVal = vLeft;
        rCode = static_cast<const bc_code>(type);
        rVal = vRight;
      }
    };


    class Interpolator1D {
    public:
      Interpolator1D() { }

      Interpolator1D(const Interpolator1D& other) { *this = other; }

      Interpolator1D(const double xStart,
                     const double xStop, /* inclusive! */
                     const std::vector<double>& fValue,
                     const BoundaryCondition& bcX = BoundaryCondition())
      {
        Ugrid grid;
        grid.start = xStart;
        grid.end = xStop;
        grid.num = fValue.size();
        fDataPtr = create_UBspline_1d_d(grid, bcX, const_cast<double*>(&fValue[0]));
      }

      ~Interpolator1D() { destroy_Bspline(fDataPtr); }

      Interpolator1D&
      operator=(const Interpolator1D& other)
      {
        if (this != &other) {
          if (fDataPtr)
            destroy_Bspline(fDataPtr);
          fDataPtr = new UBspline_1d_d;
          *fDataPtr = *other.fDataPtr;
          const int n = fDataPtr->x_grid.num + (fDataPtr->xBC.lCode == PERIODIC ? 3 : 2);
          fDataPtr->coefs = new double[n];
          for (int i = 0; i < n; ++i)
            fDataPtr->coefs[i] = other.fDataPtr->coefs[i];
        }
        return *this;
      }

      double
      operator()(const double x)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        if (x < fDataPtr->x_grid.start || x > fDataPtr->x_grid.end)
          throw utl::OutOfBoundException("argument is outside of interpolation range");

        double result;
        eval_UBspline_1d_d(fDataPtr, x, &result);
        return result;
      }

      double
      GetStart()
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        return fDataPtr->x_grid.start;
      }

      double
      GetStop()
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        return fDataPtr->x_grid.end;
      }

    protected:
      UBspline_1d_d* fDataPtr = nullptr;
    };


    class Interpolator2D {
    public:
      Interpolator2D() { }

      Interpolator2D(const Interpolator2D& other) { *this = other; }

      Interpolator2D(const double xStart,
                     const double xStop,
                     const double yStart,
                     const double yStop,
                     const boost::multi_array<double,2>& values,
                     const BoundaryCondition bcX = BoundaryCondition(),
                     const BoundaryCondition bcY = BoundaryCondition())
      {
        Ugrid grid[2];
        grid[0].num = values.shape()[0];
        grid[0].start = xStart;
        grid[0].end = xStop;
        grid[1].num = values.shape()[1];
        grid[1].start = yStart;
        grid[1].end = yStop;
        fDataPtr = create_UBspline_2d_d(grid[0], grid[1], bcX, bcY, const_cast<double*>(values.data()));
      }

      ~Interpolator2D() { destroy_Bspline(fDataPtr); }

      Interpolator2D&
      operator=(const Interpolator2D& other)
      {
        if (this != &other) {
          if (fDataPtr)
            destroy_Bspline(fDataPtr);
          fDataPtr = new UBspline_2d_d;
          *fDataPtr = *other.fDataPtr;
          const int n =
            (fDataPtr->x_grid.num + (fDataPtr->xBC.lCode == PERIODIC ? 3 : 2)) *
            (fDataPtr->y_grid.num + (fDataPtr->yBC.lCode == PERIODIC ? 3 : 2));
          fDataPtr->coefs = new double[n];
          for (int i = 0; i < n; ++i)
            fDataPtr->coefs[i] = other.fDataPtr->coefs[i];
        }
        return *this;
      }

      double
      operator()(const double x, const double y)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        if (x < fDataPtr->x_grid.start || x > fDataPtr->x_grid.end)
          throw utl::OutOfBoundException("x argument is outside of interpolation range");
        if (y < fDataPtr->y_grid.start || y > fDataPtr->y_grid.end)
          throw utl::OutOfBoundException("y argument is outside of interpolation range");

        double result;
        eval_UBspline_2d_d(fDataPtr, x, y, &result);
        return result;
      }

      double
      GetStart(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.start;
        case 1:
          return fDataPtr->y_grid.start;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

      double
      GetStop(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.end;
        case 1:
          return fDataPtr->y_grid.end;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

    protected:
      UBspline_2d_d* fDataPtr = nullptr;
    };


    class Interpolator3D {
    public:
      Interpolator3D() { }

      Interpolator3D(const Interpolator3D& other) { *this = other; }

      Interpolator3D(const double xStart,
                     const double xStop, /* inclusive! */
                     const double yStart,
                     const double yStop, /* inclusive! */
                     const double zStart,
                     const double zStop, /* inclusive! */
                     const boost::multi_array<double,3>& values,
                     const BoundaryCondition bcX = BoundaryCondition(),
                     const BoundaryCondition bcY = BoundaryCondition(),
                     const BoundaryCondition bcZ = BoundaryCondition())
      {
        Ugrid grid[3];
        grid[0].num = values.shape()[0];
        grid[0].start = xStart;
        grid[0].end = xStop;
        grid[1].num = values.shape()[1];
        grid[1].start = yStart;
        grid[1].end = yStop;
        grid[2].num = values.shape()[2];
        grid[2].start = zStart;
        grid[2].end = zStop;

        fDataPtr = create_UBspline_3d_d(grid[0], grid[1], grid[2], bcX, bcY, bcZ, const_cast<double*>(values.data()));
      }

      ~Interpolator3D() { destroy_Bspline(fDataPtr); }

      Interpolator3D&
      operator=(const Interpolator3D& other)
      {
        if (this != &other) {
          if (fDataPtr)
            destroy_Bspline(fDataPtr);
          fDataPtr = new UBspline_3d_d;
          *fDataPtr = *other.fDataPtr;
          const int n =
            (fDataPtr->x_grid.num + (fDataPtr->xBC.lCode == PERIODIC ? 3 : 2)) *
            (fDataPtr->y_grid.num + (fDataPtr->yBC.lCode == PERIODIC ? 3 : 2)) *
            (fDataPtr->z_grid.num + (fDataPtr->zBC.lCode == PERIODIC ? 3 : 2));
          fDataPtr->coefs = new double[n];
          for (int i = 0; i < n; ++i)
            fDataPtr->coefs[i] = other.fDataPtr->coefs[i];
        }
        return *this;
      }

      double
      operator()(const double x, const double y, const double z)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        if (x < fDataPtr->x_grid.start || x > fDataPtr->x_grid.end)
          throw utl::OutOfBoundException("x argument is outside of interpolation range");
        if (y < fDataPtr->y_grid.start || y > fDataPtr->y_grid.end)
          throw utl::OutOfBoundException("y argument is outside of interpolation range");
        if (z < fDataPtr->z_grid.start || z > fDataPtr->z_grid.end)
          throw utl::OutOfBoundException("z argument is outside of interpolation range");

        double result;
        eval_UBspline_3d_d(fDataPtr, x, y, z, &result);
        return result;
      }

      double
      GetStart(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.start;
        case 1:
          return fDataPtr->y_grid.start;
        case 2:
          return fDataPtr->z_grid.start;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

      double
      GetStop(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.end;
        case 1:
          return fDataPtr->y_grid.end;
        case 2:
          return fDataPtr->z_grid.end;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

    protected:
      UBspline_3d_d* fDataPtr = nullptr;
    };


    class VectorInterpolator1D {
    public:
      VectorInterpolator1D() { }

      VectorInterpolator1D(const VectorInterpolator1D& other) { *this = other; }

      VectorInterpolator1D(const double xStart,
                           const double xStop, /* inclusive! */
                           const boost::multi_array<double,2>& fValue, // first index is function vector index!
                           const BoundaryCondition& bcX = BoundaryCondition())
      {
        Ugrid grid;
        grid.start = xStart;
        grid.end = xStop;
        grid.num = fValue.shape()[1];

        fDataPtr = create_multi_UBspline_1d_d(grid, bcX, fValue.shape()[0]);
        for (size_t i = 0; i < fValue.shape()[0]; ++i)
          set_multi_UBspline_1d_d(fDataPtr, i, const_cast<double*>(fValue.data() + i*grid.num));
      }

      ~VectorInterpolator1D() { destroy_Bspline(fDataPtr); }

      VectorInterpolator1D&
      operator=(const VectorInterpolator1D& other)
      {
        if (this != &other) {
          if (fDataPtr)
            destroy_Bspline(fDataPtr);
          fDataPtr = new multi_UBspline_1d_d;
          *fDataPtr = *other.fDataPtr;
          const int n = fDataPtr->num_splines *
            (fDataPtr->x_grid.num + (fDataPtr->xBC.lCode == PERIODIC ? 3 : 2));
          fDataPtr->coefs = new double[n];
          for (int i = 0; i < n; ++i)
            fDataPtr->coefs[i] = other.fDataPtr->coefs[i];
          }
        return *this;
      }

      void
      operator()(std::vector<double>& result, const double x)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        if (x < fDataPtr->x_grid.start || x > fDataPtr->x_grid.end)
          throw utl::OutOfBoundException("argument is outside of interpolation range");
        if (result.size() != size_t(fDataPtr->num_splines))
          result.resize(fDataPtr->num_splines);
        eval_multi_UBspline_1d_d(fDataPtr, x, &result[0]);
      }

      double
      GetStart()
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        return fDataPtr->x_grid.start;
      }

      double
      GetStop()
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        return fDataPtr->x_grid.end;
      }

    protected:
      multi_UBspline_1d_d* fDataPtr = nullptr;
    };


    class VectorInterpolator2D {
    public:
      VectorInterpolator2D() { }

      VectorInterpolator2D(const VectorInterpolator2D& other) { *this = other; }

      VectorInterpolator2D(const double xStart,
                           const double xStop, /* inclusive! */
                           const double yStart,
                           const double yStop, /* inclusive! */
                           const boost::multi_array<double,3>& fValue, // first index is function vector index!
                           const BoundaryCondition& bcX = BoundaryCondition(),
                           const BoundaryCondition& bcY = BoundaryCondition())
      {
        Ugrid grid[2];
        grid[0].start = xStart;
        grid[0].end = xStop;
        grid[0].num = fValue.shape()[1];
        grid[1].start = yStart;
        grid[1].end = yStop;
        grid[1].num = fValue.shape()[2];

        fDataPtr = create_multi_UBspline_2d_d(grid[0], grid[1], bcX, bcY, fValue.shape()[0]);
        for (size_t i = 0; i < fValue.shape()[0]; ++i)
          set_multi_UBspline_2d_d(fDataPtr, i, const_cast<double*>(fValue.data() + i*grid[0].num*grid[1].num));
      }

      ~VectorInterpolator2D() { destroy_Bspline(fDataPtr); }

      VectorInterpolator2D&
      operator=(const VectorInterpolator2D& other)
      {
        if (this != &other) {
          if (fDataPtr)
            destroy_Bspline(fDataPtr);
          fDataPtr = new multi_UBspline_2d_d;
          *fDataPtr = *other.fDataPtr;
          const int n = fDataPtr->num_splines *
            (fDataPtr->x_grid.num + (fDataPtr->xBC.lCode == PERIODIC ? 3 : 2)) *
            (fDataPtr->y_grid.num + (fDataPtr->yBC.lCode == PERIODIC ? 3 : 2));
          fDataPtr->coefs = new double[n];
          for (int i = 0; i < n; ++i)
            fDataPtr->coefs[i] = other.fDataPtr->coefs[i];
        }
        return *this;
      }

      void
      operator()(std::vector<double>& result, const double x, const double y)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        if (x < fDataPtr->x_grid.start || x > fDataPtr->x_grid.end)
          throw utl::OutOfBoundException("x argument is outside of interpolation range");
        if (y < fDataPtr->y_grid.start || y > fDataPtr->y_grid.end)
          throw utl::OutOfBoundException("y argument is outside of interpolation range");
        if (result.size() != size_t(fDataPtr->num_splines))
          result.resize(fDataPtr->num_splines);
        eval_multi_UBspline_2d_d(fDataPtr, x, y, &result[0]);
      }

      double
      GetStart(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.start;
        case 1:
          return fDataPtr->y_grid.start;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

      double
      GetStop(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.end;
        case 1:
          return fDataPtr->y_grid.end;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

    protected:
      multi_UBspline_2d_d* fDataPtr = nullptr;
    };


    class VectorInterpolator3D {
    public:
      VectorInterpolator3D() { }

      VectorInterpolator3D(const VectorInterpolator3D& other) { *this = other; }

      VectorInterpolator3D(const double xStart,
                           const double xStop, /* inclusive! */
                           const double yStart,
                           const double yStop, /* inclusive! */
                           const double zStart,
                           const double zStop, /* inclusive! */
                           const boost::multi_array<double,4>& fValue, // first index is function vector index!
                           const BoundaryCondition& bcX = BoundaryCondition(),
                           const BoundaryCondition& bcY = BoundaryCondition(),
                           const BoundaryCondition& bcZ = BoundaryCondition())
      {
        Ugrid grid[3];
        grid[0].start = xStart;
        grid[0].end = xStop;
        grid[0].num = fValue.shape()[1];
        grid[1].start = yStart;
        grid[1].end = yStop;
        grid[1].num = fValue.shape()[2];
        grid[2].start = zStart;
        grid[2].end = zStop;
        grid[2].num = fValue.shape()[3];

        fDataPtr = create_multi_UBspline_3d_d(grid[0], grid[1], grid[2], bcX, bcY, bcZ, fValue.shape()[0]);
        for (size_t i = 0; i < fValue.shape()[0]; ++i)
          set_multi_UBspline_3d_d(fDataPtr, i, const_cast<double*>(fValue.data() + i*grid[0].num*grid[1].num*grid[2].num));
      }

      ~VectorInterpolator3D() { destroy_Bspline(fDataPtr); }

      VectorInterpolator3D&
      operator=(const VectorInterpolator3D& other)
      {
        if (this != &other) {
          if (fDataPtr)
            destroy_Bspline(fDataPtr);
          fDataPtr = new multi_UBspline_3d_d;
          *fDataPtr = *other.fDataPtr;
          const int n = fDataPtr->num_splines *
            (fDataPtr->x_grid.num + (fDataPtr->xBC.lCode == PERIODIC ? 3 : 2)) *
            (fDataPtr->y_grid.num + (fDataPtr->yBC.lCode == PERIODIC ? 3 : 2)) *
            (fDataPtr->z_grid.num + (fDataPtr->zBC.lCode == PERIODIC ? 3 : 2));
          fDataPtr->coefs = new double[n];
          for (int i = 0; i < n; ++i)
            fDataPtr->coefs[i] = other.fDataPtr->coefs[i];
        }
        return *this;
      }

      void
      operator()(std::vector<double>& result, const double x, const double y, const double z)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        if (x < fDataPtr->x_grid.start || x > fDataPtr->x_grid.end)
          throw utl::OutOfBoundException("x argument is outside of interpolation range");
        if (y < fDataPtr->y_grid.start || y > fDataPtr->y_grid.end)
          throw utl::OutOfBoundException("y argument is outside of interpolation range");
        if (z < fDataPtr->z_grid.start || z > fDataPtr->z_grid.end)
          throw utl::OutOfBoundException("z argument is outside of interpolation range");
        if (result.size() != size_t(fDataPtr->num_splines))
          result.resize(fDataPtr->num_splines);
        eval_multi_UBspline_3d_d(fDataPtr, x, y, z, &result[0]);
      }

      double
      GetStart(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.start;
        case 1:
          return fDataPtr->y_grid.start;
        case 2:
          return fDataPtr->z_grid.start;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

      double
      GetStop(const unsigned char dimension)
        const
      {
        if (!fDataPtr)
          throw utl::InvalidConfigurationException("no data for spline interpolator");
        switch (dimension) {
        case 0:
          return fDataPtr->x_grid.end;
        case 1:
          return fDataPtr->y_grid.end;
        case 2:
          return fDataPtr->z_grid.end;
        default:
          throw utl::DoesNotComputeException("dimension is too large");
        }
        return 0;
      }

    protected:
      multi_UBspline_3d_d* fDataPtr = nullptr;
    };

  }

}


#endif