module compute_reg
!DIR$ NOOPTIMIZE
	use compute_gen
	use draw
	implicit none

	
	contains

	function tstat result(stat)
		logical, save	:: init=.true.
		integer, save	:: rind(1000)
		real, save		:: dir(200,10)
		real, allocatable	:: stat(:),v(:)
		integer			:: i,j,c
		
		if(init) then
			rind=1234872*ranu_v(1000)
			do i=1,size(dir,2)
				call rnnoa(dir(:,i))
			enddo
			init=.false.
		endif
		stat=[log(om2),log(sigma2_F),real(cindF),muU,mt_f]
		c=1
		do i=1,24
			select case(i)
			case(1)
				v=[Xs]
			case(2)
				v=[f]
			case(3)
				v=[log(kappa2)]
			case(4)
				v=[CFs]
			case(5)
				v=[rhos]
			case(6)
				v=[log(clubs2)]
			case(7)
				v=[CoCFs]
			case(8)
				v=[CoCrhos]
			case(9)
				v=[log(CoCs2)]
			case(10)
				v=[log(clubs2)]
			case(11)
				v=[clubid]
			case(12)
				v=[CoCid]
			case(13)
				v=[cindu]
			case(14)
				v=[cindCF]
			case(15)
				v=[cindCoC]
			case(16)
				v=[rhodist]
			case(17)
				v=[CoCrhodist]
			case(18)
				v=[thdist]
			case(19)
				v=[CFthdist]
			case(20)
				v=[CoCthdist]
			case(21)
				v=[s2dist]
			case(22)
				v=[CFs2dist]
			case(23)
				v=[CoCs2dist]
			case(24)
				v=[fm]

			end select
			do j=1,2
				stat=[stat,v(mod(rind(c),size(v))+1)]
!				stat=[stat,v(j)]
				c=c+1
			enddo
		enddo
		stat=[stat,matmul(transpose(dir(1:size(stat),:)),stat)]
	end function

	
	subroutine draw_once
		call draw_X
		if(trialrun) call draw_Y
		call draw_CF
		call draw_CoCF
		call draw_rho
		call draw_CoCrho
		call draw_rhodist
		call draw_CoCrhodist
		call draw_kappa2grid
		call draw_clubs2grid
		call draw_CoCs2grid
		call draw_s2dist
		call draw_CFs2dist
		call draw_CoCs2dist
		call draw_CoCid
		call draw_clubid
		call draw_cindu
		call draw_cindCF
		call draw_cindCoC
		call draw_thdist
		call draw_CFthdist
		call draw_CoCthdist
		call draw_muU
		call draw_om2
  !
		call draw_mt_F
		call draw_F
		call draw_fm
		call draw_sigma2_F
		call draw_cindf
	end subroutine	

	subroutine iterate
		integer	:: l,leff
		real	:: Delta0(0:q,0:q)
		real, allocatable	:: tracestats(:,:)
		
		Delta0=Delta
		do l=1,nburnin
			if(mod(l,20)==1) Delta=Delta0*(1000.0)**(max(0,nburnin/2-l)/(0.5*nburnin))
			call linds(Delta,Deltainv)
			call draw_once
		enddo
		call mdisp(matmul(G,f))
		allocate(tracestats(size(tstat()),ns))
		
		do l=1,nskip*ns
			call draw_once
!			call mdisp(Us(:,5).cvr.Xs(:,5))
!			call mdisp([real(leff),sigma2_F,om2])
			if(mod(l,nskip)==0) then 
				call draw_fcsts
				call draw_fm
				leff=l/nskip
				tracestats(:,leff)=tstat()
				if(mod(leff,100)==0) then
					call mdisp([real(leff),sigma2_F,om2])
					call mdisp(sum(exp(tracestats(1:2,1:leff)),dim=2)/leff)
					call storeML(tracestats(:,1:leff),"ts")
				endif
!				call savestats(leff)
			endif
        enddo
        call printtime
		
	end subroutine
end module