*** Define subroutine to generate output layer error for one pattern,
*** as well as the derivate of the error in respect to the neuronal output
***
***
***
****Conventions:
***
*** Most of the data structures are larger than declared;
*** however, only their relevant subsection is passed to the
*** subroutine and hence considered.
***
*** nnerror:   subroutine evaluating the output layer's error.
***
*** nnweight:  system specific weight adaptions for fit
***
*** pat_err:   segment of the total error vector for one pattern.
***
*** len_in:    Actual number of input neurons
*** len_out:   Actual number of output neurons
***
*** pat_in:    Input pattern
***  pat_in(i):   value of ith input neuron for given pattern
***
*** pat_out:   Desired output pattern
***  pat_out(i):  desired value of ith output neuron for given pattern
***
*** L:         Final layer, starting at the first neuron
***************************************************************************
****Note:
*** neu_in and neu_out, if needed, should ideally not be hardcoded.
*** Instead they should be passed by the main program to ensure
*** matching dimensionalities.
***
***************************************************************************

      subroutine nnerror(pat_err,pat_in,pat_out,nn_out)
      implicit none
!     Generate error vector from adiabatic model.
!
!     nn_out:        output from neural network
!     pat_err:       error vector for single pattern
!

      include 'nnparams.incl'
      include 'nndbg.incl'
      include 'nncommon.incl'

      double precision nn_out(maxnout)
      double precision pat_in(maxnin),pat_out(maxpout),pat_err(maxpout)

      double precision adiaoutp(maxpout)

      integer j

      call nnadia(pat_in,nn_out,adiaoutp)

      do j=1,inp_out
         pat_err(j) = pat_out(j) - adiaoutp(j)
      enddo

      end

!--------------------------------------------------------------------------------------

      subroutine nnweight(wterr,pat_out)
      implicit none
!     Evaluate system specific weighting for 1 pattern.

      include 'params.incl'
      include 'nnparams.incl'
      include 'nncommon.incl'

      !include 'JTmod.incl'

      double precision wterr(maxpout),pat_out(maxpout)

      !double precision eref(nstat)
      !double precision wten(3)
       integer i 
      !integer j

      !wten=1

      !eref(1)=E0_sig
      !eref(2:3)=E0_pi
!     eref(1:3)=eref(1:3)+600.d0*icm2hart
      !eref(1:3)=eref(1:3)+1000.d0*icm2hart

!!     kill state 2 and 3
!      wten(2:3)=0

      !do j=1,3
!        weighting of energies
      !   wterr(j)=wdamp(pat_out(j)-eref(j))*wterr(j)*wten(j)
      !enddo
      do i=1,size(wterr) 
            if (wterr(i) .eq. 0.0d0) then 
                  wterr(i) =0.0d0 
            else
                  wterr(i) = 1.0d0
            endif 
      enddo  
      pat_out=pat_out
      !wterr = 1.0d0 

      contains
      double precision function wdamp(dE)
      implicit none
      include 'nnparams.incl'
      double precision dE

!     Weight decay rate
      double precision, parameter :: unit=eV2hart
!     ln(3)/2 = artanh(1/2)
!     double precision, parameter :: alpha=0.5d0*log(3.d0)/unit
!     (1+tanh(-x))/2 ~ exp(-2x)
!     double precision, parameter :: alpha=0.5d0*log(2.d0)/unit
!
      double precision, parameter :: alpha=log(10.d0)/unit

!     cutoff
      double precision, parameter :: minweight=1.d-4

!     asymptotically,
!     wdamp=(1.d0+tanh(-alpha*dE))*0.5d0 + minweight
      wdamp=1.0d0

      if (dE.lt.0) then
         wdamp=1.d0
      else if (dE.lt.(4*eV2hart)) then
         wdamp=exp(-alpha*dE)+minweight
      else
         wdamp=0
      endif

      end function wdamp
      end subroutine

!--------------------------------------------------------------------------------------

      subroutine nnoutgrad(ad_grads,pat_in,nn_out)
      implicit none
!     Compute the derivative of each output value (i.e. adiab. energies)
!     with respect to the Neural Network output-neurons.
!
!     nn_out:           output from neural network
!     ad_grads:         derivatives as described above
!       ad_grads(:,i):  gradient of energy i
!     eps:              finite differences for derivatives
!       eps(i):           finite difference for output neuron i
!     dis_out:          ANN output displaced by finite differences
!       dis_out(:,1):     equivalent to L + eps
!       dis_out(:,2):     equivalent to L - eps
!     dis_ad:           yielded adiabatic output for current displacements
!
!     ddelta:        factor used to determine eps.

      include 'nnparams.incl'
      include 'nncommon.incl'
      include 'nndbg.incl'

      double precision nn_out(maxnout)
      double precision pat_in(maxnin)
      double precision ad_grads(maxnout,maxpout)

      double precision eps(maxnout),dis_out(maxnout,2)
      double precision dis_ad(maxpout,2)
      double precision ddelta

      integer n,j,k

      parameter (ddelta = 1.0D-2) !reduce again if appropriate

!     determine appropriate finite differences for each parameter:
      do k=1,len_out
         eps(k)=abs(nn_out(k))*ddelta
         if (eps(k).lt.zero) then
            eps(k)=zero
         endif
      enddo

      do n=1,len_out
         do k=1,2
!           copy ANN-output
            do j=1,len_out
               dis_out(j,k) = nn_out(j)
            enddo
         enddo

!        apply finite difference for output-neuron n
         dis_out(n,1) = dis_out(n,1) + eps(n)
         dis_out(n,2) = dis_out(n,2) - eps(n)

         do k=1,2
!           get energies and ci-values
            call nnadia(pat_in,dis_out(1,k),dis_ad(1,k))
         enddo

!        apply finite differences to generate numerical gradient
         do k=1,inp_out
            ad_grads(n,k) = (dis_ad(k,1)-dis_ad(k,2))/(2.0d0*eps(n))
         enddo
      enddo

      end