module fgevaluation
use sep
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
use base_spline_mod
implicit none
real, allocatable                        :: taper(:)
contains
!===========================================
subroutine owfgeval_init()
end subroutine
!===========================================
subroutine taper_set()
 allocate(taper(kxm%n))
 call owtaper(tapertr,taper)
end subroutine
!===========================================
subroutine taper_close()
 deallocate(taper)
end subroutine
!===========================================
subroutine wolfe_search2(x,p,alpha0,c1,c2,f0,df0,alpha,fval,dfval,nfe,nge)

real,intent(in)            :: f,df(:),x(:),p(:),alpha0,c1,c2,f0,df0(:)
real, intent(out)          :: alpha0,fval,dfval
integer, intent(out)       :: nfe, nge

! internal variables
real                       :: slope0,slope,xplus(size(x)),old_alpha,old_fval,old_dfval,old_slope
real                       :: temp,alpha_lo,alpha_hi,f_lo,f_hi,df_lo,slope_lo,a,b,dalpha,c,alpha
integer                    :: lsiter,iter_cnt,firsttrip

! [ALPHA,FVAL,DFVAL,NFE,NGE] = wolfe_search2(F,DF,PARAMS,X,P,ALPHA0,C1,C2,F0,DF0,TRACE)
! Performs search for point satisfying the two strong Wolfe conditions:
!
!	f(x+alpha.p) <= f(x) + c1.alpha p'.grad f(0)		(WC1)
!	|p'.grad f(x+alpha.p)| <= c2|p'.grad f'(x)|		(WC2b)
!
! ALPHA0 is the initial value of alpha in the search.  If ALPHA0
!	satisfies the Wolfe conditions, then the search stops there. 
! F0 is f(x) at the initial x.
! DF0 is grad f(x) at the initial x.
! DF(PARAMS,X) is the gradient of F(PARAMS,X) with
!	respect to X
! Note that PARAMS is the parameter for F: F(PARAMS,X)
! If TRACE is non-zero then print out information about the process
!
! The algorithm follows that of Wright & Nocedal "Numerical Optimization",
! pp. 60-62, section 3.4
!
! Returns
!	ALPHA	-- step length parameter
!	FVAL	-- function value at end
!	DFVAL	-- gradient value at end
!	NFE	-- # function evaluations
!	NGE	-- # gradient evaluations

!
! Check parameters
!
alpha = 0
fval = f0
dfval = df0
nfe = 0
nge = 0

if (not(0<c1.and.c1<c2.and.c2<1)) then
  write(0,*) 'wolfe_search2: Error: Need 0 < c1 < c2 < 1'
  write(0,*) 'wolfe_search2: c1 = ',c1,' c2 = ', c2
  return
end if

! Initialization

slope0 = dot_product(p,df0)
if (slope0 >= 0) then
  write(0,*) 'wolfe_search2: Error: Need a descent direction'
  write(0,*) 'wolfe_search2: p*grad f(x) = ',slope0
  return
end

!
! Bracketing phase
!
!write(0,*) 'wolfe_search2: Bracketing phase: alpha = ', alpha0
alpha = alpha0
old_alpha = 0
old_fval = f0
old_dfval = df0
old_slope = slope0
write(0,*) 'wolfe_search2: f(x) = ',f0,' p*grad f(x) = ',slope0

! Main loop
firsttrip = 1
lsiter=1
do while (lsiter < ls_iter ) !!! forever do...

  xplus = x+alpha*p

  call taper_set()
  call owfeval(xplus,fval,B_image)
  nfe = nfe + 1
  write(0,*) 'wolfe_search2: alpha = ',alpha,' f(x+alpha*p) = ',fval
  if (fval > f0 + c1*alpha*slope0 .or. fval >= old_fval ) then
     write(0,*) 'wolfe_search2: (WC1) failed or f increased'
     exit
  end if
  write(0,*) 'wolfe_search2: (WC1) holds & f decreased'

  call owgeval(xplus,dfval,B_image)
  call taper_close()
  nge = nge + 1
  slope = p*dfval
  write(0,*) 'wolfe_search2: p''.grad f(x+alpha*p) = ', slope
  if ( abs(slope) <= c2*abs(slope0) ) then
     write(0,*)  'wolfe_search2: (WC2) holds'
    return
  end if
  if ( slope >= 0 ) then
    write(0,*) 'wolfe_search2: f''(alpha) >= 0'
    exit
  end

  ! Update variables -- note no upper limit on alpha
  temp = alpha
  ! alpha = alpha + old_alpha
  alpha = 2*alpha
  old_alpha = temp
  old_fval = fval
  old_slope = slope
  lsiter=lsiter+1
end do

write(0,*) 'wolfe_search2: Entering ''zoom'' phase.'

!
! "Zoom" phase
!
alpha_lo = old_alpha
alpha_hi =     alpha
f_lo     = old_fval
f_hi     =     fval
df_lo    = old_dfval
slope_lo = old_slope
!!df_hi   =     dfval
!!slope_hi=     slope

write(0,*) 'wolfe_search2: zoom phase: alpha_lo = ',alpha_lo,' alpha_hi = ',alpha_hi
iter_cnt = 0
lsiter=1
do while (lsiter < ls_iter ) !!! forever do...

  ! form quadratic interpolant of function values on alpha_lo, alpha_hi
  ! and the derivative at alpha_lo...
  ! and find min of interpolant within interval [alpha_lo,alpha_hi]
  a = f_lo
  b = slope_lo
  dalpha = alpha_hi-alpha_lo
  c = (f_hi - f_lo - dalpha*slope_lo)/dalpha**2
  ! iter_cnt
  if ( c <= 0 .or. mod(iter_cnt,3) == 2 ) then
    ! Use bisection
    alpha = alpha_lo + 0.5*dalpha
    write(0,*) 'wolfe_search2: using bisection: c = ', c
  else
    ! Use min of quadratic
    alpha = alpha_lo - 0.5*b/c
    write(0,*) 'wolfe_search2: using quadratic: a = ', a,' b = ',b, ' c = ', c
  end if
  write(0,*) 'wolfe_search2: alpha = ', alpha

  ! main part of loop
  xplus = x + alpha*p

  call taper_set()
  call owfeval(xplus,fval,B_image)
  call taper_close()
  nfe = nfe + 1
  write(0,*) 'wolfe_search2: f(x+alpha*p) = ', fval
  if (fval > f0 + c1*alpha*slope0 .or. fval >= f_lo ) then
     write(0,*) ('wolfe_search2: zoom: (WC1) fails or f increased')
     write(0,*) ('wolfe_search2: zoom: update alpha_hi')
     alpha_hi = alpha
     f_hi   = fval
  else
    write(0,*) 'wolfe_search2: zoom: (WC1) holds and f decreased'

    call taper_set()
    call owgeval(xplus,dfval,B_image)
    call taper_close()
    nge = nge + 1
    slope = p*dfval
    write(0,*) 'wolfe_search2: p''.grad f(x+alpha*p) = ', slope
    if (abs(slope) <= c2*abs(slope0)) then
       write(0,*) 'wolfe_search2: zoom: (WC2b) holds'
       return
    end if
    if ( slope*dalpha >= 0 ) then
       write(0,*) 'wolfe_search2: zoom: alpha_hi <- alpha_lo & update alpha_lo'
       alpha_hi = alpha_lo
       f_hi     = f_lo
       !!! df_hi  = df_lo
       !!! slope_hi= slope_lo
       alpha_lo = alpha
       f_lo     = fval
       df_lo    = dfval
       slope_lo = slope
    else
       write(0,*) 'wolfe_search2: zoom: update alpha_lo'
       alpha_lo = alpha
       f_lo     = fval
       df_lo    = dfval
       slope_lo = slope
    end if
  end if
  ! Update iteration count
  iter_cnt = iter_cnt + 1
  lsiter=lsiter+1
end do
end subroutine

!===========================================
subroutine owfeval(x,fun,B_image)

real                                          :: start, finish, start1, finish1
real                                          :: B_image(:,:,:)
integer                                       :: ni,i,ishot,j,k,iter,titer
real, allocatable                             :: deltaI(:,:,:),illum(:,:),illumination(:,:)
real, allocatable                             :: vel(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:)

allocate(deltaI(kxm%n,nsecdim,img%z%n));deltaI=0.
allocate(illum(kxm%n,img%z%n),vel(img%z%n*kxm%n));illum=0.;vel=0
allocate(illumination(img%xm%n,img%z%n));illumination=0.
start1=0.;finish1=0.
!-------------------------------------------
! Computing background image
!
       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0
          if (allocated(source)) deallocate(source); allocate(source(rec%h%n,sou%w%n))
          source=0
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n))
             source=0
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=x
             ni=1
          else
             call owgrab_vel(x,vel,.false.,ni,ishot)
          end if
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
          recv=0;deltaI=0;illum=0.

          call sreed("rec",recv,8*rec%h%n*rec%w%n)
          call owwei_init()
          call owwei(taper,recv,source,vel,deltaI,illum)

          if (areal) then
             B_image=B_image+deltaI
             illumination=illumination+illum
          else
             B_image(ni:ni+kxm%n-1,:,:)=B_image(ni:ni+kxm%n-1,:,:)+deltaI
             illumination(ni:ni+kxm%n-1,:)=illumination(ni:ni+kxm%n-1,:)+illum
          end if
          call owexec_time(finish)

          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       call owexec_time(finish1)
       write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",titer,"bkgimg")
       call srite("bkgimg",B_image,4*img%z%n*nsecdim*img%xm%n)

!-------------------------------------------
!
! Computing DSO
!
       call owdso(.false. , B_image)
       call owdso(.true.  , B_image)
       call to_history("n4",titer,"Dimage")
       call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)

!-------------------------------------------
!  Computing the objective function
!
       fun=0.
       do i=1,img%z%n
          do j=1,img%sh%n
             do k=1,img%xm%n
                fun = fun + B_image(k,j,i)*B_image(k,j,i)
             end do
          end do
       end do
       
       write(2,113) " iter = ",iter," titer= ",titer," J(x) = ", fun

if (allocated(deltaI)) deallocate(deltaI)
if (allocated(illum)) deallocate(illum)
if (allocated(illumination)) deallocate(illumination)
if (allocated(vel)) deallocate(vel)
if (allocated(recv)) deallocate(recv)
if (allocated(source)) deallocate(source)
if (allocated(source1)) deallocate(source1)

end subroutine
!===========================================
subroutine owgeval(x,g,B_image)

logical                                       :: cond,bounds,cbounds
double precision                              :: g(:)
real                                          :: B_image(:,:,:)
real                                          :: start, finish, start1, finish1, start2, finish2
integer                                       :: ni,i,ishot,iz,j,k,iter,titer,i1,i2,i3
real, allocatable                             :: scatS(:),deltaS(:),vel(:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:)
double precision,allocatable                  :: gspl(:)


       if (axkn.ne.99999 .or. azkn.ne.99999) then
          stat1 = base_spln2d_init(img%xm,img%z,axk,azk)
          allocate(gspl(size(g)));gspl=dble(0.)
!          stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
       end if

       start1=0.;finish1=0.;g=0.

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0
          if (allocated(source)) deallocate(source); allocate(source(rec%h%n,sou%w%n))
          source=0
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if

       scatS=0.
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n))
             source=0
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=x
             ni=1
          else
             call owgrab_vel(sngl(x),vel,.false.,ni,ishot)
          end if
             
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
          recv=0;deltaS=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)

          call owwemva_init(ishot)
          call owwemva(.true.,taper,recv,source,vel,deltaS,B_image(ni:ni+kxm%n-1,:,:))
          if (areal) then
             scatS=scatS+deltaS
          else
             do j=ni,ni+kxm%n-1;do k=1,img%z%n
                scatS((k-1)*img%xm%n+j)=scatS((k-1)*img%xm%n+j)+deltaS((k-1)*kxm%n+j-ni+1)
             end do;end do
          end if

          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       g=dble(scatS)

       write(0,*) "main => ",maxval(g),minval(g)

       if (axkn.ne.99999 .or. azkn.ne.99999) then
          call owexec_time(start2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, g)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, g)
          call owexec_time(finish2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
       end if

       write(0,*) "main => ",maxval(g),minval(g)
       write(3,112) "SMOOTHED GRADIENT IN ",(finish2-start2)/60

       call owexec_time(finish1)
       write(3,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",titer,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n)

       if (allocated(deltaS)) deallocate(deltaS)
       if (allocated(scatS)) deallocate(scatS)
       if (allocated(vel)) deallocate(vel)
       if (allocated(recv)) deallocate(recv)
       if (allocated(source)) deallocate(source)
       if (allocated(source1)) deallocate(source1)
       if (allocated(gspl)) deallocate(gspl)

end subroutine
!===========================================
end module
