! Copyright (c) INRIA
! License: CECILL-B
!
!> @ingroup splines
!> @brief   Derived type for non-uniform B-splines of arbitrary degree
!> @author  Yaman Güçlü  - IPP Garching
!> @author  Edoardo Zoni - IPP Garching

module sll_m_bsplines_non_uniform
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"

  use sll_m_working_precision, only: f64
  use sll_m_bsplines_base    , only: sll_c_bsplines

  implicit none

  public :: &
    sll_t_bsplines_non_uniform

  private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !> Working precision
  integer, parameter :: wp = f64

  !> Abstract type, B-splines
  type, extends(sll_c_bsplines) :: sll_t_bsplines_non_uniform

  contains
    procedure :: init                    => s_bsplines_non_uniform__init
    procedure :: free                    => s_bsplines_non_uniform__free
    procedure :: find_cell               => f_bsplines_non_uniform__find_cell
    procedure :: eval_basis              => s_bsplines_non_uniform__eval_basis
    procedure :: eval_deriv              => s_bsplines_non_uniform__eval_deriv
    procedure :: eval_basis_and_n_derivs => s_bsplines_non_uniform__eval_basis_and_n_derivs

  end type sll_t_bsplines_non_uniform

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !-----------------------------------------------------------------------------
  !> @brief      Initialize non-uniform B-splines object
  !> @param[out] self      non-uniform B-splines
  !> @param[in]  degree    spline degree
  !> @param[in]  periodic  .true. if domain is periodic, .false. otherwise
  !> @param[in]  breaks    array of breakpoints (including xmin and xmax)
  !-----------------------------------------------------------------------------
  subroutine s_bsplines_non_uniform__init( self, degree, periodic, breaks )
    class(sll_t_bsplines_non_uniform), intent(  out) :: self
    integer                          , intent(in   ) :: degree
    logical                          , intent(in   ) :: periodic
    real(wp)                         , intent(in   ) :: breaks(:)

    integer  :: i
    real(wp) :: period ! length of period for periodic B-splines

    ! Public attributes
    self % degree   = degree
    self % periodic = periodic
    self % uniform  = .false.
    self % ncells   = size(breaks) - 1
    self % nbasis   = merge( self%ncells, self%ncells+degree, periodic )
    self % offset   = merge( degree/2   , 0                 , periodic )
    self % xmin     = breaks(1)
    self % xmax     = breaks(self%ncells+1)

    ! Create the knots array from the grid points. Here take the grid points
    ! as knots and simply add to the left and the right the
    ! amount of knots that depends on the degree of the requested
    ! spline. We aim at setting up the indexing in such a way that the
    ! original indexing of 'grid' is preserved, i.e.: grid(i) = knot(i), at
    ! least whenever the scope of the indices defined here is active.
    associate( npoints => self%ncells+1 )

      allocate( self % knots(1-degree:npoints+degree) )

      do i = 1, npoints
        self%knots(i) = breaks(i)
      end do

      ! Fill out the extra nodes
      if ( periodic ) then
        period = breaks(npoints) - breaks(1)
        do i = 1, degree
          self%knots(1-i)       = breaks(npoints-i) - period
          self%knots(npoints+i) = breaks(i+1)       + period
        end do
      else ! open
        do i = 1, degree
          self%knots(1-i)       = breaks(1)
          self%knots(npoints+i) = breaks(npoints)
        end do
      end if

    end associate

  end subroutine s_bsplines_non_uniform__init

  !-----------------------------------------------------------------------------
  !> @brief        Free storage
  !> @param[inout] self non-uniform B-splines
  !-----------------------------------------------------------------------------
  subroutine s_bsplines_non_uniform__free( self )
    class(sll_t_bsplines_non_uniform), intent(inout) :: self

    deallocate( self % knots )

  end subroutine s_bsplines_non_uniform__free

  !----------------------------------------------------------------------------
  !> @brief     Find which grid cell contains the given point
  !> @param[in] self  non-uniform B-splines
  !> @param[in] x     point of interest
  !> results          cell index
  !----------------------------------------------------------------------------
  SLL_PURE function f_bsplines_non_uniform__find_cell( self, x ) result( icell )
    class(sll_t_bsplines_non_uniform), intent(in) :: self
    real(wp)                         , intent(in) :: x
    integer :: icell

    integer :: low, high

    associate( npoints => self%ncells+1 )

      ! check if point is outside of grid, set to outermost cell
      if (x > self%knots(npoints)) then; icell = self%ncells; return; end if
      if (x < self%knots(1)      ) then; icell = 1; return; end if

      ! check if point is exactly on left/right boundary
      if (x == self%knots(1)      ) then; icell = 1          ; return; end if
      if (x == self%knots(npoints)) then; icell = self%ncells; return; end if

      low   = 1
      high  = npoints
      icell = (low + high) / 2
      do while (x <  self%knots(icell) .or. x >= self%knots(icell+1))
         if (x < self%knots(icell)) then
           high = icell
         else
           low  = icell
         end if
         icell = (low + high) / 2
      end do

    end associate

  end function

  !-----------------------------------------------------------------------------
  !> Evaluate value at x of all basis functions with support in local cell
  !> values[j] = B_j(x) for jmin <= j <= jmin+degree
  !>
  !> @param[in]  self    non-uniform B-splines
  !> @param[in]  x       evaluation point
  !> @param[out] values  array of B-splines' values
  !> @param[out] jmin    index of first non-zero B-spline
  !-----------------------------------------------------------------------------
  SLL_PURE subroutine s_bsplines_non_uniform__eval_basis( self, x, values, jmin )
    class(sll_t_bsplines_non_uniform), intent(in   ) :: self
    real(wp)                         , intent(in   ) :: x
    real(wp)                         , intent(  out) :: values(0:)
    integer                          , intent(  out) :: jmin

    integer  :: icell

    integer  :: j, r
    real(wp) :: temp
    real(wp) :: saved

    ! GFortran: to allocate on stack use -fstack-arrays
    real(wp) :: left (1:self%degree)
    real(wp) :: right(1:self%degree)

    ! Check on inputs
    SLL_ASSERT( x > self%xmin - 1.0d-14 )
    SLL_ASSERT( x < self%xmax + 1.0d-14 )
    SLL_ASSERT( size(values) == 1+self%degree )

    ! 1. Compute cell index 'icell'
!    icell = f_bsplines_non_uniform__find_cell( self, x )
    icell = self % find_cell( x )
    SLL_ASSERT( icell >= 1               )
    SLL_ASSERT( icell <= self%ncells     )
    SLL_ASSERT( self%knots(icell) <= x   )
    SLL_ASSERT( x <= self%knots(icell+1) )

    ! 2. Compute index range of B-splines with support over cell 'icell'
    jmin = icell - self%offset

    ! 3. Compute values of aforementioned B-splines
    values(0) = 1.0_wp
    do j = 1, self%degree
       left (j) = x - self%knots(icell+1-j)
       right(j) = self%knots(icell+j) - x
       saved    = 0.0_wp
       do r = 0, j-1
          temp      = values(r) / (right(r+1) + left(j-r))
          values(r) = saved + right(r+1) * temp
          saved     = left(j-r) * temp
       end do
      values(j) = saved
    end do

  end subroutine s_bsplines_non_uniform__eval_basis

  !-----------------------------------------------------------------------------
  !> Evaluate derivative at x of all basis functions with support in local cell
  !> derivs[j] = B_j'(x) for jmin <= j <= jmin+degree
  !>
  !> @param[in]  self    non-uniform B-splines
  !> @param[in]  x       evaluation point
  !> @param[out] derivs  array of B-splines' derivatives
  !> @param[out] jmin    index of first non-zero B-spline
  !-----------------------------------------------------------------------------
  SLL_PURE subroutine s_bsplines_non_uniform__eval_deriv( self, x, derivs, jmin )
    class(sll_t_bsplines_non_uniform), intent(in   ) :: self
    real(wp)                         , intent(in   ) :: x
    real(wp)                         , intent(  out) :: derivs(0:)
    integer                          , intent(  out) :: jmin

    integer  :: icell

    integer  :: j, r
    real(wp) :: temp
    real(wp) :: saved

    ! GFortran: to allocate on stack use -fstack-arrays
    real(wp) :: left (1:self%degree)
    real(wp) :: right(1:self%degree)

    ! Check on inputs
    SLL_ASSERT( x > self%xmin - 1.0d-14 )
    SLL_ASSERT( x < self%xmax + 1.0d-14 )
    SLL_ASSERT( size(derivs) == 1+self%degree )

    ! 1. Compute cell index 'icell' and x_offset
!    icell = f_bsplines_non_uniform__find_cell( self, x )
    icell = self % find_cell( x )
    SLL_ASSERT( icell >= 1               )
    SLL_ASSERT( icell <= self%ncells     )
    SLL_ASSERT( self%knots(icell) <= x   )
    SLL_ASSERT( x <= self%knots(icell+1) )

    ! 2. Compute index range of B-splines with support over cell 'icell'
    jmin = icell - self%offset

    ! 3. Compute derivatives of aforementioned B-splines
    associate( degree => self%degree, degree_real => real(self%degree,wp) )

      ! Compute nonzero basis functions and knot differences
      ! for splines up to degree deg-1 which are needed to compute derivative
      ! First part of Algorithm  A3.2 of NURBS book
      derivs(0) = 1.0_wp
      do j = 1, degree-1
         left (j) = x - self%knots(icell+1-j)
         right(j) = self%knots(icell+j) - x
         saved    = 0.0_wp
         do r = 0, j-1
            ! compute and save bspline values
            temp      = derivs(r)/(right(r+1) + left(j-r))
            derivs(r) = saved + right(r+1) * temp
            saved     = left(j-r) * temp
         end do
         derivs(j) = saved
      end do

      ! Compute derivatives at x using values stored in bsdx and formula
      ! formula for spline derivative based on difference of splines of
      ! degree deg-1
      ! -------
      ! j = 0
      saved = degree_real * derivs(0) / (self%knots(icell+1)-self%knots(icell+1-degree))
      derivs(0) = -saved
      do j = 1, degree-1
         temp      = saved
         saved     = degree_real * derivs(j) / (self%knots(icell+j+1)-self%knots(icell+j+1-degree))
         derivs(j) = temp - saved
      end do
      ! j = degree
      derivs(degree) =  saved

    end associate

  end subroutine s_bsplines_non_uniform__eval_deriv

  !-----------------------------------------------------------------------------
  !> Evaluate value and n derivatives at x of all basis functions with support in local cell
  !> derivs[i,j] = (d/dx)^i B_j(x) for 0 <= i <= n and jmin <= j <= jmin+degree
  !>
  !> @param[in]  self    non-uniform B-splines
  !> @param[in]  x       evaluation point
  !> @param[in]  n       number of required derivatives
  !> @param[out] derivs  array of B-splines' (multiple) derivatives
  !> @param[out] jmin    index of first non-zero B-spline
  !-----------------------------------------------------------------------------
  SLL_PURE subroutine s_bsplines_non_uniform__eval_basis_and_n_derivs( self, x, n, derivs, jmin )
    class(sll_t_bsplines_non_uniform), intent(in   ) :: self
    real(wp)                         , intent(in   ) :: x
    integer                          , intent(in   ) :: n
    real(wp)                         , intent(  out) :: derivs(0:,0:)
    integer                          , intent(  out) :: jmin

    integer  :: icell

    integer  :: j, r
    real(wp) :: temp
    real(wp) :: saved

    integer  :: k, s1, s2, rk, pk, j1, j2
    real(wp) :: d

    ! GFortran: to allocate on stack use -fstack-arrays
    real(wp) :: left (1:self%degree)
    real(wp) :: right(1:self%degree)
    real(wp) :: ndu  (0:self%degree,0:self%degree)
    real(wp) :: a    (0:1          ,0:self%degree)

    ! Check on inputs
    SLL_ASSERT( x > self%xmin - 1.0d-14 )
    SLL_ASSERT( x < self%xmax + 1.0d-14 )
    SLL_ASSERT( n >= 0 )
    SLL_ASSERT( n <= self%degree )
    SLL_ASSERT( size(derivs,1) == 1+n           )
    SLL_ASSERT( size(derivs,2) == 1+self%degree )

    ! 1. Compute cell index 'icell' and x_offset
!    icell = f_bsplines_non_uniform__find_cell( self, x )
    icell = self % find_cell( x )
    SLL_ASSERT( icell >= 1               )
    SLL_ASSERT( icell <= self%ncells     )
    SLL_ASSERT( self%knots(icell) <= x   )
    SLL_ASSERT( x <= self%knots(icell+1) )

    ! 2. Compute index range of B-splines with support over cell 'icell'
    jmin = icell - self%offset

    ! 3. Compute nonzero basis functions and knot differences for splines
    !    up to degree (degree-1) which are needed to compute derivative
    !    Algorithm  A2.3 of NURBS book
    !
    !    21.08.2017: save inverse of knot differences to avoid unnecessary divisions
    !                [Yaman Güçlü, Edoardo Zoni]
    associate( degree => self%degree )

      ndu(0,0) = 1.0_wp
      do j = 1, degree
         left(j)  = x - self%knots(icell+1-j)
         right(j) = self%knots(icell+j) - x
         saved    = 0.0_wp
         do r = 0, j-1
            ! compute inverse of knot differences and save them into lower triangular part of ndu
            ndu(j,r) = 1.0_wp / (right(r+1) + left(j-r))
            ! compute basis functions and save them into upper triangular part of ndu
            temp     = ndu(r,j-1) * ndu(j,r)
            ndu(r,j) = saved + right(r+1) * temp
            saved    = left(j-r) * temp
         end do
         ndu(j,j) = saved
      end do
      derivs(0,:) = ndu(:,degree)

      do r = 0, degree
         s1 = 0
         s2 = 1
         a(0,0) = 1.0_wp
         do k = 1, n
            d  = 0.0_wp
            rk = r-k
            pk = degree-k
            if (r >= k) then
               a(s2,0) = a(s1,0) * ndu(pk+1,rk)
               d = a(s2,0) * ndu(rk,pk)
            end if
            if (rk > -1) then
               j1 = 1
            else
               j1 = -rk
            end if
            if (r-1 <= pk) then
               j2 = k-1
            else
               j2 = degree-r
            end if
            do j = j1, j2
               a(s2,j) = (a(s1,j) - a(s1,j-1)) * ndu(pk+1,rk+j)
               d = d + a(s2,j) * ndu(rk+j,pk)
            end do
            if (r <= pk) then
               a(s2,k) = - a(s1,k-1) * ndu(pk+1,r)
               d = d + a(s2,k) * ndu(r,pk)
            end if
            derivs(k,r) = d
            j  = s1
            s1 = s2
            s2 = j
         end do
      end do
      r = degree
      do k = 1, n
         derivs(k,:) = derivs(k,:) * r
         r = r * (degree-k)
      end do

    end associate

  end subroutine s_bsplines_non_uniform__eval_basis_and_n_derivs

end module sll_m_bsplines_non_uniform