module basicmodule
	use myfuncs
	use dotops
	use ML

	integer, parameter	:: nsim=50000
	integer, parameter	:: capN=50, Tstar=1000
	integer, parameter	:: n1=500,p1=3, nth1=2*p1-1
	integer, parameter	:: p=1, nth=2*p-1 
	real, parameter		:: hub=200,hlb=0.5
	
	integer, parameter	:: n0=200
	integer				:: nimp
	real, allocatable	:: thlistimp(:,:),thlist1(:,:)
	real				:: thlist0(nth,n0),lam(n0)
	real, allocatable	:: f1(:),LR(:,:)
	
	real				:: level=0.05
	real				:: ydata(capn)
	
end module
	
module likmod
	use basicmodule
	use polymod
	implicit none
	
	contains
	
	function getcfromh(h) result(c)
		real		:: h(:)
		complex	:: c(size(h))
		real		:: h1,h2
		complex	:: d
		integer	:: i
		do i=1,size(h)/2
			h1=h(2*i-1)
			h2=h(2*i)
			d=sqrt(cmplx(h2**2 - h1**2))
			c(2*i-1:2*i)=[h2-d,h2+d]
		enddo
		if(mod(size(h),2)==1) then
			c(size(h))=h(size(h))
		endif
	end function
	
	function getstatvar(th) result(val)
		use lincg_int
		use iercd_int
		real			:: th(nth),val(p,p)
		complex	:: lam(p), Q(p,p), Qinv(p,p)
		complex	:: a(p**2), b(p**2)
		integer	:: i,j
		lam=-getcfromh(th(1:p))
		do j=1,p
			Q(1,j)=1.0
			do i=2,p
				Q(i,j)=lam(j)*Q(i-1,j)
			enddo
			Q(:,j)=Q(:,j)/norm2(abs(Q(:,j)))
		enddo
		lam=1+lam/Tstar
		call lincg(Q,Qinv)
		if(IERCD()>0) call mdisp(real(lam).cvr.aimag(lam))
		a=1-[spread(lam,1,p)]*[spread(lam,2,p)]
		b=[spread(Qinv(:,p),dim=2,ncopies=p)*spread(Qinv(:,p),dim=1,ncopies=p)]
		val=reshape([matmul(Q,matmul(reshape(b/a,[p,p]),transpose(Q)))],[p,p])
	end function

	subroutine setstatespace(th,phi,theta,statevar)
		real			:: th(nth),Phi(p,p),theta(p)
		real			:: poly(p+1)
		real			:: statevar(p,p)
		integer		:: j
		
		Phi=0
		do j=2,p
			Phi(j-1,j)=1
		enddo
		poly=realpolyfromroots(-getcfromh(th(1:p)))
		Phi(p,:)=-poly(1:p)
		if(p>1) then
			theta=realpolyfromroots(-getcfromh(th(p+1:)))
		else
			theta=1
		endif
		Phi=eye(p)+phi/Tstar
		statevar=getstatvar(th)
	end subroutine
	
	function getKalman_invar(th,y) result(val)
		real			:: th(nth),y(capN)
		real			:: val
		real			:: Phi(p,p),theta(p)
		real			:: state(p,2), statevar(p,p), vst(p), cv, e(2)
		real			:: ssum(4)
		integer			:: s,j

		call setstatespace(th,phi,theta,statevar)
		state=0
		ssum=0
		j=1
		do s=1,Tstar
			state=matmul(Phi,state)
			statevar=matmul(Phi,matmul(statevar,transpose(Phi)))
			statevar(p,p)=statevar(p,p)+1.0
			if(s*capN>=j*Tstar) then
				vst=matmul(statevar,theta)
				cv=1.0/sum(vst*theta)
				e=[y(j),1.0]-matmul(theta,state)		
				ssum(1:3)=ssum(1:3)+cv*[e(1)**2,e(2)**2,e(1)*e(2)]
				ssum(4)=ssum(4)+log(cv)
				state=state+matmul(reshape(vst,[p,1]),reshape(cv*e,[1,2]))
				statevar=statevar-matmul(reshape(vst*cv,[p,1]),reshape(vst,[1,p]))
				j=j+1
			endif
		enddo
		val=.5*ssum(4)-.5*log(ssum(2))-0.5*(capN-1)*log(ssum(1)-ssum(3)**2/ssum(2))
	end function
	
	subroutine setyfromeps(th,y,eps)
		real			:: th(nth),y(capN),eps(capN)
		real			:: val
		real			:: Phi(p,p),theta(p)
		real			:: state(p), statevar(p,p), vst(p), cv
		integer			:: s,j
		call setstatespace(th,phi,theta,statevar)
		state=0
		j=1
		do s=1,Tstar
			state=matmul(Phi,state)
			statevar=matmul(Phi,matmul(statevar,transpose(Phi)))
			statevar(p,p)=statevar(p,p)+1.0
			if(s*capN>=j*Tstar) then
				vst=matmul(statevar,theta)
				cv=1.0/sum(vst*theta)
				y(j)=sum(theta*state)+eps(j)/sqrt(cv)
				state=state+vst*(sqrt(cv)*eps(j))
				statevar=statevar-matmul(reshape(vst*cv,[p,1]),reshape(vst,[1,p]))
				j=j+1
			endif
		enddo
	end subroutine

	function rootdist0(th) result(val)
		real		:: th(nth),val
		complex	:: dlist(p-1),clist(p)
		integer		:: i,j
		if(p==1) then
			val=1E100
			return
		endif
		clist=getcfromh(th(1:p))
		val=1E100
		do i=1,p
			do j=1,i-1
				val=min(abs(clist(i)-clist(j)),val)
			enddo
		enddo
	end function	

	function rootdist(th) result(val)
		real		:: th(nth),val
		complex	:: clist(nth)
		integer		:: i,j
		if(p==1) then
			val=1E100
			return
		endif
		clist(1:p)=getcfromh(th(1:p))
		clist(p+1:)=getcfromh(th(p+1:))
		val=1E100
		do i=1,nth
			do j=1,i-1
				val=min(abs(clist(i)-clist(j)),val)
			enddo
		enddo
	end function	

	function getsmoothprior(th) result(val)
		real		:: th(nth),val
		complex		:: glist(p-1),clist(p)
		integer		:: i,j	
		clist=getcfromh(th(1:p))
		if(p>1) glist=getcfromh(th(p+1:))
		val=0
		do i=1,p
			do j=1,p
				val=val+(clist(i)+clist(j))**-3
			enddo
		enddo
		do i=1,p
			do j=1,p-1
				val=val-2*(clist(i)+glist(j))**-3
			enddo
		enddo
		do i=1,p-1
			do j=1,p-1
				val=val+(glist(i)+glist(j))**-3
			enddo
		enddo
		val=val*8*Pi
	end function
	
	function getth(x) result(val)
		real	:: x(nth), val(nth)
		val=hlb+(hub-hlb)*x
	end function

end module