      ! File: newton_raphson_function.f90
      !
      ! Copyright (c) 2019, LEONARDO COMPANY S.p.A.  All rights reserved.
      !
      ! LEONARDO COMPANY S.p.A. and its licensors retain all intellectual property
      ! and proprietary rights in and to this software, related documentation
      ! and any modifications thereto.  Any use, reproduction, disclosure or 
      ! distribution of this software and related documentation,via any medium,
      ! without an express license agreement from LEONARDO COMPANY S.p.A. is strictly 
      ! prohibited.
      !
      ! Written by: Fabrizio Petri <fabrizio.petri@leonardocompany.com>
      ! Creation date:  
      ! 
      module NewtonRaphsonFunctionClass_mod
      
        !Class Dependency

        implicit none
        
        !Class Visibility
        private
        public :: NewtonRaphsonFunctionClass
        
        !Parameters
        double precision ,private , parameter :: 
     &                    DEFAULT_NUMERIC_PERTURBATION = 1.0d-3
        
        !Class
        type  :: NewtonRaphsonFunctionClass
          
          private
          
          !Public Field
          procedure(INRFunctionDouble),pointer,public,nopass :: m_F_ !Objective Function: Scalar Function on which minimum is searched in analityc form f:= f(x,lambda,t)
          procedure(INRGradientDouble),pointer,public,nopass :: m_G_ !Gradient Function of the scalar function  on which minimum is searched in analytic form {G}:= df/dxi  = {G(x,lambda,t)}
          procedure(INRHessianDouble ),pointer,public,nopass :: m_H_ !Hessian Function of the scalar function  on which minimum is searched in analytic form [H]:= dGj/dxi  = [H(x,lambda,t)]
          
          double precision,  public , pointer :: m_eps(:)

          
         contains
        
          !Public Methods
          procedure, public :: F => F_
          procedure, public :: G => G_
          procedure, public :: H => H_
          
          procedure         :: Dispose => DisposeNewtonRaphsonFunction !Default destructor method for the current class
          procedure, non_overridable   :: DisposeNewtonRaphsonFunction !Destructor method for inherited classes
          
        end type NewtonRaphsonFunctionClass
         


        !Class Method Interfaces  

        
        !Fortran Interface (Abstract Class)
        include 'nr_fdouble_interface.fi'      
        include 'nr_gdouble_interface.fi'   
        include 'nr_hdouble_interface.fi'          
       
        
        !Class Constructor Interface  
        public :: New_NewtonRaphsonFunction
        interface New_NewtonRaphsonFunction
          module procedure newNRF_F
        end interface
      
      contains
      
      ! = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = !
      !Class Destructor

      
      subroutine DisposeNewtonRaphsonFunction(this)
      
        implicit none
        
        class(NewtonRaphsonFunctionClass) :: this
        
         This.m_F_  => null()
         This.m_G_  => null()    
         This.m_H_  => null()  

         
      end subroutine
             
      ! = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = !
      !Class Constructor Private

      
      ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !
      
      function newNRF_F(F_, nx, h)  result(NRF)
        implicit none
        class(NewtonRaphsonFunctionClass) , pointer :: NRF
        procedure(INRFunctionDouble)   :: F_
        integer :: nx !Number of independent variables
        double precision, optional, intent(in),target :: h(:)
        
        
        allocate(NRF)
        
        NRF.m_F_ => F_
        NRF.m_G_ => NRF_GN
        NRF.m_H_ => NRF_HN
        


        
        !Set the perturbation values in a vector
        if (.NOT.Present(h)) then
          allocate(NRF.m_eps(nx))
          NRF.m_eps = DEFAULT_NUMERIC_PERTURBATION
        else
          NRF.m_eps => h
        endif  
        
        
      end function
 

      !Class Private Methods     
      
      function F_(This, x)
        implicit none
        class(NewtonRaphsonFunctionClass)                 :: This
        double precision,           intent(in) :: x(:)
        double precision :: F_

        F_ = This.m_F_(x)
        
      end function
       
      function G_(This, x)
        implicit none
        class(NewtonRaphsonFunctionClass)                 :: This
        double precision,   pointer :: G_(:)
        double precision,           intent(in) :: x(:)


        G_ => This.m_G_(x, This.m_F_, This.m_eps)

      end function

      function H_(This, x)
        implicit none
        class(NewtonRaphsonFunctionClass)                 :: This
        double precision,pointer                          :: H_(:,:)
        double precision,   intent(in)  :: x(:)

        
        H_  => This.m_H_(x, This.m_F_, This.m_G_, This.m_eps)

        
      end function      
       

      !Function for numerical evaluation of gradient
      function NRF_GN(x, F_, h)  result(G)
         implicit none
         
         !class(VectorClass), pointer :: G
         double precision, pointer :: G(:)
         double precision,intent(in) :: x(:)
         double precision, optional, intent(in) :: h(:)        
         procedure(INRFunctionDouble), optional          :: F_
         
         double precision   , allocatable :: xPlusH(:)
         double precision   , allocatable :: xMinusH(:)
         double precision             :: f_hMinus , f_hPlus

         integer                      :: n , i
         
         n = size(x)
         allocate(G(n))
         G = 0.0d0
         
         !Creates x+h and x-h as Vectors
         allocate(xPlusH(n))
         allocate(xMinusH(n))
         

         xPlusH  = x
         xMinusH = x
         
         
         !Gradient evaluation with centered finite differences
         do i = 1 , n
           
           xPlusH(i)  =  xPlusH(i)  + h(i)
           xMinusH(i) =  xMinusH(i) - h(i)
           
           f_hMinus = F_(xMinusH)
           f_hPlus  = F_( xPlusH)
           
           G(i) = (f_hPlus - f_hMinus)/(2.0d0*h(i)) 
         
           xPlusH(i) =  x(i) 
           xMinusH(i) =  x(i) 
           
         end do
         

         
         deallocate(xPlusH)
         deallocate(xMinusH)
         
      end function
      
      ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !
      
      !Function for the numerical evaluation of hessian
      function NRF_Hn(x, F_, G_, h)  result(Hn)
      
         implicit none
         double precision                        , pointer :: Hn(:,:)
         double precision,           intent(in) :: x(:)
	   double precision, optional, intent(in) :: h(:)
         procedure(INRFunctionDouble), optional :: F_
      	 procedure(INRGradientDouble), optional :: G_
         
         double precision, allocatable :: xplush(:)
         double precision, allocatable :: xminush(:)
         double precision, pointer :: gMinus(:)
         double precision, pointer :: gPlus(:)
         integer  n , i
         
         n = size(x)
         allocate(Hn(n,n))
         Hn = 0.0d0
         
         !Creates x+h and x-h as Vectors
         allocate(xPlusH(n))
         allocate(xMinusH(n))

         
         !Copy the x in Vectors
         xPlusH  = x
         xMinusH = x
         
                 
          do i = 1 , n
            
            xPlusH(i)   =  xPlush(i)  + H(i)
            xMinusH(i)  =  xMinush(i) - H(i)
                        
            gMinus =>  G_(xMinusH, F_,h)  !BUGS HERE !!!!
            gPlus  =>  G_( xPlusH, F_,h)
                       

                        
            Hn(:,i) = (gPlus - gMinus)/(2.0d0*h(i))     
            
            xPlusH(i)   =  x(i) 
            xMinusH(i)  =  x(i) 
            
          end do   
          

         
         deallocate(xPlusH)
         deallocate(xMinusH)
                  
      end function
      

      

     
          
      end module NewtonRaphsonFunctionClass_mod
      