module globals
	use myfuncs
	use dotops
	use ML
   	integer, parameter	:: nsim=250000				! number of monte Carlo draws
	integer, parameter	:: mkSig_T=5000				! number of divisions of unit interval for Riemann integration of Sigmas
	integer				:: q=12, qmax=12			! default values for q
	real				:: level=0.1				! default value of alpha
	real, parameter		:: lengthfudge=0.01			! epsilon of epsilon ALFD
	logical				:: xgridflag=.true.			! flag to generate extended grid for size checking
	logical				:: cleanLamflag=.false.		! flag whether lambda vectors are rounded down to zero if small 
	logical				:: saveSigflag=.false.		! if true, save all Sigma matrices (input to empirical work), and corresponding list of thetas
	

	real, parameter		:: levellist(2)=[.1,.33]
	real, parameter		:: clist(12)=[0.00, exp(-3+7.0*([(i,i=0,10)]/10.0))]
	real, parameter		:: dlist(8)=[(real(i)/10,i=-4,10,2)]
	real, parameter		:: blist(12)=[0.0,exp(-5+10.0*([(i,i=0,10)]/10.0))]

	real, allocatable	:: rlist(:)
	real, parameter		:: rlist0(1)=[0.5]
	real, parameter		:: rlistx(16)=[.075, .1, .15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4,1.6,1.8]

	real, parameter		:: clistx(33)=[0.00, exp(-3+7.0*([(i,i=-1,30)]/20.0))]
	real, parameter		:: dlistx(15)=[(real(i)/10,i=-4,10,1)]
	real, parameter		:: blistx(28)=[0.00, exp(-5+10.0*([(i,i=-1,25)]/20.0))]

	real, parameter		:: niceblist(size(blist))=((8*pi)**2)*[0.0,[.5,1.0,2.0,5.0,10.0,20.0,50.0,100.0,200.0,500.0,1000.0]**(-2)]
	logical				:: useniceblist=.false.							! flag whether to use niceblist (for Table 2), or standard blist
	integer, parameter	:: n=size(dlist)+size(dlist)*((size(blist)-1)+(size(clist)-1))+(size(clist)-1)*(size(blist)-1)	! number of theta points in ALFD
	integer				:: nws=size(dlist)								! number of theta points that appear in length calculation
	integer,parameter	:: nx=n+size(clistx)*size(dlistx)*size(blistx)	! number of theta points for coverage check
	
	real, allocatable	:: Sigs(:,:,:), SigXinvs(:,:,:),Siginvs(:,:,:)	
	real				:: SigXdets(nx),Sigdets(nx),thlist(3,nx)
	real, allocatable	:: w_prior(:), lengths0(:)

	real, allocatable	:: Ys(:,:)
	real				:: densavWs(nsim)
	logical				:: bayesrej(nsim)
	real, allocatable	:: densWs(:,:),densgs(:,:),densXs(:,:),ISw(:,:)

end module



module compute
	use globals
    implicit none

	contains

	
	subroutine mkSigs(q,ir)
		implicit none
		integer		:: q, ir
		
		integer, parameter		:: nbase=size(dlist)*size(clist)+size(dlistx)*size(clistx)
		integer		:: ic,id,ib, lc,ibase,i
		
		real, save, allocatable		:: Sigx(:,:,:)
		real, save		:: thbaselist(2,nbase)
		logical, save	:: SigxComputed=.false.
		logical, save	:: Sigs_alreadysaved=.false.
		integer			:: ind(q+1)
		
		if(SigxComputed) deallocate(Sigs, SigXinvs,Siginvs,Ys)
		allocate(Sigs(q+1,q+1,nx), SigXinvs(q,q,nx),Siginvs(q+1,q+1,nx),Ys(q+1,nsim))

		if(.not.saveSigflag) Sigs_alreadysaved=.true.

		if(.not.SigxComputed) then
			allocate(Sigx(qmax+size(rlist),qmax+size(rlist),nbase))
			do id=1,size(dlist)
				call mkSig(dlist(id),0.0,Sigx(:,:,id))
				thbaselist(:,id)=[dlist(id),0.0]
			enddo
			ibase=size(dlist)
!$omp parallel do private(ic)				
			do id=1,size(dlist)
				do ic=2,size(clist)
					call mkSig(dlist(id),clist(ic),Sigx(:,:,ibase+(size(clist)-1)*(id-1)+ic-1))
					thbaselist(:,ibase+(size(clist)-1)*(id-1)+ic-1)=[dlist(id),clist(ic)]
				enddo
			enddo
			ibase=ibase+size(dlist)*(size(clist)-1)
			print *,ibase
			if(xgridflag) then
!$omp parallel do private(id)				
				do ic=1,size(clistx)
					do id=1,size(dlistx)
						call mkSig(dlistx(id),clistx(ic),Sigx(:,:,ibase+size(dlistx)*(ic-1)+id))
						thbaselist(:,ibase+size(dlistx)*(ic-1)+id)=[dlistx(id),clistx(ic)]
					enddo
				enddo
			endif
			call VSprint("Sigx computation done")
			SigxComputed=.true.
		endif
				
		ind=[(i,i=1,q),qmax+ir]
		lc=1
		do id=1,size(dlist)
			call setSigetc(id,blist(1),clist(1),dlist(id))	
		enddo
		do id=1,size(dlist)
			do ib=2,size(blist)
				if(id==size(dlist).and.useniceblist) then
					call setSigetc(id,niceblist(ib),0.0,1.0)
				else
					call setSigetc(id,blist(ib),clist(1),dlist(id))	
				endif
			enddo
		enddo
		ibase=size(dlist)+1
		do id=1,size(dlist)
			do ic=2,size(clist)
				call setSigetc(ibase,blist(1),clist(ic),dlist(id))
				ibase=ibase+1
			enddo
		enddo
		ibase=size(dlist)+(size(dlist)-1)*(size(clist)-1)+1
		id=size(dlist)
		do ic=2,size(clist)
			do ib=2,size(blist)
				if(clist(ic)==0.0 .and. blist(ib)==0.0) cycle
				call setSigetc(ibase,blist(ib),clist(ic),dlist(id))
			enddo
			ibase=ibase+1
		enddo
		if(.not.Sigs_alreadysaved) then
			call storeML(transpose([w_prior,(0.0,i=1,n-nws)].cud.thlist(:,1:n)),'wthlist')
			call execinML("save('c:/out/wthlist.asc','wthlist','-ascii')")
			Sigs_alreadysaved=.true.
		endif
		if(xgridflag) then
			do ic=1,size(clistx)
				do id=1,size(dlistx)
					do ib=1,size(blistx)
						call setSigetc(ibase,blistx(ib),clistx(ic),dlistx(id))
					enddo
					ibase=ibase+1
				enddo
			enddo
		endif
	contains 
		subroutine setSigetc(ibase,b,c,d)
			implicit none
			integer	:: ibase
			real	:: b,c,d
			real	:: Sig(q+1,q+1),cSig(qmax+size(rlist),qmax+size(rlist)),bfac

			if(sum(abs(thbaselist(:,ibase)-[d,c]))>1E-5) then
				print *,'thbaselist error'
				stop
			endif
			cSig=Sigx(:,:,ibase)
			bfac=2*pi*b*((8*pi)**2+c**2)**(-d)
			do i=1,qmax
				cSig(i,i)=cSig(i,i)+bfac
			enddo
			do i=1,size(rlist)
				cSig(qmax+i,qmax+i)=cSig(qmax+i,qmax+i)+bfac*(1/rlist(i)+1)
			enddo
			cSig=cSig/cSig(1,1)
			if(.not.Sigs_alreadysaved) then
				call storeML(cSig,'Sig')
				call VSprint("save('c:/out/Sig"//itostring(lc)//".asc','Sig','-ascii')")
				call execinML("save('c:/out/Sig"//itostring(lc)//".asc','Sig','-ascii')")
			endif
			Sig=cSig(ind,ind)
			Sigs(:,:,lc)=Sig
			call linds(Sig,Siginvs(:,:,lc))
			call linds(Sig(1:q,1:q),SigXinvs(:,:,lc))
			Sigdets(lc)=detpd(Sig)
			SigXdets(lc)=detpd(Sig(1:q,1:q))
			thlist(:,lc)=[b,c,d]
			lc=lc+1
			
		end subroutine			
	endsubroutine

	subroutine mkSig(d0,c0,Sig) 
		use bsks_int
		implicit none
		
		real	:: d0,d,c0,c,Sig(:,:)
					
		integer				:: Tmax
		real, allocatable	:: v(:,:),vqd(:,:),ga(:)
		real		:: bseval(1),x
		integer		:: i,j
		
		Tmax=nint((1+maxval(rlist))*mkSig_T)+2	
		allocate(v(Tmax,qmax+size(rlist)),vqd(Tmax,qmax+size(rlist)),ga(-Tmax+1:Tmax-1))
		
		d=d0; if(d<0.5) d=d+1
		if(d==0.5) d=0.51
		c=c0; if(c0==0) c=0.01
		v=0
		do i=1,Tmax
			x=(i-0.5)/mkSig_T
			if(x<1) then
				do j=1,qmax
					v(i,j)=sqrt(2.0)*cos(pi*x*j)
				enddo
				v(i,qmax+1:)=-1
			else
				do j=1,size(rlist)
					v(i,qmax+j)=boole(x<1+rlist(j))/rlist(j)
				enddo
			endif
		enddo
		
		if(d0<0.5) then
			do i=1,Tmax-1
				vqd(i,:)=(c/mkSig_T)*v(i,:)-(v(i+1,:)-v(i,:))
			enddo
			vqd(1,:)=vqd(1,:)-v(1,:)
			vqd(Tmax,:)=0
			v=vqd*mkSig_T
		endif

		ga(0)=(c**(1 - 2*d)*Sqrt(Pi)*Gamma(-0.5 + d))/Gamma(d)
		do i=1,Tmax-1
			x=real(i)/mkSig_T
			if(abs(c*x)<700) then
				call bsks(-0.5 + d,c*x,1,bseval)
			else
				bseval=0
			endif
			ga(i)=(2**(1.5 - d)*Sqrt(Pi)*(x/c)**(-0.5 + d)*bseval(1))/Gamma(d)
			ga(-i)=ga(i)
		enddo
		do i=1,Tmax
			vqd(i,:)=matmul(ga(-i+1:Tmax-i),v)
		enddo
		Sig=matmul(transpose(v),vqd)/mkSig_T**2
	end subroutine
	
	
	function qform(Sigi,Y) result(val)
		implicit none
		real	:: Sigi(:,:), Y(:),val
		integer	:: i,j,k
		
		k=size(Y)
		val=0.0
		do i=1,k
			do j=i+1,k
				val=val+Y(i)*Sigi(i,j)*Y(j)
			enddo
		enddo
		val=2.0*val
		do i=1,k
			val=val+Y(i)**2*Sigi(i,i)
		enddo
	end function
	
	subroutine mkstats()
		implicit none
		
		real	:: SigsChol(q+1,q+1,n)
		real	:: densWfac,densXfac,densgfac
		real	:: cqf
		integer	:: i,l,j1,j2
		
		do i=1,n
			SigsChol(:,:,i)=choleski(Sigs(:,:,i))
		enddo
		
		densWfac=.5*pi**(-0.5*(q+1))*gamma(.5*(q+1))
		densXfac=.5*pi**(-0.5*q)*gamma(.5*q)
		densgfac=pi**(-.5*q)*gamma(.5*(q+1))/sqrt(2.0)
		do l=1,nsim
			call rnnoa(Ys(:,l))
		enddo
!$omp parallel do private(i,cqf) num_threads(6)			
		do l=1,nsim
			Ys(:,l)=matmul(SigsChol(:,:,(l-1)*n/nsim+1),Ys(:,l))
			Ys(:,l)=Ys(:,l)/norm2(Ys(1:q,l))
			
			do i=1,nws
				densWs(i,l)=denswfac*qform(Siginvs(:,:,i),Ys(:,l))**(-0.5*(q+1))/sqrt(Sigdets(i))
				cqf=qform(SigXinvs(:,:,i),Ys(1:q,l))
				densgs(i,l)=densgfac*cqf**(-0.5*(q+1))/sqrt(SigXdets(i))
				densXs(i,l)=densxfac*cqf**(-0.5*q)/sqrt(SigXdets(i))
			enddo
			do i=nws+1,n
				densWs(i,l)=denswfac*qform(Siginvs(:,:,i),Ys(:,l))**(-0.5*(q+1))/sqrt(Sigdets(i))
			enddo
			densavWs(l)=sum(densWs(:,l))/n
			ISw(:,l)=densWs(:,l)/densavWs(l)
		enddo	
	end subroutine
			
	function getavlength(indw,rej) result(val)
		implicit none
		integer	:: indw
		logical	:: rej(nsim)
		real	:: val
		val=sum(densgs(indw,:)/densavWs,mask=.not.rej)/nsim
	end function
	
	subroutine setlengths0() 
		implicit none
		real	:: w_prior0(nws)
		logical	:: b0rej(nsim)
		integer	:: j
		do j=1,nws
			w_prior0=0
			w_prior0(j)=1.0
			b0rej=getbayesrej(w_prior0)
			lengths0(j)=getavlength(j,b0rej)
		enddo
	end subroutine
		

	function getbayesrej(w_prior) result(val)
		implicit none
		
		real	:: w_prior(nws)
		logical	:: val(nsim)
		integer	:: l,j
		real	:: densyfac,num,den,C0,B0,A0,post

		densyfac=.5*pi**(-0.5*q)*gamma(.5*q)
		
!$omp parallel do private(num,den,A0,B0,C0,post)
		do l=1,nsim
			num=0;den=0
			do j=1,nws
				if(w_prior(j)==0) cycle
				B0=dot_product(Siginvs(1:q,q+1,j),Ys(1:q,l))
				A0=Siginvs(q+1,q+1,j)
				C0=qform(Siginvs(1:q,1:q,j),Ys(1:q,l))-B0**2/A0
				post=C0**(-.5*q)*A0**(-.5)*tdf(sqrt(q*A0/C0)*(Ys(q+1,l)+B0/A0),real(q))/sqrt(Sigdets(j))
				num=num+w_prior(j)*post
				den=den+w_prior(j)*densXs(j,l)
			enddo
			num=num*densyfac
			post=num/den
			val(l)=(post<level/2).or.(post>1-level/2)
		enddo
	end function


	subroutine setlam(w_length,lam,MNrej)
		implicit none
		real	:: lam(n),w_length(nws)
		logical	:: MNrej(nsim)
		real	:: RP(n),wlength0,wlengthc,avdensgs(nsim), avdgsnotbayes(nsim), avls(nws), lamstep
		real	:: oldRP(n), lamfac(n), cv
		integer	:: l,j

		avdensgs=matmul(w_length/lengths0,densgs(1:nws,:))
		avdgsnotbayes=avdensgs*boole(bayesrej)
		lamfac=4
		oldRP=level
		do j=1,500
			RP=0.0
!$omp parallel do reduction(+:RP) num_threads(6)
			do l=1,nsim
				if(dot_product(densWs(:,l),lam)<avdgsnotbayes(l)) RP=RP+ISw(:,l)
			enddo
			RP=RP/nsim
			lam=lam*exp(lamfac*(RP-level))
			where((RP-level)*(oldRP-level)>0)
				lamfac=lamfac*1.03
			elsewhere
				lamfac=lamfac/2.0
			endwhere
			oldRP=RP
			if(j>400) lamfac=0.1
			call setbounds(lamfac,0.01,20.0)
		enddo
		
		if(abs(maxval(RP-level))>0.002) then
			print '(a,f8.4,f)','problem in lambda determination: maximum of RP after lambda iterations',maxval(RP)
		endif

		if(cleanLamflag) then
			where(lam<0.001) lam=0.0
		endif

		cv=sum(lam)
		lam=lam/cv
		
		lamstep=0.1
		do j=1,12
			RP=0.0
!$omp parallel do reduction(+:RP) num_threads(6)			
			do l=1,nsim
				if(dot_product(densWs(:,l),cv*lam)<avdgsnotbayes(l)) RP=RP+ISw(:,l)
			enddo
			RP=RP/nsim
			if(sum(lam*RP)>level) then
				cv=cv*exp(lamstep)
			else
				cv=cv*exp(-lamstep)
			endif
			lamstep=lamstep/2
		enddo
		if(abs(sum(lam*RP)-level)>0.0001) print '(a,f12.4)','problem in cv determination: average RP',sum(lam*RP)
		lam=lam*cv

		wlength0=(1+lengthfudge)*sum(avdensgs/densavWs,mask=matmul(lam,densWs)>avdgsnotbayes)/nsim
		lamstep=0.5
		do j=1,15
			wlengthc=sum(avdensgs/densavWs,mask=matmul(lam,densWs)>avdgsnotbayes)/nsim
			if(wlength0>wlengthc) then
				lam=lam*exp(lamstep)
			else
				lam=lam*exp(-lamstep)
			endif
			lamstep=lamstep/2
		enddo
		if(abs(wlengthc/wlength0-1)>0.0005) then
			call VSprint("problem in length-correction--target and achieved:")
			call outML([wlength0,wlengthc])
		endif
		MNrej=matmul(lam,densWs)<avdgsnotbayes		
	end subroutine
	
	function getRPx(MNrej) result(RP)
		implicit none
		integer		:: l,i
		real		:: RP(nx), densWfac
		logical		:: MNrej(nsim)
		
		densWfac=.5*pi**(-0.5*(q+1))*gamma(.5*(q+1))
		RP=0
!$omp parallel do private(l) 
		do i=1,nx
			do l=1,nsim
				if(MNrej(l)) RP(i)=RP(i)+qform(Siginvs(:,:,i),Ys(:,l))**(-0.5*(q+1))/densavWs(l)
			enddo
			RP(i)=RP(i)/sqrt(Sigdets(i))
		enddo
		RP=RP*densWfac/nsim
	end function
	
	subroutine chksize(MNrej)
		implicit none
		logical		:: MNrej(nsim)
		integer		:: l,i
		real		:: RP(nx), RP2(nx), densWfac, ind
		
		densWfac=.5*pi**(-0.5*(q+1))*gamma(.5*(q+1))

		RP=getRPx(MNrej)
		
		i=maxloc(RP,dim=1)
		RP2=0
!$omp parallel do reduction(+:RP2)
		do l=1,nsim
			if(MNrej(l)) RP2(i)=RP2(i)+(qform(Siginvs(:,:,i),Ys(:,l))**(-0.5*(q+1))/densavWs(l))**2
		enddo
		RP2(i)=RP2(i)/Sigdets(i)
		RP2=RP2*densWfac**2/nsim
		call VSprint("largest RP, MCstdev and corresponding theta")
		call outML([RP(i),sqrt((RP2(i)-RP(i)**2)/nsim),thlist(:,i)])
		if(RP(i)>level+sqrt((RP2(i)-RP(i)**2)/nsim)) then 
			print *,'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
			print *,'XXXXXXXXXXXXXXX  Size Distortion !!! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
			print *,'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
		endif
	end subroutine
	
	function getavls(MNrej,inds) result(val)
		implicit none
		integer		:: inds(:),l
		real		:: val(size(inds))
		logical		:: MNrej(nsim)

		val=0
		do l=1,nsim
			if(.not.MNrej(l)) val=val+densgs(inds,l)/densavWs(l)
		enddo
		val=(val/nsim)/lengths0(inds)
	end function

	subroutine	finlengthq()
		implicit none
		integer		:: qlist(4)=[12,6,24,48]
		real		:: lam(n), outtab(size(levellist),size(qlist),2),lengths0q12(nws)
		integer		:: iq, ia,i
		logical		:: MNrej(nsim)
		
		rlist=rlist0
		outtab=-1
		qmax=48
		do ia=1,size(levellist)
			level=levellist(ia)
			do iq=1,size(qlist)
				q=qlist(iq)
				call mkSigs(q,1)
				call rnset(15)
				call mkstats
				print '(a,i4,f6.3,a)','XXXXXXXXX q,level, =',q,level,"  XXXXXXXXXXXXXX"
				if(iq==1) then
					call setlengths0
					lengths0q12=lengths0
				endif
				lengths0=lengths0q12
				bayesrej=getbayesrej(w_prior)
				lam=1.0
				call setlam(w_prior,lam,MNrej)
				call chksize(MNrej)
				outtab(ia,iq,1)=sum(w_prior*getavls(MNrej,[(i,i=1,nws)]))
				call setlengths0
				outtab(ia,iq,2)=sum(w_prior*lengths0/lengths0q12)
				print *,'theta unknown'
				call outML(outtab(:,:,1))
				print *,'theta known'
				call outML(outtab(:,:,2))
			enddo
		enddo

	end subroutine
	
	
	subroutine finlengthenv()
		implicit none
		integer	:: ia,l,i,ibayes,iout
		real	:: lam(n)
		integer, parameter		:: nenv=size(dlist)+size(blist)+size(clist)
		real, allocatable	:: w_env(:)
		integer	:: indenv(nenv)
		logical	:: MNrej(nsim)
		real	:: outtab(5,nenv,size(levellist))


		nws=size(dlist)*(1+size(blist)-1+size(clist)-1)
		deallocate(lengths0)
		allocate(lengths0(nws),w_env(nws))

		rlist=rlist0
		qmax=12
		useniceblist=.true.
		call mkSigs(q,1)
		call VSprint("Sigs computed")

		w_prior=[(1.0/size(dlist),i=1,size(dlist)),(0.0,i=size(dlist)+1,nws)]

		indenv(1:size(dlist))=[(i,i=1,size(dlist))]
		do i=1,size(blist)
			indenv(size(dlist)+i)=findind(niceblist(i),0.0,1.0)
		enddo
		do i=1,size(clist)
			indenv(size(dlist)+size(blist)+i)=findind(0.0,clist(i),1.0)
		enddo
		
		call mkstats
		do ia=1,size(levellist)
			level=levellist(ia)
			print '(a,f6.3,a)','XXXXXXXXX level=',level,"  XXXXXXXXXXXXXX"
			call setlengths0
			iout=1
			do ibayes=1,2
				if(ibayes==1) then
					bayesrej=getbayesrej(w_prior)
				else
					bayesrej=.true.
				endif
				lam=1.0
				call setlam(w_prior,lam,MNrej)
				call chksize(MNrej)
				outtab(iout,:,ia)=getavls(MNrej,indenv)
				iout=iout+1
			
				if(ibayes==1) then
					print *,'envelope with regular prior '
				else
					print *,'envelope with no Bayes interior set'
				endif
				do i=1,nenv
					w_env=0
					w_env(indenv(i))=1.0
					lam=1.0
					call setlam(w_env,lam,MNrej)
					call chksize(MNrej)
					outtab(iout,i:i,ia)=getavls(MNrej,[indenv(i)])
					print *,i,nenv
				enddo
				
				iout=iout+1
				if(ibayes==2) exit
				print *,'envelope with envelope prior '
				do i=1,nenv
					w_env=0
					w_env(indenv(i))=1.0
					bayesrej=getbayesrej(w_env)
					lam=1.0
					call setlam(w_env,lam,MNrej)
					call chksize(MNrej)
					outtab(iout,i:i,ia)=getavls(MNrej,[indenv(i)])	
					print *,i,nenv
				enddo
				iout=iout+1
			enddo
			
			call outML(thlist(:,indenv(1:size(dlist))).cud.outtab(:,1:size(dlist),ia))
			call outML(thlist(:,indenv(size(dlist)+1:size(dlist)+size(blist))).cud.outtab(:,size(dlist)+1:size(dlist)+size(blist),ia))
			call outML(thlist(:,indenv(size(dlist)+size(blist)+1:)).cud.outtab(:,size(dlist)+size(blist)+1:,ia))
		enddo
		
		contains
		
		function findind(b,c,d) result(val)
			implicit none
			real		:: b,c,d
			integer		:: val,i
		
			val=n+1
			do i=1,n
				if (all(thlist(:,i)-[b,c,d]==0)) then
					val=i
					exit
				endif
			enddo
			if(val>n) then 
				print *,"couldn't find corresponding index"
				stop
			endif
		end function
	end subroutine
	
	function getperf(trej) result(val)
		implicit none
		logical	:: trej(nsim)
		real	:: val(3),RPd(size(dlist)),RPx(nx)
		integer	:: i
		
		do i=1,size(dlist)
			RPd(i)=sum(ISw(i,:),mask=trej)/nsim
		enddo
		RPx=getRPx(trej)
		val=[sum(w_prior*RPd),maxval(RPd),maxval(RPx)]
	end function

	subroutine findistortions()
		implicit none
		integer	:: ia
		real	:: lam(n),outtab(3,3,size(levellist))
		logical	:: MNrej(nsim)
		
		rlist=rlist0
		qmax=12
		call mkSigs(q,1)
		call VSprint("Sigs computed")
		call mkstats

		do ia=1,size(levellist)
			level=levellist(ia)
			call setlengths0
			bayesrej=getbayesrej(w_prior)
			outtab(1,:,ia)=getperf(bayesrej)
			lam=0
			lam(1:nws)=1.0
			call setlam(w_prior,lam,MNrej)
			outtab(2,:,ia)=getperf(MNrej)
			lam=1.0
			call setlam(w_prior,lam,MNrej)
			outtab(3,:,ia)=getperf(MNrej)
			call outML(outtab(:,:,ia))
		enddo
		
	end subroutine
	
	subroutine fincomputeall()
		implicit none
		integer		:: qlist(5)=[6,9,12,24,48]
		real		:: lam(n),outlam(size(rlistx),n+1,size(levellist)), outlength(size(rlistx),nws+1,size(levellist))
		integer		:: iq, ir, ia
		logical		:: MNrej(nsim)
		
		saveSigflag=.true.
		rlist=rlistx
		cleanLamflag=.true.
		qmax=48
		do iq=1,size(qlist)
			q=qlist(iq)
			do ir=1,size(rlist)
				call mkSigs(q,ir)
				call rnset(15)
				call mkstats
				do ia=1,size(levellist)
					level=levellist(ia)
					print '(a,i4,2f6.3,a)','XXXXXXXXX q, r, level, =',q,rlist(ir),level,"  XXXXXXXXXXXXXX"
					call setlengths0
					outlength(ir,2:,ia)=lengths0
					outlength(:,1,ia)=rlistx
					bayesrej=getbayesrej(w_prior)
					lam=1.0
					call setlam(w_prior,lam,MNrej)
					call chksize(MNrej)
					outlam(ir,2:,ia)=lam
					outlam(:,1,ia)=rlistx
					call printtime
				enddo
			enddo
			do ia=1,size(levellist)
				call storeML(transpose(outlam(:,:,ia)),"lamtab")
				call VSprint("save('c:/out/int"//itostring(q)//"_"//itostring(nint(100*levellist(ia)))//".asc','lamtab','-ascii')")
				call execinML("save('c:/out/int"//itostring(q)//"_"//itostring(nint(100*levellist(ia)))//".asc','lamtab','-ascii')")

				call storeML(transpose(outlength(:,:,ia)),"ltab")
				call VSprint("save('c:/out/length0vec"//itostring(q)//"_"//itostring(nint(100*levellist(ia)))//".asc','ltab','-ascii')")
				call execinML("save('c:/out/length0vec"//itostring(q)//"_"//itostring(nint(100*levellist(ia)))//".asc','ltab','-ascii')")
			enddo
		enddo
	end subroutine		
		
	
end module

program mfort
	use globals
	use compute
	implicit none
	
	integer	:: i
	
	call openML
	call rnopt(9)
	call rnset(15)
	call inittime
	
	rlist=rlist0
	allocate(densWs(n,nsim),densgs(n,nsim),densXs(n,nsim),ISw(n,nsim)) 

	w_prior=[(1.0/nws,i=1,nws)]
	allocate(lengths0(nws))

!  Select desired computation by commenting in one of the next four
	
!	call findistortions				! Table 1 computation
!	call finlengthenv				! Table 2 computation
!	call finlengthq					! Table 3 computation
	call fincomputeall				! compute MN sets for many horizons and values of q
end program
