program Ultimatch_row
  
	use sep
	use cartesian
	use row_findnsmatch_mod
	use row_findmatch_mod

	implicit none

	real,pointer, dimension(:)			:: muls,prims,target,out1,out2
	type (filter)						:: match1,match2,laplac,laplac1,laplac2
	type (filter), pointer, dimension(:):: nmatch
	type (nfilter)						:: nsmatch1,nsmatch2,nsmatch

	integer :: ndim				! dimension of dataset
	integer :: mode             ! 1=single filter, 2=NS filtering
	integer :: niter
	integer :: ntimes
	integer :: rows				! =1 row operation =0 normal operation
	real    :: eps,rbar,mu
	logical :: twodim
	character (len=20) :: targetname
	integer, dimension(:), pointer :: nfilt  ! no. of coefficients in each filter
	integer, dimension(:), pointer :: n      ! data-size, dimension(ndim)
	real,    dimension(:), pointer :: d,o
	integer, dimension(:), pointer :: npatch ! no. of patches, dimension(ndim)
	integer, dimension(:), pointer :: nwind  ! overlapping patch-size
	integer, dimension(:), pointer :: idx    ! index counter (ndim)

	integer, dimension(:), pointer :: nh    ! length of each filter in the table
	integer, dimension(:), pointer :: pch   ! ns filter indexing
	integer, dimension(:), pointer :: psize ! non-overlapping patch-size

	integer :: i,stat,ilag,ipatch,idim,ipanels,npanels ! junk variables
	integer	:: nhmax
  
	call sep_init(SOURCE)

	ndim=sep_dimension()
write(0,*)'NDIM=',ndim
	allocate(n(ndim), o(ndim), d(ndim))
	
	!! parameters from history file (target)
	call from_history('n',n)
	call from_history('o',o)
	call from_history('d',d)

write(0,*)'n=',n,'o=',o,'d=',d

	!! if primaries requested as output, write out the history file parameters
	if (exist_file("primaries_out")) then
        call to_history('n1',n(1),"primaries_out")
        call to_history('d1',d(1),"primaries_out")
        call to_history('o1',o(1),"primaries_out")
		if (ndim>1) then
        	call to_history('n2',n(2),"primaries_out")
        	call to_history('d2',d(2),"primaries_out")
        	call to_history('o2',o(2),"primaries_out")
		end if
		if (ndim>2) then
        	call to_history('n3',n(3),"primaries_out")
        	call to_history('d3',d(3),"primaries_out")
        	call to_history('o3',o(3),"primaries_out")
		end if
	end if

	!! read matching parameters from the command line and set defaults
	call from_param('mode',mode,1)
	call from_param('2D',twodim,.false.)
	call from_param('ntimes',ntimes,1)
	call from_param('verb',verb,10)
	call from_param('row',rows,0)
	call from_param('mu',mu,1.)
	call from_param('eps',eps,1.) ! for the adaptive subtraction
write(0,*)'2D=',twodim,'ntimes=',ntimes,'verb=',verb,'row=',rows,'mu=',mu,'eps=',eps

	!! check number of dimensions and initialize correspondingly
	npanels=1
	if (twodim) then
		allocate(nfilt(2))
		nfilt=n(1:2)
write(0,*)'NFILT=',nfilt
		npanels=1
		if (ndim>2) npanels=n(3)
		deallocate(n)
		ndim=2
write(0,*)'NPANELS=',npanels
		allocate(n(2))
		n=nfilt
		deallocate(nfilt)
	end if

	allocate(nfilt(ndim),idx(ndim))
	nfilt=1; nfilt(1)=21				!! ????

	call from_param('nfilt',nfilt,nfilt)
write(0,*)'nfilt=',nfilt
	call from_param('niter',niter,product(nfilt)-1)
write(0,*)'niter=',niter

	!! Patching parameters
	if (mode /= 1) then
		allocate(npatch(ndim),nwind(ndim))
		npatch=1
		call from_param('npatch',npatch,npatch)
	end if

	allocate(muls(product(n)),prims(product(n)))
	allocate(target(product(n)),out1(product(n)),out2(product(n)))

	! Create TWO filters, one to match primaries and one to match multiples 
	call allocatehelix(match1,product(nfilt))
	call allocatehelix(match2,product(nfilt))
	do ilag=1,product(nfilt)
		call line2cart(nfilt,ilag,idx)
		match1%lag(ilag)=idx(1)-nfilt(1)/2-1
		match2%lag(ilag)=idx(1)-nfilt(1)/2-1
		do idim=2,ndim
			match1%lag(ilag)=match1%lag(ilag)+(idx(idim)-nfilt(idim)/2-1)*product(n(1:idim-1))
			match2%lag(ilag)=match2%lag(ilag)+(idx(idim)-nfilt(idim)/2-1)*product(n(1:idim-1))
		end do
	end do

	select case (mode)
	
		case (1)
			if (verb>0) write(0,*) '############################'
			if (verb>0) write(0,*) '# Single matched-filter...  '
			if (verb>0) write(0,*) '#---------------------------'
			if (verb>1) write(0,*) '#    nfilt  = ',nfilt
			if (verb>0) write(0,*) '############################'

		case (2)	
			if (verb>0) write(0,*) '#####################################'
			if (verb>0) write(0,*) '# Non-stationary matched-filtering...'
			if (verb>0) write(0,*) '#------------------------------------'
			if (verb>1) write(0,*) '#    ndim   = ',ndim
			if (verb>1) write(0,*) '#    n      = ',n
			if (verb>1) write(0,*) '#    npatch = ',npatch
			if (verb>1) write(0,*) '#    nfilt  = ',nfilt

			!! Create a non-stationary filter
			allocate(nh(product(npatch))) !! nh: length of each filter in table
			nh=product(nfilt)             !! Set them all constant
			call psize_init(n,npatch,psize) !! psize(ndim) is the patch-size
     
			if (verb>0) write(0,*) '#    patch-size = ',psize
			if (verb>0) write(0,*) '#####################################'
     
			call pch_init(n,npatch,psize,pch)  !! pch(ndata) is the index
			call nallocate(nsmatch1,nh,pch)
			call nallocate(nsmatch2,nh,pch)

			!! initialize the non-stationary filters
			do ipatch=1,product(npatch)
				nsmatch1%hlx(ipatch)%lag=match1%lag
				nsmatch2%hlx(ipatch)%lag=match2%lag
			end do

			!! compute the helix version of the laplacian regularizer
			nhmax=size(nsmatch1%hlx(1)%lag)+size(nsmatch2%hlx(1)%lag)
			call laplac1_init(laplac,nhmax,npatch) !! laplac is the regularizer

		case default
			call seperr('mode not implemented')
	end select	

	!! compute the filters panel by panel
	panels:do ipanels=1,npanels

		if (verb>1) write(0,*) '###################################'
		if (verb>1) write(0,*) '   Panel = ',ipanels,' out of ',npanels
		if (verb>1) write(0,*) '-----------------------------------'

		out1=0.
		out2=0.
     
		! Read in data
        call sep_read(target)
		call sep_read(muls,'multiples_in')
		call sep_read(prims,'primaries_in')

		if ((maxval(muls).eq.minval(muls)).or.(maxval(prims).eq.minval(prims))) then
			!call sep_write(out1)
			cycle
		end if

		select case(mode)

			case (1) !! one filter

				if (verb>2) write(0,*) ' -----------------------------------'
				if (verb>2) write(0,*) ' Calling findmatch with no weight...'
				if (verb>2) write(0,*) ' -----------------------------------'

				call row_findmatch(muls,prims,target,match1,match2,niter,mu)

				call dumpfilter(match1)
				call dumpfilter(match2)

				!! apply the first filter to the data to get the matched multiples
				call applymatch(muls,out1,match1)

				!! if requested, apply the second filter to the data to get 
				!!	the matched primaries
				if (exist_file("primaries_out")) then
					call applymatch(prims,out2,match2) 
				end if

			case (2) !! NS filters
           
				if (verb>2) write(0,*) ' -------------------------------------'
				if (verb>2) write(0,*) ' Calling row findnsmatch with no weight...'
				if (verb>2) write(0,*) ' -------------------------------------'
           
				!! the 1 is for regularized solver.
				call row_findnsmatch(muls,prims,target,nsmatch1,nsmatch2,&
					& laplac,1,niter,eps,mu)
           
				!! apply the first filter to the data to get the matched multiples
				call applynsmatch(muls,out1,nsmatch1) 
        
				!! if requested, apply the second filter to the data to get 
				!!	the matched primaries
				if (exist_file("primaries_out")) then
					call applynsmatch(prims,out2,nsmatch2) 
				end if 

			case default
				call seperr('mode not implemented')
		end select

		out1=mu*out1
		call sep_write(out1)
			
		if (exist_file("primaries_out")) then
			call srite("primaries_out",out2,4*size(out2))
		end if
	end do panels
  
	write(0,*) '---- The End ----'
  
	call sep_close()
  
	contains
     
	!! Calculate non-overlapping patch-size given a certain number of patches
	subroutine psize_init(n,npatch,psize)
		integer, dimension(:), intent(in)  :: n,npatch
		integer, dimension(:), pointer     :: psize
		allocate(psize(ndim))
		do i=1,ndim
			psize(i)=1+(n(i)-1)/npatch(i)
		end do
	end subroutine psize_init

	!! Divides data-space between rectangular patches
	subroutine pch_init(n,npatch,psize,pch)
		integer, dimension(:), intent(in)  :: n,npatch,psize
		integer, dimension(:), pointer     :: pch,ii
		integer :: i,j

		allocate(pch(product(n)),ii(size(n)))
		do i=1,product(n)
			pch(i)=0
			call line2cart(n,i,ii)
			do j=1,size(n)
				ii(j)=(1+(ii(j)-1)/psize(j))
			end do
			call cart2line(npatch,ii,pch(i))
		end do
		deallocate(ii)
	end subroutine pch_init
  
	!! Create a roughener to dampen NS filter estimation
	subroutine laplac_init(laplac,nsfilt,npatch)
		type  (filter)        :: laplac
		type (nfilter)        :: nsfilt
		integer, dimension(:) :: npatch
		real    :: e1,e2
		integer :: nhmax

		nhmax=size(nsfilt%hlx(1)%lag)

		write(0,*) 'size',size(npatch)
		select case(size(npatch))
			case (1)
				call allocatehelix(laplac,1)  
				laplac%lag(1)=nhmax;   laplac%flt(1)=-.99
			case (2)
				!! 1-D smoothing only

				! crap 2-D smoothing
				call allocatehelix(laplac,3)  
				laplac%lag=(/ nhmax, npatch(1)*nhmax, nhmax*(npatch(1)+1)  /)
				e1=0.2; e2=0.1
				laplac%flt=(/ e1-1.,      e2-1.,     (1.-e1)*(1.-e2)  /)

				if (verb>5) then
					write(0,*) '  Roughening filter...'
					call dumpfilter(laplac)
				end if

			case default
				write(0,*) 'Testmatch'
				call seperr('dimensionality not implemented')
		end select

	end subroutine laplac_init

	!! Create a roughener to dampen NS filter estimation
	subroutine laplac1_init(laplac,nhmax,npatch)
		type  (filter)			:: laplac
		integer					:: nhmax
		integer, dimension(:)	:: npatch
		real    				:: e1,e2

		write(0,*) 'size',size(npatch)
		select case(size(npatch))
			case (1)
				call allocatehelix(laplac,1)  
				laplac%lag(1)=nhmax;   laplac%flt(1)=-.99
			case (2)
				!! 1-D smoothing only

				! crap 2-D smoothing
				call allocatehelix(laplac,3)  
				laplac%lag=(/ nhmax, npatch(1)*nhmax, nhmax*(npatch(1)+1)  /)
				e1=0.2; e2=0.1
				laplac%flt=(/ e1-1.,      e2-1.,     (1.-e1)*(1.-e2)  /)

				if (verb>5) then
					write(0,*) '  Roughening filter...'
					call dumpfilter(laplac)
				end if

			case default
				write(0,*) 'Testmatch'
				call seperr('dimensionality not implemented')
			end select

		end subroutine laplac1_init

		subroutine row_laplac_init (laplac1,laplac2,nsmatch1,nsmatch2,npatch)
			type  (filter)        :: laplac1
			type  (filter)        :: laplac2
			type (nfilter)        :: nsmatch1
			type (nfilter)        :: nsmatch2
			integer, dimension(:) :: npatch

			call laplac_init(laplac1,nsmatch1,npatch)
			call laplac_init(laplac2,nsmatch2,npatch)
	
		end subroutine row_laplac_init

end program Ultimatch_row
