module compute
	use likmod
	use altmod
	implicit none
	
	contains
	
	subroutine setfs	! draw ys and store importance sampling, null and alternative densities
		real	:: y(capn),eps(capn),th(nth),h,hs(nimp),lfimp
		integer	:: l,i
		allocate(f1(nsim),LR(n0,nsim))
		do i=1,n0
			thlist0(:,i)=getth([(i-1.0)/(n0-1)])
		enddo
		call mdisp(thlist0)
!$omp parallel do private(eps,th,y,i,hs,lfimp,h)
		do l=1,nsim
			call rnnoa(eps)
			th=thlistimp(:,mod(l,nimp)+1)
			call setyfromeps(th,y,eps)
			do i=1,nimp
				th=thlistimp(:,i)
				hs(i)=getKalman_invar(th,y)
			enddo
			lfimp=maxval(hs)
			hs=exp(hs-lfimp)
			lfimp=lfimp+log(sum(hs)/nimp)
			h=0
			do i=1,n1
				h=h+exp(getKalman_invar1(thlist1(:,i),y)-lfimp)
			enddo
			h=h/n1
			f1(l)=h
			do i=1,n0
				LR(i,l)=exp(getKalman_invar(thlist0(:,i),y)-lfimp)
			enddo
		enddo
		call mdisp(sum(LR,dim=2)/nsim)
	end subroutine
	
	subroutine setlam			! find non-negative weights so that maximum rejection probability is close to level
		real	:: rp(n0)
		integer	:: i,l
		
		lam=1
		do i=1,1000
			rp=0
			do l=1,nsim
				if(sum(lam*LR(:,l))<f1(l)) rp=rp+LR(:,l)
			enddo
			rp=rp/nsim
			lam=lam*exp(rp-level)
			if(mod(i,500)==0) then
				call mdisp(rp.cvr.lam)
			endif
		enddo
	end subroutine
	
	subroutine loaddata
		real, allocatable	:: mdata(:)
		integer	:: s,j,T,Tdata
		mdata=loadvec("c:/dropbox/mystuff/beyondltu/USUK_RER_long.txt")
		T=size(mdata)
		Tdata=T
		print *,"lenght of data in years:", Tdata
		j=1
		do s=1,T
			if(s*capN>=j*T) then
				ydata(j)=mdata(s)
				j=j+1
			endif
		enddo
		call mdisp(mdata)
		call mdisp(ydata)
	end subroutine

	subroutine testydata			! compute test on actual data
		real	:: h0,h1,h1s(n1)
		integer	:: i
		
		h0=0
		do i=1,n0
			h0=h0+lam(i)*exp(getKalman_invar(thlist0(:,i),ydata))
		enddo
		
		h1=0
		do i=1,n1
			h1=h1+exp(getKalman_invar1(thlist1(:,i),ydata))
			h1s(i)=exp(getKalman_invar1(thlist1(:,i),ydata)-log(h0))
		enddo
		h1=h1/n1
		call mdisp(sort(h1s))
		print *,h1/h0
	
	end subroutine
		
end module