module analyticalIS
	use likmod
	use iomod
	implicit none

	integer, parameter	:: nconst=0, nz=capN+nth

	contains
		
	subroutine setanaIS
		use bconf_int
		use nnlpf_int
		use nnlpg_int
		real		:: th(nth)
		external	:: getminusISratio, aconstr_func, aconstr_grad
		real		:: zlb(nz),zub(nz),z(nz),zguess(nz),zo(nz)
		real		:: r0,r1
		integer, parameter	:: ntries=50
		integer		:: i
		
		if(allocated(thlistimp)) deallocate(thlistimp)
		allocate(thlistimp(nth,1))
		call rnun(thlistimp(:,1))
		thlistimp(:,1)=getth(thlistimp(:,1))
		nimp=1
		zlb=0.001; zub=1-zlb
		call erset(0,0,0)
		i=0
		do
			call rnun(zguess)
			call nnlpg(aconstr_func,aconstr_grad,nconst,0,0,zlb,zub,x=z,xguess=zguess,maxitn=100,iprint=0)
			if(.not. isvalid(z)) cycle
			i=i+1
			th=getth(z(1:nth))
			if(rootdist(th)<0.5) cycle
			if( getISratio(z)>5 .and. all(isfinite(th))) then
				call mdisp([real(nimp),real(i),getISratio(z)])
				call mdisp(th)
				call mdisp(getYfromz(z))
				call plotspecdens(th)		! plot the implied spectral density in matlab
				zo=z
				r0=getISratio(zo)
				nimp=nimp+1
				thlistimp=reshape([thlistimp,th],[nth,nimp])
				r1=getISratio(zo)
				print *,"before",r0,"after",r1
				i=0
				cycle
			endif
			if(i>ntries) exit
		enddo
		print *,"done with IS determination"
		call mdisp(thlistimp)
		call printtime
	end subroutine
	
	function getYfromz(z) result(y)
		real	:: th(nth),z(nz), y(capn),eps(capn)
		th=getth(z(1:nth))
		eps=z(nth+1:)
		call correps(eps)
		call setyfromeps(th,y,eps)
		
	contains
	
		subroutine correps(eps)
			real	:: eps(capn),e2,co
			eps=gausscdfinv(eps)
			e2=sum(eps**2)
			co=capn+3*sqrt(2.0*capn)
			if(e2>co) eps=eps/sqrt(e2/co)
		end subroutine

	end function
		
	function getISratio(z) result(val)
		real	:: z(nz), val
		real	:: th(nth), Y(capn), rs(nimp)
		integer	:: i
		
		th=getth(z(1:nth))
		Y=getYfromz(z)
		val=0
!$omp parallel do
		do i=1,nimp
			rs(i)=getKalman_invar(thlistimp(:,i),Y)
		enddo
		rs=rs-getKalman_invar(th,y)
		val=1.0/(sum(exp(rs))/nimp)
	end function
	
	function isvalid(z) result(val)
		real	:: z(nz), c
		logical	:: val, ierr
		integer	:: i
		external	:: getminusISratio, aconstr_func, aconstr_grad

		do i=1,nconst
			call aconstr_func(z,i,c,ierr)
			if(c<-.01) then
				val=.false.
				return
			endif
		enddo
		val=.true.

	end function
end module
	
subroutine getminusISratio(n,z,val)
	use analyticalIS
	implicit none
	
	integer	:: n
	real	:: z(n), val
	val=-getISratio(z)
	if(val/=val) val=huge(1.0)
end subroutine
	
subroutine aconstr_func(z,iact,val,ierr)
	use analyticalIS
	implicit none
	integer	:: iact
	real	:: z(nz), val
	logical	:: ierr
	real	:: Y(capn),th(nth),rd
	
	rd=rootdist(getth(z(1:nth)))
	if(rd<0.5) then
		ierr=.true.
		val=1E150
		return
	endif

	select case(iact)
		case(0)
			val=-getISratio(z)
		case(1)
			val=rootdist(getth(z(1:nth)))-1.0
		case(2)
			val=.05-getsmoothprior(getth(z(1:nth)))
	end select
		
	ierr=.not. isfinite(val)
	if(ierr) val=1E50
end subroutine	
	
subroutine aconstr_grad(z,iact,val)
	use analyticalIS
	implicit none
	integer	:: iact
	real	:: z(nz), val(nz)
	logical	:: ierr
	integer	:: i
	real	:: zc(nz), v, vvec(nz)
	real, parameter	:: eps=epsilon(1.0)**(0.33)
	do i=1,nz
		zc=z
		zc(i)=zc(i)+eps
		call aconstr_func(zc,iact,vvec(i),ierr)
		zc(i)=zc(i)-2*eps
		call aconstr_func(zc,iact,v,ierr)
		vvec(i)=vvec(i)-v
		zc(i)=zc(i)+eps
	enddo
	val=vvec/(2*eps)
	if(any(.not. isfinite(val))) val=0
end subroutine	

	