module Markov
	use globals
	implicit none
	
	contains
	
	subroutine balanceSig(Sig)
		real	:: Sig(:,:)
		real	:: x(size(Sig,1))
		integer	:: i,j,n
		x=1
		do i=1,20
			x=1.0/matmul(Sig,x)
		enddo
		n=size(Sig,1)
		do i=1,n
			do j=1,n
				Sig(j,i)=x(j)*x(i)*Sig(j,i)
			enddo
		enddo
		Sig=n*Sig/sum(sig)
		do i=1,n
			Sig(i,:)=Sig(i,:)/sum(Sig(i,:))
		enddo
	end subroutine
	
	function avcfromPi(u,Pi,n) result(val)
		real	:: u(:),Pi(:,:),val
		integer	:: n
		real	:: A(size(u),size(u))
		integer	:: i,j
		A=Pi
		val=0
		do i=1,n-1
			val=val+(n-i)*sum(u*matmul(A,u))
			A=matmul(Pi,A)
		enddo
		val=val*2.0/((n-1)*n*sum(u**2))
	end function
	
	function getMarkov_Pi(avc,u,distmat,n) result(Pi)
		use rootfinder
		real	:: avc,distmat(:,:),u(:),Pi(size(u),size(u))
		integer	:: n,i
		real	:: c,c0,c1,tol,f
		integer(4)	:: rstat
		
		
		c0=2; c1=20; i=0
		do
			i=i+1
			if(getavcd(c0)*getavcd(c1)<0 .or. i>40) exit
			c0=c0/2; c1=c1*2
		enddo
		if(i>40) then
			print *,"issue in getMarkov_Pi 1"
			print *,c0, getavcd(c0), c1, getavcd(c1)
			Pi=-1
			return
		endif
		rstat=0; tol=0.001
		call zero_rc(c0,c1,tol,c,rstat,f)
		i=0
		do
			i=i+1
			f=getavcd(c)
			call zero_rc(c0,c1,tol,c,rstat,f)
			if(rstat==0.or.i>40) exit
		enddo
		if(i>40) then
			print *,"issue in getMarkov_Pi 2"
			Pi=-1
			return
		endif
		
	contains
		function getavcd(c) result(val)
			real	:: c,val
			integer	:: j
			Pi=exp(-c*distmat)
			do j=1,size(Pi,1)
				Pi(j,j)=0
			enddo
			call balanceSig(Pi)
			val=avcfromPi(u,Pi,n)		
			val=val-avc	
		end function
	end function
	
	subroutine set_sumPi(Pi)
		real	:: Pi(:,:)
		integer	:: i,j
		do i=1,size(Pi,1)
			do j=2,size(Pi,2)
				Pi(i,j)=Pi(i,j-1)+Pi(i,j)
			enddo
			Pi(i,:)=Pi(i,:)/Pi(i,size(Pi,2))
		enddo
	end subroutine
	
	function getMarkovSample(SPi,n) result(val)
		real	:: SPi(:,:)
		integer	:: n,val(n)
		real	:: u(n)
		integer	:: capN,i
		call rnun(u)
		capN=size(Spi,1);
		val(1)=capN*u(1)+1
		do i=2,n
			val(i)=max(min(count(u(i)>SPi(val(i-1),:))+1,capN),1)
		enddo
	end function

end module