module non_lin_pick_rand_mod
use sep
use random
use pick_mod
use max_pick_mod
use sep3d_struct_mod
use io_nfilt
use solver_reg_mod
use nhelicon
use weight
use cgstep_mod
type(sep3d) :: d
real, private,pointer,save :: data(:,:),wt(:,:),range(:),eps(:),s(:,:,:)
integer, private,save :: niter,nl_iter
real, private,save :: pow


contains

logical function init_non_lin_pick(data_s,max_pick,weight,semb)
type(sep3d):: data_s
integer :: i
real,pointer :: max_pick(:,:),weight(:,:),semb(:,:,:)
real :: range_min,range_max,maxd,eps_min,eps_max,t1,t2

init_non_lin_pick=.false.

  d=data_s
  wt=>weight
  data=>max_pick
  s=>semb

  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))
  t2=max(nl_iter-1,1)
  do i=1,nl_iter
     t1=i-1.
     maxd=t1/t2
     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

integer function i_op(adj,add,model,data)
  logical,intent(in) :: adj,add
  real :: model(:),data(:)
  call adjnull(adj,add,model,data)
  if(.not. adj) data=data+model
  if(adj) model=model+data
  i_op=0
end function

logical function non_lin_pick(mod)
  real,pointer :: mod(:)
  real, pointer :: w(:),dat(:),rm0(:),r0(:),var(:)
  type(nfilter) :: filt
  real :: scale
  integer :: iter_nl,seed_it(2),seed
  non_lin_pick=.false.

  allocate(dat(size(data)))
  allocate(w(size(data)))
  allocate(mod(size(data)))
  allocate(var(size(data)))
  dat=reshape(data,(/size(data)/))
  w=reshape(wt,(/size(wt)/))
  mod=1.

  if(.not. calc_var(s,mod,var,d))&
   call seperr("trouble calculating variance")


  
  if(.not. read_nfilt("nfilt",filt)) then
    write(0,*) "trouble reading filter"
    return
  end if

  call from_param("seed",seed,-1)
  if(seed/=-1) then
    seed_it(1)=seed
    if(seed_it(1)==0) seed_it(1)=2000
    seed_it(2)=seed_it(1)/31328
    seed_it(1)=seed_it(1)-seed_it(2)*31328
    call random_seed(put=seed_it)
  end if
  call from_param("scale",scale,0.)

  allocate(r0(size(mod))) 
  allocate(rm0(size(mod))) 
  do i=1,size(mod)
     r0(i)=random_normal()*scale
  end do

  call nhelicon_init(filt)

  mod=0
  do  iter_nl=1,nl_iter
    write(0,*) "ITER",iter_nl," of ", nl_iter,minval(dat),maxval(dat),minval(mod),maxval(mod)
    call weight_init(w)
    dat=dat-1.
    rm0=r0*var
    call solver_reg(mod,dat,i_op,nhelicon_lop,cgstep,size(mod),niter,&
     eps(iter_nl),verb=.false.,wop=weight_lop,rm0=rm0)
    mod=mod+1.
    call cgstep_close()
    if(iter_nl/=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. calc_var(s,mod,var,d))&
       call seperr("trouble calculating variance")
    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 do

  deallocate(dat,w)
  non_lin_pick=.true.
end function





end module
