module npef {                                   # Estimate non-stationary PEF
  use nhconest
  use npolydiv
  use polydiv
  use cgstep_mod
  use solver_mod
	use dottest
	use hconest
	use wilson
	use autocorr
	use helix
contains
  subroutine find_pef( dd, aa, rr, niter, eps, nh) {
    integer,          intent( in)    :: niter, nh      # number of iterations
    real,             intent( in)    :: eps	       # epsilon
    type( nfilter)                   :: aa             # estimated PEF output.
    type( nfilter),   intent( in)    :: rr             # roughening filter.
		type(filter)                     :: dummy
    real,    dimension(:), pointer   :: dd             # input data
    real,    dimension(:,:), pointer   :: rmov             # input data
    integer                          :: ip, ih, np, nr # filter lengths
    real, dimension (:), allocatable :: flt	       # np*na filter coefs
		real, dimension(:), allocatable  :: temp1,temp2
		real :: dot1(2),dot2(2),scale,tempr
		integer :: i1,ierr
    np = size( aa%hlx)			# data   length
    nr = np*nh
    allocate( flt( nr));flt=0.
    call nhconest_init( dd, aa, nh)
    call npolydiv_init( nr, rr)


		write(0,*) "dying 1"
		dd=dd-sum(dd)/size(dd)

#		 aa%hlx(1)%mis=>aa%mis
#		call hconest_init(dd,aa%hlx(1))
#		 call solver(hconest_lop,cgstep,aa%hlx(1)%flt,dat=-dd,niter=niter)

!    call solver_prec( nhconest_lop, cgstep, x= flt, dat= -dd, niter= niter,
!		      prec= npolydiv_lop, nprec= nr, eps= eps,verb=.true.)
    call solver( nhconest_lop, cgstep, x= flt, dat= -dd, niter= niter,verb=.true.)

    call   cgstep_close()
    call npolydiv_close()
    call nhconest_close()
    do ip = 1, np{
       do ih = 1, size( aa%hlx(ip)%lag)	{
          aa%hlx( ip)%flt( ih) = flt( (ip-1)*nh + ih)
			}
		}
		return;


		#TEST FOR BAD FILTER
		allocate(temp1(size(dd)),temp2(size(dd)))
		temp1=0	
		do i1=1,19
			temp1(floor(size(temp1)/19.*i1))=1.

		call npolydiv_init(size(dd),aa)
		ierr=npolydiv_lop(.false.,.false.,temp1,temp2)
		if(maxval(abs(temp2))>maxval(abs(temp1))*1.1) {
			write(0,*) "TRYING TO STABILIZE"
			scale=1.
			call srite("filter_warn",dd,size(temp1)*4)
			call srite("filter_warn",temp1,size(temp1)*4)
			call srite("filter_warn",temp2,size(temp2)*4)

			call stabilize(aa,temp1,temp2)
#			call wilson_init(maxval(aa%hlx(1)%lag)*10)
#			dummy%mis=>aa%mis
#			do ip=1,np{
#				tempr=1.
#				write(0,*) "here we go",ip,aa%hlx(ip)%flt
#				allocate(dummy%flt(size(aa%hlx(ip)%flt)))
#				allocate(dummy%flt(20));dummy%flt=0.
#				allocate(dummy%lag(20))
#				dummy%lag=(/1,2,3,4,5,6,7,8,9,10,191,192,193,194,195,196,197,198,199,200/)
##				dummy%lag=>aa%hlx(1)%lag
#				call wilson_factor(20,2.0,aa%hlx(ip),tempr,dummy)
#				deallocate(aa%hlx(ip)%flt)
#				deallocate(aa%hlx(ip)%lag)
#				allocate(aa%hlx(ip)%flt(20))
#				allocate(aa%hlx(ip)%lag(20))
#				aa%hlx(ip)%flt=dummy%flt*tempr
#				aa%hlx(ip)%lag=dummy%lag
#				write(0,*) "done we go",tempr,ip,aa%hlx(ip)%flt
#				deallocate(dummy%flt)
#			}
#			call wilson_close()
   		call npolydiv_close()
			call npolydiv_init(size(dd),aa)
			ierr=npolydiv_lop(.false.,.false.,temp1,temp2)
			call srite("filter_warn",temp2,size(temp2)*4)
			call npolydiv_close()
			call srite("filter_warn",dd,size(temp1)*4)
    	deallocate( flt,temp1,temp2)
		}
		else write(0,*) "GOOD F",maxval(temp2)
			
#			while(scale>.95 && maxval(abs(temp2)) > maxval(abs(temp1))*3.1){
#				scale=scale*.99
#				do ip=1,np
#					aa%hlx(ip)%flt=aa%hlx(ip)%flt*.99
#    		call npolydiv_close()
#				call npolydiv_init(size(dd),aa)
#				ierr=npolydiv_lop(.false.,.false.,temp1,temp2)
#			}
#			if(scale<.95) {
#				call auxputch("n1","d",size(dd),"filter_warn")
#				call auxputch("n2","d",1,"filter_warn")
##				call srite("filter_warn",dd,size(temp1)*4)
##				call srite("filter_warn",temp1,size(temp1)*4)
#				call srite("filter_warn",temp2,size(temp2)*4)
#				call erexit("Instable filter (filter_warn) separate will fail")
#			}
#			else{
#				write(0,*) "WARNING, FILTER COEF SCALED BY ", scale, " TO REACH STABILITY"
#			}
#		}
}

subroutine stabilize(filt,temp1,temp2){
type(nfilter) :: filt
type(filter)  :: autocor,temp_filt
real          :: a0,temp1(:),temp2(:)
integer       :: i1,i,ierr,nd


do i1=1,size(filt%hlx){
	if(1==i1-floor(i1/10.)*10.) write(0,*) i1
	write(0,*) "i1",i1,size(filt%hlx)
	write(0,*) "inpf",filt%hlx(i1)%flt
	write(0,*) "inpl",filt%hlx(i1)%lag
	call polydiv_init(size(temp1),filt%hlx(i1))
	ierr=polydiv_lop(.false.,.false.,temp1,temp2)
 	call polydiv_close()
	write(0,*) i1,maxval(abs(temp2))
	if(all(abs(temp2)<=1.01)) next  #stabil filter
	write(0,*) "failed",i1
	call allocatehelix(temp_filt,size(filt%hlx(i1)%lag))
	temp_filt%flt=filt%hlx(i1)%flt
	temp_filt%lag=filt%hlx(i1)%lag
	autocor=autocorr1 (filt%hlx(i1),  a0,1.)
#	write(0,*) "auto it",autocor%lag,autocor%flt
	call wilson_init(maxval(autocor%lag)*10)
	call wilson_factor(25,a0*1.01,autocor,a0,temp_filt)
	call wilson_close()
#	deallocate(filt%hlx(i1)%lag)
#	deallocate(filt%hlx(i1)%flt)
#	allocate(filt%hlx(i1)%lag(200))
#	allocate(filt%hlx(i1)%flt(200))
	call polydiv_init(size(temp1),temp_filt)
	ierr=polydiv_lop(.false.,.false.,temp1,temp2)
 	call polydiv_close()
	if(i1==435 || i1==678 || i1==840){
		write(0,*) "checking 435",abs(maxval(temp2))
		write(0,*) "filt in",filt%hlx(i1)%flt
		write(0,*) "filt out",temp_filt%flt
	}
	if(any(abs(temp2)>2.)) write(0,*) "BAD BAD BAD",maxval(abs(temp2))
	else filt%hlx(i1)%flt=temp_filt%flt

#	filt%hlx(i1)%lag=temp_filt%lag
	call deallocatehelix (temp_filt)
	call deallocatehelix (autocor)
	call polydiv_init(size(temp1),filt%hlx(i1))
	ierr=polydiv_lop(.false.,.false.,temp1,temp2)
 	call polydiv_close()
	if(any(temp2>1.2)){
		write(0,*) "BADF",i1,maxval(temp2)
		write(0,*) "FILT",filt%hlx(i1)%flt
		write(0,*) "LAGS",filt%hlx(i1)%lag
	}
}

}
}