module tomo_chain_mod
use tomo_wt_mod
use tomo_eps_mod
use dottest
use tomo_mask_mod
use tomo_op_mod
use tomo_reg_prec_mod
use tomo_eps_mod
use tomo_reg_prec_mod
implicit none
integer, private, save :: nfit,nstyle1
logical, private, save :: prec,lbfgs
real, pointer, private, save :: tempm(:),temp1(:),tempm2(:),tempx(:)
real, private, save :: eps_dampen

contains

logical function init_chains(gam0,reg0,dat0)
  real, pointer :: gam0(:),reg0(:),dat0(:),junk(:)
  integer :: ierr,i
  integer :: nm 
  init_chains=.false.

  call from_param("lbfgs",lbfgs,.false.)
  call from_param("prec",prec,.true.)

  nstyle1=size(reg0)
  nfit=wt_range()
  nm=reg_domain()
  allocate(tempm(nm))
  allocate(tempm2(nm),tempx(nm))
  allocate(dat0(nstyle1+nfit+nstyle1))
  dat0=0

  allocate(junk(nfit))
  ierr=wt_op(.false.,.false.,gam0,junk)
  dat0(1:nfit)=junk

  write(0,*)  "XXX",maxloc(abs(dat0)),maxval(junk)
  deallocate(junk)
  call from_param("eps_dampen",eps_dampen,1.)

  allocate(junk(size(reg0)))
 ! write(0,*) "BEFORE STUFF",sum(reg0)
  ierr=eps_op(.false.,.false.,reg0,junk)
  write(0,*) "ALP",maxval(abs(dat0)),maxval(abs(junk)),maxloc(abs(junk)),maxloc(abs(reg0))
  dat0(nfit+1:nfit+nstyle1)=junk
  write(0,*) "BEFORE STUFF",sum(reg0),sum(junk),maxval(abs(dat0))
  write(0,*)  "XXX",maxloc(abs(dat0)),maxval(junk),nfit+1,nfit+nstyle1
 ! call seperr("")
 deallocate(junk)


  allocate(temp1(max(size(gam0),nm)))
  deallocate(gam0,reg0)
  init_chains=.true.
end function

logical function lbfgs_val()
lbfgs_val=lbfgs
end function

integer function m_p(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)
m_p=0
!write(0,*) "IN MP",adj
!if (prec .and. .not. lbfgs) then
  call chain(mask_op,prec_op,adj,add,model,data,tempx)
!else
!  m_p=mask_op(adj,add,model,data)
!end if
end function

logical function chain_dot()
chain_dot=.false.
if(.not. my_dot(m_p,size(tempm),size(tempm))) then
   write(0,*) "failed m_p"
   return
end if
if(.not. my_dot(wt_op,wt_domain(),wt_range())) then
   write(0,*) "failed wt_op"
   return
end if
if(.not. my_dot(tomo_op,size(tempm),wt_domain())) then
   write(0,*) "failed tomo_op"
   return
end if
if(.not. my_dot(eps_op,size(tempm),size(tempm))) then
   write(0,*) "failed eps_op"
   return
end if
if(.not. my_dot(reg_op,size(tempm),size(tempm))) then
   write(0,*) "failed reg_op"
   return
end if
if(.not. my_dot(reg_term_op,size(tempm),size(tempm))) then
   write(0,*) "failed reg_term_op"
   return
end if
if(.not. my_dot(bound_op,size(tempm),size(tempm))) then
   write(0,*) "failed bound_op"
   return
end if
if(.not. my_dot(split_style_op,size(tempm),size(tempm))) then
   write(0,*) "failed split_style_op"
   return
end if
if(.not. my_dot(m_p,size(tempm),size(tempm))) then
   write(0,*) "failed m_p"
   return
end if
if(.not. my_dot(model_styling_op,size(tempm),size(tempm))) then
   write(0,*) "failed model_styling_op"
   return
end if

chain_dot=.true.
end function

logical function my_dot(oper,nmod,ndat)
  integer :: nmod,ndat
   interface 
     integer function oper(adj,add,model,data)
       logical, intent(in) :: adj,add
       real :: model(:),data(:)
     end function 
   end interface
   real :: dot1(2),dot2(2)

   my_dot=.true.
   call dot_test(oper,nmod,ndat,dot1,dot2)
   if(abs(dot1(1)/dot1(2)-1.)>.02) my_dot=.false.
   if(abs(dot2(1)/dot2(2)-1.)>.02) my_dot=.false.
   if(.not. my_dot) write(0,*) dot1,dot2
  

end function




integer function w_t(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)
w_t=0
!write(0,*) "W_T ",adj
call chain(wt_op,tomo_op,adj,add,model,data,temp1)
!write(0,*) "OUT W_T",adj
end function

integer function data_fitting_op(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)
data_fitting_op=0
!write(0,*) "IN DATA FITTING",adj
call chain(w_t,m_p,adj,add,model,data,tempm)
!write(0,*) "OUT DATA FITTING",adj
!data_fitting_op=tomo_op(adj,add,model,data);
end function






integer function bound_op(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)

bound_op=0
call adjnull(adj,add,model,data)
if(.not. adj) then
   data=data+model*eps_dampen
else 
  model=model+eps_dampen*data
end if
end function



integer function reg_term_op(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)

!write(0,*) "IN REG TERM",adj
reg_term_op=0
if(.not.prec) then
  call chain(eps_op,reg_op,adj,add,model,data,tempm)
else
   reg_term_op=eps_op(adj,add,model,data)
end if
!write(0,*) "OUT REG TERM",adj
end function



integer function split_style_op(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)
!write(0,*) "SPLIT STYLE",adj
!call array(reg_term_op,bound_op,adj,add,model,data(:nstyle1),data(nstyle1+1:))
call chain(reg_term_op,bound_op,adj,add,model,data,temp1)
!write(0,*) "OUT STLE",adj
end function



integer function model_styling_op(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)
!write(0,*) "MODEL STYLE",adj

call chain(split_style_op,m_p,adj,add,model,data,tempm2)
!write(0,*) "OUT MODEL STYLE",adj
model_styling_op=0
end function





integer function op(adj,add,model,data)
logical,intent(in) :: adj,add
real               :: model(:),data(:)

call array(data_fitting_op,model_styling_op,adj,add,model,data(1:nfit),&
  data(nfit+1:size(data)))
!op=data_fitting_op(adj,add,model,data(1:nfit))
op=0
end function







end module
