# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!err_mod
!
!Author: Robert Clapp, ESMB 463, 7230253
!
! Sun Jan 12 14:02:42 PST 2003
!
!
!Purpose:
!
! Handle the eps regularization
!
!To add a new error type:
!
! -add to check_pars
! -add to initializer
! -add to op
!
!
! Currently
!
! -scalar
! -gradient
!
module eps_reg_mod
  use adj_mod
  use param_check_mod
  use basic_util_mod
  use model_covar_mod
  implicit none
  character(len=128),private,save :: eps_method
  real, private,pointer,save :: eps_prec_scale(:),eps_reg_scale&
    &(:)
  real, private ,save :: eps_reg,grad_reg,eps_prec&
    &,grad_prec
  logical, private, save :: eps_diagonal_op,doing_prec&
    &,ignore_pad
  contains
  logical function init_eps_reg()
    logical :: no_rerurn,is_2d
    type(sep3d) :: tomo_s
    init_eps_reg=.false.
    if (.not. vel_tomo_structure(tomo_s)) then
      write(0,*) "trouble grabbing tomo structure"
      return
    end if
!check to see if parameters are different
    if (eps_method(1:6).eq."scalar" .or. eps_method(1:8).eq."gradient"&
      &) then
      if (.not. init_grad_eps(tomo_s)) then
        write(0,*) "trouble initializing gradient_eps"
        return
      end if
    else
      write(0,*) "Unknown resid method=",trim(eps_method)
      return
    end if
    init_eps_reg=.true.
  end function
  logical function get_model_styling_dat0(dat)
    real, pointer :: dat(:),junk(:)
    integer ::n
    get_model_styling_dat0=.false.
    if (.not. get_filt_d0(junk)) then
      write(0,*) "trouble getting filt d0"
      return
    end if
    allocate(dat(get_eps_reg_nrange(size(junk))))
    n=eps_reg_op(.false.,.false.,junk,dat)
    deallocate(junk)
    get_model_styling_dat0=.true.
  end function
  integer function get_eps_reg_nrange(n)
    integer :: n
    get_eps_reg_nrange=0
    if (eps_method(1:6).eq."scalar") then
      get_eps_reg_nrange=n
    else if (eps_method(1:4).eq."gradient") then
      get_eps_reg_nrange=n
    else
      write(0,*) "eps_model_method",trim(eps_method)
      call seperr("unknown eps_model_method")
    end if
  end function
  logical function init_grad_eps(struct)
    type(sep3d) :: struct
    integer :: n(3),i3,i1,i2,i,ierr
    real :: o(3),d(3),loc,a
    integer :: beg(3),end(3)
    init_grad_eps=.false.
    call set_n_o_d_3d(struct,n,o,d)
    allocate(eps_prec_scale(product(n)),stat=ierr)
    allocate(eps_reg_scale(product(n)),stat=ierr)
    if (ierr.ne.0) then
      write(0,*) "trouble allocating eps_scale"
      return
    end if
    call add_mem(size(eps_reg_scale)+size(eps_prec_scale))
    if (ignore_pad) then
      if (.not. return_padding(beg,end)) then
        write(0,*) "trouble getting padding info"
        return
      end if
    else
      beg=0
      end=0
    end if
    i=0
    do i3=1,n(3)
      do i2=1,n(2)
        do i1=1,n(1)
          loc=o(1)+d(1)*(i1-1)
          i= i+1
          if (i3 > beg(3) .and. i3 <= n(3)-end(3) .and. i2 > beg(2&
            &) .and. i2 <= n(2)-end(2) .and.i1 > beg(1) .and. i1 <= n&
            &(1)-end(1)) then
            eps_reg_scale(i)=eps_reg+loc*grad_reg
            eps_prec_scale(i)=eps_prec+loc*grad_prec
          else
            eps_reg_scale(i)=0.
            eps_prec_scale(i)=0.
          end if
        end do
      end do
    end do
    init_grad_eps=.true.
  end function
  integer function eps_reg_op(adj,add,model,data)
    logical,intent(in) :: adj,add
    real :: model(:),data(:),eps
    real, pointer :: eps_scale(:)
    call adjnull(adj,add,model,data)
!BY DEFAULT DIAGONAL 1
    if (eps_method(1:8).eq."gradient" .or. eps_method(1:6).eq."scalar"&
      &) then
      if (doing_prec) then
        eps_scale=>eps_prec_scale
      else
        eps_scale=>eps_reg_scale
      end if
      if (.not. adj) then
        data= data+model*eps_scale
      else
        model= model+data*eps_scale
      end if
    else
      write(0,*) "eps_reg_method",trim(eps_method)
      call seperr("Unknown eps_reg_method")
    end if
    eps_reg_op=0
  end function
  logical function model_eps_diagonal()
    model_eps_diagonal=eps_diagonal_op
  end function
  subroutine set_model_eps_prec(myp)
    logical :: myp
    doing_prec=myp
  end subroutine
  logical function check_eps_reg_params()
    character(len=128) :: tag
    logical :: stat(100)
    stat=.true.
    check_eps_reg_params=.false.
!#########################################
! IMPORTANT #
! #
! THIS WHERE WE DESCRIBE HOW TO #
! INITIALIZE/CHECK OBJECT PARAMS #
! #
!#########################################
    call from_param("eps_reg_method",eps_method,"scalar")
    eps_diagonal_op=.false.
    if (eps_method(1:6).eq."scalar" .or. eps_method(1:4).eq."grad")&
      & then
      eps_diagonal_op=.true.
      call from_param("eps_grad",grad_reg,0.)
      call from_param("eps_0",eps_reg,1.)
      eps_prec=eps_reg
      grad_prec=grad_reg
    else
      write(0,*) "    ERROR-Unknown eps_reg_method=",trim(eps_method)
      return
    end if
    call from_param("ignore_pad",ignore_pad,.false.)
    check_eps_reg_params=.true.
  end function
end module
