program steepest_ls

! line_search = Option switch for line seaching (1=2-pt,0=simple) 
!      = 0 = none 
!      = 1 = simple 
!      = 2 = 2-point parabolic 
!      = 3 = 2-point parabolic 
 
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 steepest_mod
use base_spline_mod
use fx2d
use lloyd_vel

implicit none

character*80                                  :: hor1
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: vmin,vmax,p_up,p_lo,tol(2)
integer                                       :: ni,i,ishot,j,k,iters,titer,stat1
real, allocatable                             :: velint(:),taper(:),pick1(:)
real, allocatable                             :: R_image(:,:,:)
real, allocatable                             :: vel0(:)
logical                                       :: cond,bounds,cbounds
double precision, allocatable                 :: x(:), g(:)
real, allocatable                             :: l(:), u(:)


integer                                       :: max_iters, line_search, a
real                                          :: err, err1, epsteep, factor, lambda_min, tola, nrm, hh
real                                          :: lambda,lam1,lam2,lam3,lamc,lam4,gamma
double precision                              :: ft,fn,f2,gt,gn,g2,h20,h320,h32
double precision, allocatable                 :: xn(:), dx(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()

 call from_param("line_search",line_search,3)
 call from_param("max_iters",max_iters,10)
 call from_param("gamma",gamma,1.)

allocate(R_image(img%xm%n,nsecdim,img%z%n));R_image=0.
allocate(taper(kxm%n));taper=0.
allocate(g(img%z%n*img%xm%n));g=dble(0.)
allocate(velint(v%z%n*v%xm%n),vel0(img%z%n*img%xm%n))
allocate(x(img%z%n*img%xm%n));x=dble(0.)
allocate(xn(img%z%n*img%xm%n));xn=dble(0.)
allocate(dx(img%z%n*img%xm%n));dx=dble(0.)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 call owtaper(tapertr,taper)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Read velocities
!
 call sreed("vel",velint,v%z%n*v%xm%n*4)
 call owgrab_vel(velint,vel0,.true.,ni)
 deallocate(velint)

!----------------------------------------------
! Initial velocity model = background velocity
 x=dble(vel0)

!----------------------------------------------
! Set bounds
 call from_param("bounds",bounds,.false.)
 if (bounds) then
    allocate(l(img%z%n*img%xm%n));l=0.
    allocate(u(img%z%n*img%xm%n));u=0.
!    call from_param("cbounds",cbounds,.false.)
!    if (cbounds) then
!       call from_param("vmin",vmin)
!       call from_param("vmax",vmax)
!       l=1/vmax
!       u=1/vmin
!    else
       call from_param("p_up",p_up)
       call from_param("p_lo",p_lo)
       if (p_up >= 1.) p_up=1+p_up/100
       if (p_lo >= 1.) p_lo=1-p_lo/100
       l=p_lo*x
       u=p_up*x
!    end if
 end if 
!----------------------------------------------
! Horizon with no velocity update
if (exist_file("horv0")) then
   call from_param('horv0',hor1)
   open(51,file=hor1,status='old')
   allocate(pick1(img%xm%n));pick1=9999.
   do i=1,img%xm%n
      read(51,*,end=1) pick1(i)
   end do
1 close(51)
end if

! ********************************************** 
!  Initialization of code variables 
 iters=1
 err=1
 epsteep=1e-8 
 tol(1)=1e-8
 tol(1)=1e-5

! Line search variables 
 factor=0.5; 
 lambda_min=0.00001; 
! ********************************************** 

 xn=dble(vel0)
 x=xn
 call to_history("n4",iters,"Velout")
 call srite("Velout",sngl(x),4*img%z%n*img%xm%n)
 call fgeval_init(sngl(x),taper,iters,pick1)
 call feval(fn,R_image)
 call fgeval_close()
!   Set tolerance for stopping criterium 
 tola=tol(1)*norm(g)+tol(2)
 lam3=1
 err1=tiny(1.)
! ********************************************** 
 do while (err > err1 .and. iters <= max_iters) ! Main code loop starts here  
 
    nrm=norm(xn) 
    call fgeval_init(sngl(x),taper,iters,pick1)
    call geval(g,R_image)
    call fgeval_close()
!  Error 
    err=norm(g); 
    if (iters==1) err1=err*0.05
 
!  Update the Newton search direction 
    dx=-g/err*nrm*gamma
    call to_history("n4",iters,"grad")
    call srite("grad",sngl(dx),4*img%z%n*img%xm%n)


! ******************************************************************  
! Line Searching below 
! ******************************************************************  
    if (line_search==0) then ! No line search 
       lambda=1; 
    else if (line_search==1.or.line_search==2) then ! Simple OR 2-point parabolic line search 
       lam1=1; 
       gn=fn*fn; 
       x=xn+dx

!       call ow_smoovel(x,vel0,pick1)
       call fgeval_init(sngl(x),taper,iters,pick1)
       call feval(ft)
       call fgeval_close()
       gt=ft*ft; 
   
       do while (gt > gn.and.lam1 > lambda_min) 
          lam1=factor*lam1; 
          x=xn+lam1*dx
!          call ow_smoovel(x,vel0,pick1)
          call fgeval_init(sngl(x),taper,iters,pick1)
          call feval(ft)
          call fgeval_close()
          gt=ft*ft
       end do     
 
       if (line_search==2) then ! 2-point parabolic line search 
          lamc=(gn*lam1**2)/(gt*(2*lam1-1));  
          x=xn+lamc*dx
!          call ow_smoovel(x,vel0,pick1)
          call fgeval_init(sngl(x),taper,iters,pick1)
          call feval(f2)
          call fgeval_close()
          g2=f2*f2
    
          if (g2 < gt) then
             lambda=lamc
             fn=f2
             gn=g2
          else 
             lambda=lam1
             fn=ft
             gn=gt
          end if
          write(0,*) "02) fn = ",fn," ft = ",ft," f2 = ",f2
       else  ! Line_search==1 => simple line search 
 
          lambda=lam1; 
          fn=ft; 
          gn=gt; 
     
       end if 
 
    else if (line_search==3) then ! 3-point parabolic line search 
 
       lam3=1
999    x=xn+lam3*dx
       if (any(mask=(x <= l.or.x >= u))) then
          lam3=factor*lam3
          write(0,*) "1) Bound reached - recomputing lam3", lam3
          goto 999
       end if
       gn=fn*fn
       x=xn+dx
!       call ow_smoovel(x,vel0,pick1)
       call fgeval_init(sngl(x),taper,iters,pick1)
       call feval(ft)
       call fgeval_close()
       gt=ft*ft
       write(0,*) " gt = ", gt," gn = ", gn
       do while (gt > gn.and.lam3 > lambda_min) 
1000      lam3=factor*lam3
          x=xn+lam3*dx
          if (any(mask=(x <= l.or.x >= u))) then
             write(0,*) "2) Bound reached - recomputing lam3", lam3
             goto 1000
          end if
!          call ow_smoovel(x,vel0,pick1)
          call fgeval_init(sngl(x),taper,iters,pick1)
          call feval(ft)
          call fgeval_close()
          gt=ft*ft
          write(0,*) "01) fn = ",fn," ft = ",ft
       end do     
     
       lam2=factor*lam3; 
       x=xn+lam2*dx
!       call ow_smoovel(x,vel0,pick1)
       call fgeval_init(sngl(x),taper,iters,pick1)
       call feval(f2)
       call fgeval_close()
       g2=f2*f2
       h20=(g2-gn)/lam2; 
       h32=(gt-g2)/(lam3-lam2); 
       h320=(h32-h20)/lam3; 

       lamc=0.5*(lam2-h20/h320); 
       x=xn+lamc*dx
!       call ow_smoovel(x,vel0,pick1)
       call fgeval_init(sngl(x),taper,iters,pick1)
       call feval(f2)
       call fgeval_close()
       g2=f2*f2
   
       write(0,*) "02) fn = ",fn," ft = ",ft," f2 = ",f2
       if (g2 < gt) then
          lambda=lamc; 
          fn=f2; 
          gn=g2; 
       else 
          lambda=lam3; 
          fn=ft; 
          gn=gt; 
       end if   
 
    end if 
! ****************************************************************** 
! Line search complete 
! ****************************************************************** 
    xn=xn+lambda*dx;       
! Update Newton iterate 
    x=xn
!    call ow_smoovel(x,vel0,pick1)
    write(0,*) 'Iteration = ',iters,' error=',err
    write(0,*) "03) fn = ",fn

    iters=iters+1;  !Update iterations counter 
    R_image=0.
    call fgeval_init(sngl(x),taper,iters,pick1)
    call feval(fn,R_image)
    call fgeval_close()

    call to_history("n4",iters,"Velout")
    call srite("Velout",sngl(x),4*img%z%n*img%xm%n)
 
 end do
 
! Success statement based on error and iters ******** 
 if (err > err1.or.iters == max_iters) then
    write(0,*) 'Convergence may not have occurred; because either the error is greater than 20% grad_0 or max_iters was reached' 
 else if (err < err1) then
    write(0,*)  'The error is less than 20% grad_0; convergence has likely been achieved' 
 end if 

 if (allocated(R_image)) deallocate(R_image)
 if (allocated(taper)) deallocate(taper)
 if (allocated(g)) deallocate(g)
 if (allocated(x)) deallocate(x)
 if (allocated(xn)) deallocate(xn)
 if (allocated(dx)) deallocate(dx)

stop
end program !End of main program
