module globals
	use myfuncs
	use dotops
	use ML
    integer             :: nsim, k, q
    integer, parameter 	:: nsim0=2500, nsimf=10000
    real, parameter     :: level=0.05
    
    real, allocatable   :: Ga(:,:,:),sumGainv(:,:),sGair(:,:),Vs(:,:),eps(:,:)
    
    logical             :: DinDflag, BCHflag
 
end module



module compute
	use globals
    implicit none

    contains
    
    subroutine setGa(ik)
		use svrgp_int

        implicit none
        integer     :: ik,i,isort(q)
        real        :: v(2*k*k), A(2*k,k),vsort(q),Gax(k,k,q)
        
        if(allocated(Ga)) deallocate(Ga,sumGainv,sGair,Vs)
        allocate(Ga(k,k,q), sumGainv(k,k), sGair(k-1,k-1),Vs(k,q))    

        sumGainv=0.0
        Ga=0
        do i=1,q
            call rnnoa(v)
            if(DinDflag) then
                if(i<=ik+1) then 
                    Ga(2,2,i)=v(1)**2+v(2)**2
                else 
                    Ga(:,:,i)=v(1)**2+v(2)**2
                endif
            else
                A=reshape(v,[2*k,k])
                Ga(:,:,i)=matmul(transpose(A),A)
            endif
            vsort(i)=Ga(1,1,i)
            sumGainv=sumGainv+Ga(:,:,i)
        enddo

        if(k>1) call linds(sumGainv(2:,2:),sGair)
		call linds(sumGainv,sumGainv)
        isort=[(i,i=1,q)]
        call svrgp(vsort,vsort,isort)
        Gax=Ga
        do i=1,q
            Ga(:,:,i)=Gax(:,:,isort(i))
        enddo
        
    end subroutine
   
  
    function gettstat(Y) result(out)
        implicit none
        real    :: Y(k,q), out
        real    :: beta(k),s2
        integer :: i
        
        beta=matmul(sumGainv,sum(Y,2))
        s2=0.0
        do i=1,q
            s2=s2+sum(sumGainv(1,:)*(Y(:,i)-matmul(Ga(:,:,i),beta)))**2
        enddo
        out=beta(1)/sqrt(q*s2/(q-1))
    end function      
        
    function getRP() result(RP)
        implicit none
        real        :: RP
        integer     :: lc,i,j, jx, rejcount
        real        :: Y(k,q),betar(k-1),e(k,q),e0(k,q),t0,qs(2), tstats(2**q)
        real        :: cv
       rejcount=0
       cv=tin(1-level/2,real(q-1))
!$omp parallel do private (i,Y,t0,betar,e,e0,j,jx,tstats,qs) reduction(+:rejcount)
       do lc=1,nsim
            do i=1,q
                Y(:,i)=eps(i,lc)*Vs(:,i)
            enddo
            t0=gettstat(Y)
            if(BCHflag) then
                if (abs(t0)>cv) rejcount=rejcount+1
            else
                betar=matmul(sGair,sum(Y(2:,:),2))
                do i=1,q
                    e(:,i)=Y(:,i)-matmul(Ga(:,2:,i),betar)
                enddo
                e0=e
                do j=1,2**q
                   jx=j
                    do i=1,q
                        if(mod(jx,2)==1) then 
                            e(:,i)=-e0(:,i) 
                        else  
                            e(:,i)=e0(:,i) 
                        endif
                        jx=jx/2
                    enddo
                    tstats(j)=gettstat(e)
                enddo
                qs=quantile(tstats,[level/2,1-level/2])
                if ((t0<qs(1)) .or. (t0>qs(2))) rejcount=rejcount+1
            endif
        enddo
       RP=real(rejcount)/nsim
                
    end function
    
    subroutine setVs(czero,cV)
        implicit none
        integer :: czero,cV,i
        Vs=0

        if(czero>q) then
            Vs(:,q:czero-q+1:-1)=1
        else
            Vs(:,1:czero)=1
        endif
        
        if(DinDflag) then
            do i=1,q
                if(Ga(1,1,i)==0.0) Vs(1,i)=0.0
            enddo
            return
        endif

        do i=2,k
            select case(mod(cV/3**(i-2),3))
				case (1)
					Vs(i,:)=0
				case (2)
					Vs(i,:)=-Vs(i,:)
            end select
        enddo
    endsubroutine

    subroutine seteps()
        implicit none

        integer :: i
        
        if (allocated(eps)) deallocate(eps)
        allocate(eps(q,nsim))
        do i=1,nsim
            call rnnoa(eps(:,i))
        enddo
    endsubroutine
        
end module

program mfort
	use globals
	use compute
	implicit none
    integer     :: lc, ik, id, iM, iq, czero, cV, czeromax, cVmax, ncV, nk
    integer, parameter  :: nlc=100, qlist(3)=[4,8,12]
    
    real        :: maxRP, cRP, outtab(8,10,3)
	real        :: maxRPs(nlc)
    real, allocatable        :: epsf(:,:), eps0(:,:)
    
	call openML()
    call rnopt(5)
    call rnset(14)
	call inittime()

    outtab=-1
    do iq=1,3
		q=qlist(iq)

		allocate(eps0(q,nsim0),epsf(q,nsimf))
		nsim=nsim0
		call seteps()
		eps0=eps
		nsim=nsimf
		call seteps()
		epsf=eps
    

		do id=1,2
			if (id==1) then 
				DinDflag=.true.
				nk=q/2-1
			else 
				DinDflag=.false.
				nk=3
			endif
			do ik=1,nk
				if(DinDflag) then 
					k=2
					nCV=1
				else 
					k=ik
					nCV=3**(k-1)
				endif
    
				do iM=1,2
					if(iM==1) then 
						BCHflag=.false.
					else 
						BCHflag=.true. 
					endif
					do lc=1,nlc
						call setGa(ik)
						maxRP=-1.0
						nsim=nsim0
						eps=eps0
						do czero=1,2*q-1
							do cV=0,ncV-1
								call setVs(czero,cV)
								cRP=getRP()
								if(cRP>maxRP) then
									maxRP=cRP
									czeromax=czero
									cVmax=cV
								endif
							enddo
						enddo    
						nsim=nsimf
						eps=epsf
						call setVs(czeromax,cVmax)
						cRP=getRP()
						maxRPs(lc)=cRP
					enddo
					outtab(ik+5*(id-1),5*(iM-1)+1:,iq)=[minval(maxRPs),quantile(maxRPs,[.25,.5,.75]),maxval(maxRPs)]
					do lc=1,3
						call outML(outtab(:,:,lc))
					enddo
					call printtime()
				enddo
			enddo
		enddo
		deallocate(eps0,epsf)
    enddo
    
end program

