program InterpRMS
!--------------------------------------------------------------------------------
! Program InterpRMS
!--------------------------------------------------------------------------------
! Usage:  InterpRMS.x o1= d1= n1= weight= < pickfile.H > interp.H
!--------------------------------------------------------------------------------
! Description:  Takes a file of time coordinates/velocity pick pairs, and uses
!               trivial regularized least squares optimization to estimate RMS
!               velocity on a regular grid defined by o1, d1, n1.   
!               Regularization is 1-D gradient.
!--------------------------------------------------------------------------------
! Parameters:  n1 - Size of output model grid.
!              o1 - Origin of output model grid.
!              d1 - Grid spacing of output model grid.
!              nframe - Number of equiprobable realizations.
!              randscale - approx variance of random vector
!              eps - epsilon for missing data problem.
!              niter - number of iterations.
!--------------------------------------------------------------------------------
! Input: (pickfile.H) A file of time coordinates/velocity pick pairs. n2=2
!--------------------------------------------------------------------------------
! Output: (interp.H) RMS velocity on a regular grid defined by o1, d1, n1.
!           (weight) By default, puts picks on regular grid (nearest neighbor) and
!                     outputs.
!--------------------------------------------------------------------------------
  use reg_solver_new
  use cgstep_mod
  use random_mod
  use sep
  use bin1
  use igrad1

  implicit none
  integer                        :: n2, nd, m1, niter, nframe, iframe, stat 
  real, dimension (:),   pointer :: mm, coord, ord, mnull, mbar, variance
  real                           :: o1, d1, scale, eps






  call sep_init ()

  call from_param( "n1", m1);  
  call from_param( "nframe", nframe, 1);  
  call from_param( "randscale", scale, 0.0);  
  call from_param( "niter", niter, m1);  
  call from_param( "eps", eps, 1.0);  
  call from_either( "o1", o1)
  call from_either( "d1", d1)

  allocate ( mm(m1), mnull(m1), mbar(m1), variance(m1) )

  call from_history (nd, n2)
  if (n2 /= 2) call erexit ('bad n2')
  allocate( coord(nd), ord(nd) )
  !--------------------------------------------------------------------------------
  ! Read time coordinates (coord) and picked RMS slowness (ord) from disk.
  !--------------------------------------------------------------------------------
  call sep_read (coord)
  call sep_read (ord)

  call to_history( "n1", m1);  call to_history( "n2", nframe+1)
  call to_history( "o1", o1);  call to_history( "d1", d1)
  call to_history( "n1", m1, "weight");  call to_history( "n2", 1, "weight")
  call to_history( "o1", o1, "weight");  call to_history( "d1", d1, "weight")
  if( exist_file("refl") ) then
    call to_history( "n1", m1, "refl");  call to_history( "n2", 1, "refl")
    call to_history( "o1", o1, "refl");  call to_history( "d1", d1, "refl")
  end if
  call sep_close()

  !--------------------------------------------------------------------------------
  ! 1-D nearest-neighbor binning (adjoint) maps ord to regular grid (mm).
  ! Write this to a file of reflections, "refl", for later use.
  !--------------------------------------------------------------------------------
  call bin1_init( m1, o1, d1, coord )
  stat = bin1_lop( .true., .false., mm, ord )

  mm=1.3   !!!!!delete this once you edit this!!!!!

  if( exist_file("refl") ) where(mm/=0) mm=1; call sep_write( mm, "refl")
  mnull=0.0;  mbar=mm;  variance=0.0

!--------------------------------------------------------------------------------
! Compute minimum variance solution here, store in array mbar.  
!--------------------------------------------------------------------------------
! !!!!!!!!!!!!!!!! COMPUTE MINIMUM VARIANCE SOLUTION HERE !!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------

  call sep_write(mbar)

!--------------------------------------------------------------------------------
! Compute multiple realizations (and variance) within this loop.
!--------------------------------------------------------------------------------
! !!!!!!!!!!!!!!!!!!!! COMPUTE MULTIPLE REALIZATIONS HERE !!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------

  do iframe=1,nframe
     call sep_write(mm)
  end do


!--------------------------------------------------------------------------------
! This is where you write "weight" to disk.  Default writes ones at pick 
! locations, zeros elsewhere.  You will probably want to remove the first line
! in your implementation, but keep the write statement.
!--------------------------------------------------------------------------------
! !!!!!!!!!!!!!!!!!!!!!!! COMPUTE RESIDUAL WEIGHT HERE !!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------

  call sep_write(mm, "weight")

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


  deallocate( mm, ord, coord, mbar, mnull, variance)

  call exit(0) 

end program InterpRMS


