# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module non_lin_pick_mod
  use sep
  use max_pick_mod
  use sep3d_struct_mod
  use igrad2
  use solver_mod
  use cgstep_mod
  implicit none
  type(sep3d) :: d
  real, private,pointer,save :: data(:,:),wt(:,:),range(:),eps(:)
  integer, private,save :: niter,nl_iter
  real, private,save :: pow
  contains
  logical function init_non_lin_pick(data_s,max_pick,weight)
    type(sep3d):: data_s
    integer :: i
    real,pointer :: max_pick(:,:),weight(:,:)
    real :: range_min,range_max,maxd,eps_min,eps_max
    init_non_lin_pick=.false.
    d=data_s
    wt=>weight
    data=>max_pick
    call from_param("niter",niter,200)
    call from_param("pow",pow,1.)
    call from_param("nl_iter",nl_iter,1)
    maxd=max(abs(1.-d%o(3)),abs(1.-d%o(3)-(d%n(3)-1)*d%d(3)))
    call from_param("range_max",range_max,maxd/3.)
    call from_param("range_max",range_min,maxd/6.)
    call from_param("eps_max",eps_max,1.)
    call from_param("eps_min",eps_min,1.)
    allocate(range(nl_iter),eps(nl_iter))
    do i=1,nl_iter
      maxd=real(i-1)/real(nl_iter-1)
      eps(i)=(1.-maxd)*eps_max+maxd*eps_min
      range(i)=(1.-maxd)*range_max+maxd*range_min
    end do
    init_non_lin_pick=.true.
  end function
  logical function non_lin_pick(mod)
    real,pointer :: mod(:)
    real, pointer :: w(:),dat(:)
    integer :: iter_nl
    non_lin_pick=.false.
    allocate(dat(size(data)))
    allocate(w(size(data)))
    allocate(mod(size(data)))
    dat=reshape(data,(/size(data)/))
    w=reshape(wt,(/size(wt)/))
    mod=1.
    call igrad2_init(d%n(1),d%n(2))
    do iter_nl=1,nl_iter
      write(0,*) "ITER",iter_nl," of ", nl_iter
      call srite("aa.H",dat,size(dat)*4)
      call solver_reg(i_op,cgstep,igrad2_lop,size(mod)*2,mod,dat,niter&
        &,eps(iter_nl),verb=.false.,x0=mod,wt=w)
      call cgstep_close()
      if (iter_nl.ne.nl_iter) then
        if (.not. reset_bounds(mod,range(iter_nl))) then
          write(0,*) "trouble reseting bounds after ",iter_nl,"&
            & iterations"
          return
        end if
        if (.not. max_pick(data,wt)) then
          write(0,*) "trouble picking max after ",iter_nl,"&
            & iterations"
          return
        end if
        dat=reshape(data,(/size(data)/))
        w=reshape(wt,(/size(wt)/))
      end if
    end do
    deallocate(dat,w)
    non_lin_pick=.true.
  end function
  integer function i_op(adj,add,model,data)
    logical :: adj,add
    real :: model(:),data(:)
    call adjnull(adj,add,model,data)
    if (.not. adj) then
      data= data+model
    else
      model= model+data
    end if
    i_op=0
  end function
end module
