module output
!	!DIR$ NOOPTIMIZE
	use compute_reg
	use dgp
	use prep

	implicit none
	
	contains
	
	function getpwcorrs() result(val)
		real	:: val(3,(n-1)*n/2),corr(3,n,n)
		real	:: v(2),var(n)
		integer	:: ir1,ir2,ir,j
		integer	:: i1,i2,k1,k2
		integer	:: ic
!$omp parallel do private(ir2,i1,i2,k1,k2,v)
		do ir1=1,n
			do ir2=1,n
				i1=clubid(ir1); i2=clubid(ir2)
				k1=CoCid(i1); k2=CoCid(i2)
				v=0
				if(i1==i2) v=v+om2*rhos(ir1)*rhos(ir2)*(1-CoCrhos(i1)**2)*clubs2(i1)*[1.0,Suss(:,cindCF(i1))]
				if(k1==k2) v=v+om2*CoCrhos(i1)*CoCrhos(i2)*rhos(ir1)*rhos(ir2)*CoCs2(k1)*[1.0,Suss(:,cindCoC(k1))]
				if(ir1==ir2) then
					v=v+om2*(1-rhos(ir1)**2)*kappa2(ir1)*[1.0,Suss(:,cindu(ir1))]
				endif
				corr(1,ir1,ir2)=v(2)+sigma2_F(1)*sfmss(1,cindf)+sigma2_F(2)*sfass(1)
				corr(2,ir1,ir2)=v(2)
				corr(3,ir1,ir2)=v(1)
			enddo
		enddo
		do ic=1,3
			var=diagonal(corr(ic,:,:))
			j=1
			do ir1=1,n
				do ir2=1,ir1-1
					val(ic,j)=corr(ic,ir1,ir2)/sqrt(var(ir1)*var(ir2))
					j=j+1
				enddo
			enddo
		enddo
	end function
	
	function getpaths result(val)
		real	:: val(npath,0:n)
		integer	:: ir
		val(1:T,0)=matmul(G(:,0:q0),f(0:q0))
		val(T+1:,0)=fcstf(1:nh)+val(T,0)
!$omp parallel do
		do ir=1,n
			val(1:T,ir)=matmul(G(:,0:q0),Us(0:q0,ir))
			val(T+1:,ir)=fcstUs(1:nh,ir)+val(T,ir)
		enddo
	end function
	
	subroutine savestats(leff,postflag)
		integer	:: leff
		logical	:: postflag
		
		integer	:: l,ir,j,i,k,ir1,ir2
		real	:: gvec(size(gammas,1))

		real, save, allocatable	:: rfvals(:,:),rcindudist(:,:)
		real, allocatable, save	:: rcorr(:,:,:),rpaths(:,:,:),rfm(:,:)
		
		real, save	:: rhl(n,ns),rfdist(nf), rvola(2,n,ns), rfstats(6,ns)
		real, save	:: rcstats(2,n,ns)
		real, save	:: rrvolaU0(2,ns)
		real, save	:: rrhodist(2,nrhogrid)
		real, save	:: rmuUs(ns)
		real, save	:: rss(ntdel,3,n,ns),ccorr(n,n,ntdel)
		real, save	:: rdelfss(ns)
		
		real, save	:: rcindu(n,ns),rcindCF(nclubs,ns),rcindCoC(nCoCs,ns),rkappa2(n,ns),rclubs2(nclubs,ns),rCoCs2(nCoCs,ns),rrhos(n,ns),rCoCrhos(nclubs,ns),rclubid(n,ns),rCoCid(nclubs,ns)
		real, save	:: rkappas(n,ns)
		
		if(leff==1) then
			if(allocated(rcorr)) deallocate(rcorr,rpaths,rfvals,rcindudist,rfm)
			allocate(rcorr(3,(n-1)*n/2,ns),rpaths(npath,0:n,ns),rfvals(0:q,ns),rcindudist(nth,n),rfm(0:q0,ns))
			rfdist=0; rrhodist=0	
			rcindudist=0
		endif			

		rpaths(:,:,leff)=getpaths()
		rmuUs(leff)=muU
		do ir=1,n
			rcindudist(cindu(ir),ir)=rcindudist(cindu(ir),ir)+1
		enddo
		rkappas(:,leff)=kappa2
		
		rcindu(:,leff)=cindu
		rcindCF(:,leff)=cindCF
		rcindCoC(:,leff)=cindCoC
		rkappa2(:,leff)=kappa2
		rclubs2(:,leff)=clubs2
		rCoCs2(:,leff)=CoCs2
		rrhos(:,leff)=rhos
		rCoCrhos(:,leff)=CoCrhos
		rclubid(:,leff)=clubid
		rCoCid(:,leff)=CoCid

		rfdist(cindf)=rfdist(cindf)+1
		rfstats(:,leff)=[frhotab(cindf),sigma2_F,mt_f,om2]
		rfvals(:,leff)=f(0:q0)
		rfm(:,leff)=fm(0:q0)
				
		rdelfss(leff)=sqrt(sigma2_F(1)*sfmss(1,cindf))
		
!$omp parallel do private(gvec,i,k)
		do ir=1,n
			gvec=(1-rhos(ir)**2)*kappa2(ir)*gammas(:,cindu(ir))
			i=clubid(ir)
			gvec=gvec+rhos(ir)**2*(1-CoCrhos(i)**2)*clubs2(i)*gammas(:,cindcf(i))
			k=CoCid(i)
			gvec=gvec+rhos(ir)**2*CoCrhos(i)**2*CoCs2(k)*gammas(:,cindCoC(k))
			gvec=gvec/gvec(1)
			rhl(ir,leff)=count(gvec(2:)>.5)
		enddo
		block
			real	:: vx(2,n), vy(2,n)
			do ir=1,n
				vx(:,ir)=[1.0,Suss(:,cindu(ir))]*(1-rhos(ir)**2)*kappa2(ir)
				i=clubid(ir)
				vx(:,ir)=vx(:,ir)+[1.0,Suss(:,cindcf(i))]*rhos(ir)**2*(1-CoCrhos(i)**2)*clubs2(i)
				k=CoCid(i)
				vx(:,ir)=vx(:,ir)+[1.0,Suss(:,cindCoC(k))]*rhos(ir)**2*CoCrhos(i)**2*CoCs2(k)
				vx(:,ir)=sqrt(vx(:,ir))
				vy(1,ir)=sum(G(1,0:q0)*Us(0:q0,ir))
				vy(2,ir)=sum(G(T,0:q0)*Us(0:q0,ir))-vy(1,ir)
			enddo
			rcstats(:,:,leff)=vy(:,:)
					
			rvola(:,:,leff)=vx
			do i=1,2
				rrvolaU0(1,leff)=getcorr(vx(2,:),vy(1,:))
				rrvolaU0(2,leff)=getcorr(vx(2,:),vy(2,:))
			enddo
		end block
				
		rcorr(:,:,leff)=getpwcorrs()

		rrhodist(1,:)=rrhodist(1,:)+rhodist
		rrhodist(2,:)=rrhodist(2,:)+CoCrhodist

		if(leff==ns) then
		block			
			real	:: growth(5,n)
			integer, parameter	:: h=50
			do ir=1,n
				growth(:,ir)=quantile_v(rpaths(T+h,ir,1:leff)-rpaths(T,ir,1:leff),[0.05,0.16,0.5,0.84,0.95])/h
			enddo
			call mydispnames(names)
			print *,"stdev del50(U)"
			call mdisp(sum(rvola(2,:,1:leff),dim=2)/leff)
			print *,"U growth quantiles"
			call mdisp(growth)
			do ir=1,n
				growth(:,ir)=quantile_v(rpaths(T+h,ir,1:leff)-rpaths(T,ir,1:leff)+rpaths(T+h,0,1:leff)-rpaths(T,0,1:leff),[0.05,0.16,0.5,0.84,0.95])/h
			enddo
			print *,"Y growth quantiles"
			call mdisp(growth)
			print *,"f growth quantiles"
			call mdisp(quantile_v(rpaths(T+h,0,1:leff)-rpaths(T,0,1:leff),[0.05,0.16,0.5,0.84,0.95])/h)
			call mdisp(sum(rkappas(:,1:leff),dim=2)/leff)
		end block
!		stop
		print *,"starting pushing to Matlab"
		call execinML("clear all")
		call storeML(fweights,"fweights")			
		call storeML(rfvals(0:q0,:),"F")
		call storeML(rfstats(1,:),"frho") 
		call storeML(rfstats(2:3,:),"sigma2_F")
		call storeML(rfstats(4:5,:),"mtau_F")
		call storeML(rfstats(6,:),"om2")
		call storeML(rfm,"m_of_f")
		call storeML(rdelfss,"stddev_del50_m")
		call storeML(G(:,0:q0),"g_matrix")
		call storeML(rhl,"halflife")
		call storeML(rvola(1,:,:),"stddev_U")
		call storeML(rvola(2,:,:),"stddev_del50U")
		call storeML(rcstats(1,:,:),"U0")
		call storeML(rcstats(2,:,:),"UT_U0")
		call storeML(rcorr(1,:,:),"pwcorr_del50Y")
		call storeML(rcorr(2,:,:),"pwcorr_del50U")
		call storeML(rcorr(3,:,:),"pwcorr_U")
		call storeML(rcindu,"cindu")
		call storeML(rcindCF,"cindCF")
		call storeML(rcindCoC,"cindCoC")
		call storeML(rkappa2,"kappa2")
		call storeML(rclubs2,"clubs2")
		call storeML(rCoCs2,"CoCs2")
		call storeML(rrhos,"rhos")
		call storeML(rCoCrhos,"CoCrhos")
		call storeML(rclubid,"clubid")
		call storeML(rCoCid,"CoCid")
		call storeML(rmuUs,"muU")
		call storeML(hlpriordist,"hl_prior")
		call storeML(tscorrprop,"r1r2w1")
		call storeML(rpaths(:,0,:),"path_f")			
		call storeML(rpaths(:,1:n,:),"paths_U")	
		postmeanf=sum(rfvals(:,:),dim=2)/ns
!	print *,"done"
!	stop
			!block
			!	integer		:: ind(T+2)
			!	ind(1:T)=[(i,i=1,T)]
			!	ind(T+1:)=T+[50,100]
   !
			!	call storeML(rpaths(ind,0,:),"path_f")
			!	call storeML(rpaths(ind,1:n,:),"paths_U")
			!end block
			if(postflag) then
				call execinML("save('"//trim(dir)//"postdraws_"//trim(cname)//".mat')")
			else
!				call execinML("save('"//trim(dir)//"priordraws.mat')")
			endif
			call mydispnames(names)
		endif
	end subroutine
	
	subroutine savepostdraws
		integer	:: l,leff
		real	:: Delta0(0:q,0:q)
		real, allocatable	:: tracestats(:,:)
		
		Delta0=Delta
		do l=1,nburnin
			if(mod(l,20)==1) then
				Delta=Delta0*(1000.0)**(max(0,nburnin/2-l)/(0.5*nburnin))
				call linds(Delta,Deltainv)
			endif
			if(2*l<nburnin) sigma2_F(2)=max(sigma2_F(2),1E-5)
			call draw_once
		enddo
		call mdisp(matmul(G,f))
		allocate(tracestats(size(tstat()),ns))

		do l=1,nskip*ns
			call draw_once
			if(mod(l,nskip)==0) then 
				call draw_fcsts
				leff=l/nskip
				call savestats(leff,.true.)
!				tracestats(:,leff)=tstat()
				if(mod(leff,1000)==0) print *,leff
			endif
		enddo
!		call savemat(trim(dir)//"tracestats.csv",tracestats)
        call printtime
	end subroutine

	subroutine savepriordraws
		integer	:: l
		
		do l=1,ns
!			call dgp_sigma2_F
			call dgp_cindf
!			call dgp_mt_F
			call dgp_F
			call dgp_s2dist
			call dgp_thdist
			call dgp_CFthdist
			call dgp_CoCthdist
			call dgp_clubid
			call dgp_CoCid
			call dgp_rhodist
			call dgp_CoCrhodist
			call dgp_cindu
			call dgp_cindCF
			call dgp_cindCoC
			call dgp_rho
			call dgp_CoCrho
			call dgp_CoCs2grid
			call dgp_clubs2grid
			call dgp_kappa2grid
!			call dgp_om2
			call dgp_CoCF
			call dgp_CF
!			call dgp_muU
			call dgp_U

			call draw_fcsts
			call savestats(l,.false.)
			if(mod(l,1000)==0) print *,l
		enddo
        call printtime
	end subroutine
end module