#ifndef _SplineSolver_Spline_h_
#define _SplineSolver_Spline_h_

#include <cmath>
#include <vector>
#include <boost/multi_array.hpp>

#include <utl/EquationSolver.h>

namespace Spline {

  typedef unsigned char dim_t;

  enum BoundaryConditionType {
    eNatural          = 0, /* minimize total curvature          */
    ePeriodic         = 1, /* for a polar coordinate            */
    eFirstDerivative  = 2, /* set 1st derivatives at endpoints  */
    eSecondDerivative = 3  /* set 2nd derivatives at endpoints  */
  };

  struct BoundaryCondition
  {
    BoundaryCondition(const BoundaryConditionType type = eNatural,
                      const double vLeft = 0.0,
                      const double vRight = 0.0) :
      fType(type), fLeft(vLeft), fRight(vRight)
    {
    }

    inline
    bool
    operator==(const BoundaryConditionType type)
      const
    {
      return fType == type;
    }

    BoundaryConditionType fType;
    double fLeft;
    double fRight;
  };

  /// Computes B-spline coefficients
  template<dim_t ADimension, typename AKnotVector, typename ABasisFunction>
  class Solver
  {
  public:

    /// Prepare inversion matrix
    void
    Configure(const dim_t idim, const AKnotVector& xs, const std::vector<ABasisFunction>& bs, const BoundaryCondition& bc)
    {
      fBCs[idim] = bc;
      fSizes[idim] = bs.size();
      if (fBuffer.size() < bs.size())
        fBuffer.resize(bs.size());

      const size_t n = bs.size();
      fInv[idim].resize(boost::extents[n][n]);
      for (size_t i=0;i<n;++i)
        for (size_t j=0;j<n;++j)
          fInv[idim][i][j] = 0;

      for (size_t i = 0; i < n-2; ++i)
        for (size_t k = i; k < i+3; ++k)
          fInv[idim][i][k] = bs[k](xs[i]);

      if (bc == eFirstDerivative)
      {
        // set first derivative at x0 and xn
        for (size_t der = 0; der < 2; ++der)
          for (size_t k = 0; k < 3; ++k)
        {
          fInv[idim][n-2][k] = bs[k](xs[0],1);
          fInv[idim][n-1][n-3+k] = bs[n-3+k](xs[n-3],1);
        }
      }

      if (bc == eNatural || bc == eSecondDerivative)
      {
        // second derivative zero at x0 and xn
        for (size_t k = 0; k < 3; ++k)
        {
          fInv[idim][n-2][k] = bs[k](xs[0],2);
          fInv[idim][n-1][n-3+k] = bs[n-3+k](xs[n-3],2);
        }
      }

      if (bc == ePeriodic)
      {
        // first and second derivative at x0 and xn equal
        for (size_t der = 0; der < 2; ++der)
          for (size_t k = 0; k < 3; ++k)
        {
          fInv[idim][n-2+der][k] = bs[k](xs[0],1+der);
          fInv[idim][n-2+der][n-3+k] = -bs[n-3+k](xs[n-3],1+der);
        }
      }

      utl::InvertMatrix(n, fInv[idim]);
    }

    void
    operator()(double* cs)
      const
    {
      size_t strides[ADimension];
      strides[0] = 1;
      for (short idim = 0; idim < short(ADimension-1); ++idim)
        strides[idim+1] = strides[idim]*fSizes[idim];

      for (dim_t idim = 0; idim < ADimension; ++idim)
      {
        size_t strides2[ADimension+1];
        strides2[0] = 1;
        for (dim_t jdim = 0; jdim < ADimension; ++jdim)
        {
          const size_t n = jdim == idim ? 1 :
            (jdim < idim ? fSizes[jdim] : fSizes[jdim]-2);
          strides2[jdim+1] = strides2[jdim]*n;
        }

        for (size_t i = 0; i < strides2[ADimension]; ++i)
        {
          size_t j = 0;
          for (dim_t jdim = 0; jdim < ADimension; ++jdim)
          {
            if (idim != jdim)
              j += ((i % strides2[jdim+1])/strides2[jdim]) * strides[jdim];
          }

          // solve along a 1d slice in solution space
          const size_t m = fSizes[idim];
          for (size_t k = 0; k < m; ++k)
          {
            fBuffer[k] = 0;
            for (size_t l = 0; l < m-2; ++l)
              fBuffer[k] += fInv[idim][k][l]*cs[j + l*strides[idim]];
            fBuffer[k] += fInv[idim][k][m-2]*fBCs[idim].fLeft;
            fBuffer[k] += fInv[idim][k][m-1]*fBCs[idim].fRight;
          }

          for (size_t l = 0; l < m; ++l)
            cs[j + l*strides[idim]] = fBuffer[l];
        }
      }
    }

  protected:
    BoundaryCondition fBCs[ADimension];
    size_t fSizes[ADimension];
    mutable std::vector<double> fBuffer;
    boost::multi_array<double,2> fInv[ADimension];
  };

} // NS Spline

#endif