#ifndef _UniformBSpline_Spline_Uniform_h_
#define _UniformBSpline_Spline_Uniform_h_

#include <cmath>
#include <vector>
#include <algorithm>
#include <limits>
#include <iostream>
#include <utl/AugerException.h>

namespace Spline {

  typedef unsigned char dim_t;

  template <dim_t ADimension, typename AKnotVector, typename ABasisFunction>
  class Function;

  namespace Uniform {

    class KnotVector {
    public:
      KnotVector()
        : fN(0), fStart(0), fDelta(1)
      {}

      KnotVector(const size_t n, const double start, const double stop)
        : fN(n), fStart(start), fDelta(stop-start)
      {
        if (n<2)
          throw utl::InvalidConfigurationException("KnotVector: size < 2");

        if (stop < start)
          throw utl::InvalidConfigurationException("KnotVector: stop < start");
      }

      inline
      void
      PushBack(const double x)
      {
        if (fN == 0)
          fStart = x;
        fDelta = x - fStart;
        ++fN;
      }

      inline
      size_t
      Locate(const double x)
        const
      {
        if (x < fStart)
          return 0;

        if (x >= fStart + fDelta)
          return fN-2;

        return size_t((x - fStart)/fDelta*(fN-1));
      }

      inline
      double
      operator[](const size_t i)
        const
      {
        return fStart + double(i)/(fN-1)*fDelta;
      }

      inline
      size_t
      GetSize()
        const
      {
        return fN;
      }

      inline
      double
      GetStart()
        const
      {
        return fStart;
      }

      inline
      double
      GetStop()
        const
      {
        return fStart + fDelta;
      }

    protected:
      size_t fN;
      double fStart;
      double fDelta;
    };

    class BasisFunction {
    public:
      BasisFunction()
        : fKnotPtr(NULL)
      {}

      BasisFunction(const KnotVector* xknot, int j)
        : fKnotPtr(xknot),
          fIndex(j-3) // shift to internal index
      {}

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

      BasisFunction&
      operator=(const BasisFunction& other)
      {
        fKnotPtr = other.fKnotPtr;
        fIndex = other.fIndex;
        return *this;
      }

      inline
      double
      operator()(const double x, const char derivative = 0)
        const
      {
        if (x < GetStart())
          return 0.0;

        if (x > GetStop())
          return derivative == -1 ? Internal(fKnotPtr->Locate(x), GetStop(), -1) : 0.0;

        return Internal(fKnotPtr->Locate(x), x, derivative);
      }

      inline
      double
      GetStart()
        const
      {
        return (*fKnotPtr)[std::max(fIndex,0)];
      }

      inline
      double
      GetStop()
        const
      {
        return (*fKnotPtr)[std::min(fIndex+4,int(fKnotPtr->GetSize())-1)];
      }

    protected:

      double
      Internal(const size_t i, const double x, const char derivative)
        const
      {
        const KnotVector& xknot = *fKnotPtr;

        const double dx = xknot[1] - xknot[0];
        const double z = (x - xknot[i])/dx;

        double t0,t1,t2,t3;
        switch (i-fIndex)
        {
          case 3:
            t0 =  1.0/6.0; t1 = -0.5    ; t2 =  0.5    ; t3 = -1.0/6.0;
            break;
          case 2:
            t0 =  2.0/3.0; t1 =  0.0    ; t2 = -1.0    ; t3 =  0.5    ;
            break;
          case 1:
            t0 =  1.0/6.0; t1 =  0.5    ; t2 =  0.5    ; t3 = -0.5    ;
            break;
          case 0:
            t0 =  0.0    ; t1 =  0.0    ; t2 =  0.0    ; t3 =  1.0/6.0;
            break;
          default:
            t0 =  0.0    ; t1 =  0.0    ; t2 =  0.0    ; t3 =  0.0    ;
        }

        switch(derivative)
        {
          case -1:
          {
            double result = z*(t0 + z*(t1/2 + z*(t2/3 + z*t3/4)))*dx;
            for (size_t j=std::max(fIndex,0); j<i; ++j)
            {
              switch (j-fIndex)
              {
                case 3:
                  t0 =  1.0/6.0; t1 = -0.5    ; t2 =  0.5    ; t3 = -1.0/6.0;
                  break;
                case 2:
                  t0 =  2.0/3.0; t1 =  0.0    ; t2 = -1.0    ; t3 =  0.5    ;
                  break;
                case 1:
                  t0 =  1.0/6.0; t1 =  0.5    ; t2 =  0.5    ; t3 = -0.5    ;
                  break;
                case 0:
                  t0 =  0.0    ; t1 =  0.0    ; t2 =  0.0    ; t3 =  1.0/6.0;
                  break;
                default:
                  t0 =  0.0    ; t1 =  0.0    ; t2 =  0.0    ; t3 =  0.0    ;
              }
              result += (t0 + (t1/2 + (t2/3 + t3/4)))*dx;
            }
            return result;
          }
          case 0:
            return t0 + z*(t1 + z*(t2 + z*t3));
          case 1:
            return (t1 + z*(2*t2 + z*3*t3))/dx;
          case 2:
            return (2*t2 + z*6*t3)/(dx*dx);
          case 3:
            return 6*t3/(dx*dx*dx);
          default:
            return 0;
        }
      }

      const KnotVector* fKnotPtr;
      int fIndex;

      template <dim_t ADimension, typename AKnotVector, typename ABasisFunction>
      friend class ::Spline::Function;
    };

  } // NS Uniform

} // NS Spline

#endif