module runmod
	use compute
	use scpcmod
	use Markov
	implicit none
	
	real, parameter		:: avc0=0.03
	integer, parameter	:: ntab=12, ncase=14, ntabmis=6, ntablength=8
	integer, parameter	:: nbgrid=20
	
	contains
	
	function getquants(tab) result(val)
		real	:: tab(:,:), val(6,size(tab,2))
		real, allocatable	:: ctab(:)
		integer	:: j
		do j=1,size(tab,2)
			ctab=pack(tab(:,j),tab(:,j)>=0)
			val(:,j)=quantile_v(ctab,[.05,.25,.5,.75,.95])
			val(6,j)=sum(ctab)/size(ctab)
		enddo
	end function
	
	function getMClength(m) result(val)
		type(mtype) :: m
		real, parameter	:: avc0s(ntablength/2)=[0.1,0.03,0.01,0.003]
		real		:: c0,val(ntablength),len
		real, allocatable	:: Sig(:,:),W(:,:),Om(:,:)
		real		:: cv,ccv
		integer		:: q,i
		do i=1,size(avc0s)
			c0=getavc_c(avc0s(i),0,m%s)
			Sig=getcSig(c0,m%s)
			W=getWeigen(Sig,50+1)
			call setSCPCqcv(W,m%s,c0,q,cv)
			m%W=W
			call setWx(m)
			ccv=max(cv,getetacv(m,c0))
			Om=matmul(transpose(m%Wx),m%Wx)
			len=getlength(Om)
			val(2*i-1)=len*cv/(2*1.96*sqrt(real(m%m)))
			val(2*i)=len*ccv/(2*1.96*sqrt(real(m%m)))
		enddo
	end function
			
	function getMCstat(m,c0) result(val)
		type(mtype) :: m
		real		:: c0,val(ntab)
		real, allocatable	:: Om(:,:),Sig0(:,:),Sig(:,:),W(:,:),WKV(:,:,:)
		real		:: cv,ccv,brpgrid(nbgrid)
		integer		:: q,i
		val=-1
		Sig=getcSig(c0,m%s)
		W=getWeigen(Sig,qmax+1)
		call setSCPCqcv(W,m%s,c0,q,cv)
		m%W=W
		call setWx(m)
		ccv=max(cv,getetacv(m,c0))
		
		Sig0=getcSig(c0,m%s)
		
		call setWHR(m)		
		allocate(WKV(m%m,m%n+1,nbgrid))
!$omp parallel do		
		do i=1,nbgrid
			WKV(:,:,i)=getWKV(m,((i-1.0)/(nbgrid-1))**2)
		enddo
		
		Sig=geteSig(m,Sig0)

		Om=matmul(transpose(m%Wx),matmul(Sig,m%Wx))
		val(1:2)=[getrp(Om,cv),getrp(Om,ccv)]
		
!$omp parallel do 
		do i=1,nbgrid	
			brpgrid(i)=getrp(matmul(transpose(WKV(:,:,i)),matmul(Sig,WKV(:,:,i))),1.96)
		enddo
		val(3)=minval(brpgrid)
		val(4)=brpgrid(1)

		m%Wx(:,2:)=m%Wx(:,2:)*ccv
		val(5)=robust_avcx(m,avc0,.false.)
! now eta		
		Om=matmul(transpose(m%Weta),matmul(Sig0,m%Weta))
		val(7:8)=[getrp(Om,cv),getrp(Om,ccv)]
		
		Sig=getetaSig(m,Sig0)
!$omp parallel do 
		do i=1,nbgrid
			brpgrid(i)=getrp(matmul(transpose(WKV(:,:,i)),matmul(Sig,WKV(:,:,i))),1.96)
		enddo
		val(9)=minval(brpgrid)
		val(10)=brpgrid(1)	
		m%Weta(:,2:)=m%Weta(:,2:)*ccv
		val(11)=robust_avcx(m,avc0,.true.)
	end function
	
	function getMCstat_mis(m,c0) result(val)
		type(mtype) :: m
		real		:: c0,val(ntabmis)
		real, allocatable	:: Om(:,:),Sig0(:,:),Sig(:,:),W(:,:)
		real		:: cv,ccv
		integer		:: q,i
		val=-1
		Sig=getcSig(c0,m%s)
		W=getWeigen(Sig,qmax+1)
		call setSCPCqcv(W,m%s,c0,q,cv)
		m%W=W
		call setWx(m)
		ccv=max(cv,getetacv(m,c0))
		Sig0=getcSig(c0,m%s)

		Sig=geteSig(m,Sig0)

		Om=matmul(transpose(m%Wx),matmul(Sig,m%Wx))
		val(1:2)=[getrp(Om,cv),getrp(Om,ccv)]
		
		Sig=geteSig_het(m,Sig0)

		Om=matmul(transpose(m%Wx),matmul(Sig,m%Wx))
		val(3:4)=[getrp(Om,cv),getrp(Om,ccv)]
		
		Sig0=getcSig(getavc_c(0.1,0,m%s),m%s)
		Sig=geteSig(m,Sig0)

		Om=matmul(transpose(m%Wx),matmul(Sig,m%Wx))
		val(5:6)=[getrp(Om,cv),getrp(Om,ccv)]
	end function
	
	subroutine setm_case(icase,panelflag,m,c0)
		use rlse_int
		integer, parameter	:: capT=4
		integer	:: icase
		type(mtype) :: m
		logical	:: panelflag
		real	:: c0
		real	:: ux(m%n*capT),uZ(m%n*capT,3),b(m%n+5)
		integer	:: j,k,i

		call rnnoa(ux)
		do j=1,size(uZ,2)
			call rnnoa(uZ(:,j))
		enddo
		if(allocated(m%nclust)) deallocate(m%nclust,m%clustid)
		allocate(m%nclust(m%n),m%clustid(m%n))
		
		if(panelflag) then
			m%nclust=capT
			m%m=m%n*capT
			m%clustid=[(capT*(i-1)+1,i=1,m%n)]
			
		else
			m%m=m%n
			m%nclust=1
			m%clustid=[(i,i=1,m%n)]
		endif
		c0=getavc_c(avc0,0,m%s)
		select case(icase)
			case(1)
				m%xZ=reshape(ones(m%m),[m%m,1])
			case(2)
				m%xZ=drawx(c0).cvc.ones(m%m)
			case(3)
				m%xZ=reshape(drawx(c0),[m%m,1])
			case(4)
				m%xZ=drawx(c0).cvc.ones(m%m).clr.drawZ(c0)
			case(5)
				m%xZ=ux(1:m%m).cvc.ones(m%m)
			case(6)
				m%xZ=ux(1:m%m).cvc.ones(m%m).clr.drawZ(c0)
			case(7)
				m%xZ=step(0.5).cvc.ones(m%m)
			case(8)
				m%xZ=step(0.15).cvc.ones(m%m)
			case(9)
				m%xZ=step(0.5).cvc.ones(m%m).clr.drawZ(c0)
			case(10)
				m%xZ=step(0.15).cvc.ones(m%m).clr.drawZ(c0)
			case(11)
				m%xZ=((step(0.5)+1.0)*drawx(c0)).cvc.ones(m%m)
			case(12)
				m%xZ=((step(0.15)+1.0)*drawx(c0)).cvc.ones(m%m)
			case(13)
				m%xZ=drawx(.1).cvc.ones(m%m)
			case(14)
				m%xZ=drawx(.1).cvc.ones(m%m).clr.drawZ(c0)
			case default
				print *,"error, no such case"
				stop
		end select
		if(allocated(m%x)) deallocate(m%x,m%Z,m%dxsi)	
		allocate(m%x(m%m),m%Z(m%m,merge(size(m%xZ,2)-1,size(m%xZ,2)-1+m%n-1,.not. panelflag)),m%dxsi(m%n))	
		m%Z=0
		if(panelflag) then
			do j=1,m%n
				m%z(m%clustid(j):m%clustid(j)+m%nclust(j)-1,1:size(m%xZ,2)-1)=m%xZ(m%clustid(j):m%clustid(j)+m%nclust(j)-1,2:)	
				if(j<m%n) m%z(m%clustid(j):m%clustid(j)+m%nclust(j)-1,size(m%xZ,2)+j-1)=1.0
			enddo
			m%xZ=m%xZ(:,1).clr.m%Z
		endif
		k=size(m%xZ,2)
		if(k>1) then
			call rlse(m%xZ(:,1),m%xZ(:,2:),b(1:k-1),intcep=0)
			m%xZ(:,1)=m%xZ(:,1)-matmul(m%xZ(:,2:),b(1:k-1))
			if(sum(m%xZ(:,1)**2)<1E-8) then
				print *,"collinearity in case", icase, panelflag
				call mdisp(m%xZ)
				stop
			endif
		endif
		m%xscale=sqrt(sum(m%xZ(:,1)**2)/m%m)
		m%xZ(:,1)=m%xZ(:,1)/m%xscale
		do j=1,m%n
			m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)=m%xZ(m%clustid(j):m%clustid(j)+m%nclust(j)-1,1)
			if(k>1 .and. (.not. panelflag)) m%Z(m%clustid(j):m%clustid(j)+m%nclust(j)-1,:)=m%xZ(m%clustid(j):m%clustid(j)+m%nclust(j)-1,2:)
			m%dxsi(j)=sqrt(sum(m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)**2))
			if(m%dxsi(j)==0) m%dxsi(j)=1
		enddo
		
	contains
	
		function drawx(c0) result(val)
			real :: c0
			real :: val(m%m)
			val=matmul(choleski(geteSig(m,getcSig(c0,m%s))),ux(1:m%m))
		end function
		
		function drawZ(c0) result(val)
			real :: c0
			real :: val(m%m,size(uZ,2))
			val=matmul(choleski(geteSig(m,getcSig(c0,m%s))),uZ(1:m%m,:))
		end function
		
		function step(fr) result(val)
			real	:: fr,val(m%m)
			integer	:: inds(m%n),j
			real	:: ss(merge(capT,1,panelflag))
			real	:: s(m%n)
			inds=[(j,j=1,m%n)]
			call svrgp(m%s(1,:),s,inds)
			if(panelflag) then
				ss=[(merge(1,0,2*j<=capT),j=1,capT)]
			else
				ss=1
			endif
			do j=1,m%n
				val(m%clustid(inds(j)):m%clustid(inds(j))+m%nclust(inds(j))-1)=ss*merge(1.0,0.0,j<fr*m%n)
			enddo
		end function
	end subroutine
	
	subroutine MCstates
		integer, parameter	:: ndraws=5,nseries=nstates*ndraws,n=250
		type(mtype) :: m
		real	:: tab(ntab,nseries,ncase),tabmis(ntabmis,nseries,ncase),tablength(ntablength,nseries,ncase)
		integer	:: i,istate,ip,icase
		logical	:: panelflag
		real	:: c0

		allocate(m%s(2,n))
		do ip=1,2
			panelflag=ip==2
			tab=-1; tabmis=-1
			do i=1,nseries
				istate=(i-1)/ndraws+1
				m%n=n
				m%s=getstatelocs(istate,1.0,m%n)
				do icase=merge(2,1,panelflag),ncase
					call setm_case(icase,panelflag,m,c0)
					tab(:,i,icase)=getMCstat(m,c0)
					tabmis(:,i,icase)=getMCstat_mis(m,c0)
					tablength(:,i,icase)=getMClength(m)
				enddo
			enddo		
			do icase=1,ncase
				print *,"light MC: panel=",panelflag, ", case=",icase
				print *,"default panel of tests:"
				call mdisp(getquants(transpose(tab(:,:,icase))))
				print *,"heteroskedasticity and misspecification of c0:"
				call mdisp(getquants(transpose(tabmis(:,:,icase))))
				print *,"excess length in i.i.d. model"
				call mdisp(getquants(transpose(tablength(:,:,icase))))
				print *,"XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
			enddo

		enddo
	end subroutine

	subroutine MCWB
		type(mtype) :: m
		real	:: c0,cv,ccv
		real, allocatable	:: Sig(:,:),W(:,:),Om(:,:)
		real	:: tab(ntab,nseries),tabmis(ntabmis,nseries)
		integer	:: q,i,k,l,ip,ife
		logical	:: panelflag=.false.

		call setWBpanel(m,4,1,.true.)
		call setWB(m,4)
		do ip=1,2
			panelflag=ip==2
			do k=1,4,3
				do ife=1,merge(2,1,panelflag)
					tab=-1; tabmis=-1
					do i=1,nseries
						if(panelflag) then
							call setWBpanel(m,k,i,ife==2)
						else
							call setWB(m,k,i)
						endif
						if(m%n<=0) cycle
						c0=getavc_c(avc0,0,m%s)
!						tab(:,i)=getMCstat(m,c0)
						tabmis(:,i)=getMCstat_mis(m,c0)
					enddo
					print *,"MCWB: k=",k,"panel=",panelflag,"time fe=",ife==2
!					print *,"default panel of tests:"
!					call mdisp(getquants(transpose(tab)))
					print *,"heteroskedasticity and misspecification of c0:"
					call mdisp(getquants(transpose(tabmis)))
					print *,"XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
				enddo
			enddo		
		enddo
		call printtime
	end subroutine	
	
	
	function KVrej(m,ls) result(val)
		use rlse_int
		use iercd_int
		type(mtype)	:: m
		integer		:: ls(:)
		real		:: val(nbgrid)
		real		:: ehat(m%m),beta(size(m%xZ,2)),KVSig(m%n,m%n),u(m%n),num
		integer		:: j,i
		
		call rlse(m%y,m%xZ,beta,intcep=0)
		if(iercd()>0) print *,"issue in HRref y xZ reg",iercd()
		ehat=m%y-matmul(m%xZ,beta)
		do j=1,m%n
			u(j)=sum(m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)*ehat(m%clustid(j):m%clustid(j)+m%nclust(j)-1))
		enddo
		num=sum(m%y*m%x)
		do i=1,nbgrid
			val(i)=num**2>1.96**2*sum(u*matmul(getKVSig(m%s,((i-1.0)/(nbgrid-1))**2),u))
		enddo
	end function
	
	subroutine MarkovWB
		integer, parameter	:: ns=50,nsim=5000,ndraws=200
		real, parameter		:: c0max=50
		type(mtype)	:: m0,m
		integer		:: ls(ns)
		real, allocatable	:: Pi(:,:),u(:),us(:),W(:,:),Sig(:,:)
		real		:: cr(2+nbgrid),c0,cv,t2
		real :: tab(4,ndraws)
		integer		:: q,l,lc,lb,imax,cc,j,k,ik,i,ip
		logical	:: panelflag=.false.

		do ip=1,2
		panelflag=ip==2
		do k=1,4,3
			lc=0
			do
				if(.not. panelflag) then
					call setWB(m0,k)
				else
					call setWBpanel(m0,k)
				endif
			
				if(allocated(u)) deallocate(u)
				allocate(u(m0%n))
				do j=1,m0%n
					u(j)=sum(m0%y(m0%clustid(j):m0%clustid(j)+m0%nclust(j)-1)*m0%x(m0%clustid(j):m0%clustid(j)+m0%nclust(j)-1))
				enddo
				if(sum(u*matmul(getKVSig(m0%s,0.3),u))<3*sum(u**2)) cycle
				lc=lc+1
				Pi=getMarkov_Pi(avc0,u,distmat(m0%s),ns)
				if(Pi(1,1)==-1) cycle
				call set_sumPi(Pi)
				cr=0
				cc=0
!$omp parallel do private(m,ls,c0,q,W,cv,Sig,t2,us,j) reduction(+:cr,cc) 		
				do l=1,nsim
					ls=getMarkovSample(Pi,ns)
					call setmfromls(m,m0,ls)
					if(m%n<0) then
						cc=cc+1
						cycle
					endif
					cr(3:)=cr(3:)+KVrej(m,ls)
					if(getavcorr(getcSig(c0max,m%s))>avc0) then
						c0=c0max
					else
						c0=getavc_c(avc0,0,m%s)
					endif
					m%W=getWeigen(getcSig(c0,m%s),min(qmax,m%n-k))
					call setSCPCqcv(m%W,m%s,c0,q,cv)

					call setWx(m)

					if(m%n<=0 .or. all(sum(m%Wx(:,2:)**2,dim=1)<1E-8) .or. all(sum(m%Weta(:,2:)**2,dim=1)<1E-8)) then
						cr=cr-1E100
						cycle
					endif					
					t2=sum(m%Wx(:,1)*m%y)**2/sum(matmul(m%y,m%Wx(:,2:))**2)
					cr(1)=cr(1)+boole(t2>cv**2)
					cv=max(cv,getetacv(m,c0))
					cr(2)=cr(2)+boole(t2>cv**2)
				enddo
				cr=cr/nsim
				tab(:,lc)=[cr(1),cr(2),minval(cr(3:)),cr(3)]
				if(lc==ndraws) exit
			enddo	
			print *,"Markov: panel=",panelflag,", k=",k
			call mdisp(getquants(transpose(tab)))
			call printtime
		enddo
		enddo
	end subroutine	
	
	function KVstats(m,WKV) result(val)
		use rlse_int
		use iercd_int
		type(mtype)	:: m
		real		:: WKV(:,:,:)
		real		:: val(nbgrid)
		real		:: ehat(m%m),beta(size(m%xZ,2)),KVSig(m%n,m%n),u(m%n),v(m%n),num
		integer		:: j,i
		
		call rlse(m%y,m%xZ,beta,intcep=0)
		if(iercd()>0) print *,"issue in HRref y xZ reg",iercd()
		ehat=m%y-matmul(m%xZ,beta)
		do j=1,m%n
			u(j)=sum(m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)*ehat(m%clustid(j):m%clustid(j)+m%nclust(j)-1))
		enddo
		do i=1,nbgrid
			v=matmul(u,WKV(:,:,i))
			val(i)=sum(v**2)
		enddo
	end function
	
end module
	
	
program mfort
	use runmod
	implicit none
	integer	:: i

	call initmaterngrid
	call rnset(14)
!$omp parallel do
	do i=1,100
		call rnset(49184*i)
	enddo
	call inittime
	
	call mkGQxw(GQxw)
	print *,"now entering MCstates"
	call MCstates

	call printtime
	print *,"now entering MCWB"
	call MCWB

	call printtime
	print *,"now entering MarkovWB"
	call MarkovWB
	call printtime
end program