module Sigmod
	use globals
	implicit none
	
	integer		:: disttype
	integer, parameter	:: nBa=300
	real		:: materngrid(nBa,2)
	
	contains
	
	subroutine initmaterngrid
		use bskes_int
		real, parameter	:: nu=0.25
		integer	:: j
		real	:: x,ks(1)
		do j=1,nBa
			x=-8+14*(real(j)/nBa)
			materngrid(j,1)=x
			x=exp(x)
			call bskes(nu,x,1,ks)
			materngrid(j,2)=exp(nu*log(x)-x)*ks(1)
		enddo
		materngrid(:,2)=materngrid(:,2)/materngrid(1,2)
	end subroutine
	
	function materncov(r) result(val)
		use qdval_int
		real	:: r, val
		if(r<exp(-8.0)) then
			val=1
		else
			val=qdval(log(r),materngrid(:,1),materngrid(:,2),check=.false.)
		endif
	end function
	
	function getOmfromW(W,Sig) result(val)
		real	:: W(:,:),Sig(:,:),val(size(W,2),size(W,2))
		val=matmul(transpose(W),matmul(Sig,W))
	end function
	
	function get_dist(s1,s2) result(val)	! distance on surface of sphere given latitude and longitude
		real	:: s1(2),s2(2),val
		real	:: t1(2),t2(2),d(2),x(3)
		integer	:: i
		if(disttype==0) then
			val=norm2(s1-s2)
			return
		endif
		if(disttype==1) then
			t1=s1*pi/180
			t2=s2*pi/180
			d=t2-t1
			val = 2*asin(sqrt(sin(d(1) / 2)**2 + cos(t1(1)) * cos(t2(1)) * sin(d(2)/ 2)**2))
			return
		endif
!x = r cos(long) cos(lat) y = r sin(long) cos(lat) z = r sin(lat)
		t1=s1*pi/180
		x(1)= cos(t1(2))*cos(t1(1))
		x(2)= sin(t1(2))*cos(t1(1)) 
		x(3)= sin(t1(1))
		t2=s2*pi/180
		x(1)= x(1)-cos(t2(2))*cos(t2(1))
		x(2)= x(2)-sin(t2(2))*cos(t2(1)) 
		x(3)= x(3)-sin(t2(1))
		val=.5*norm2(x)		
	end function
	
	function distmat(s) result(val)
		real	:: s(:,:)
		real, allocatable	:: val(:,:)
		integer	:: n,i,j
		real	:: cmax
		n=size(s,2)
		allocate(val(n,n))
		cmax=0
		do i=1,n
			do j=1,i
				val(j,i)=get_dist(s(:,j),s(:,i))
				cmax=max(cmax,val(j,i))
			enddo
		enddo
		do i=1,n
			do j=1,i
				val(j,i)=val(j,i)/cmax
				val(i,j)=val(j,i)
			enddo
		enddo
	end function
				
	
	function getcSig(c,s) result(val)
		real	:: c,s(:,:)
		real, allocatable	:: val(:,:)
		val=exp(-c*distmat(s))
	end function	

	function getcpSig(c,p,s) result(val)
		real	:: c,s(:,:)
		real, allocatable	:: val(:,:),d(:,:)
		integer	:: p,i,j
		d=distmat(s)
		select case(p)
			case(-1)
				allocate(val,mold=d)
				do i=1,size(d,1)
					do j=1,size(d,2)
						val(j,i)=materncov(c*d(j,i))
					enddo
				enddo
			case(0)
				val=exp(-c*d)
			case(1)
				val=(1+c*d)*exp(-c*d)
			case(2)
				val=(1+c*d+(c*d)**2/3.0)*exp(-c*d)
			case default
				val=exp(-c*d**2)
		end select
	end function	
	
	function getetaSig(m,Sig0) result(Sig1)
		type(mtype)	:: m
		real		:: Sig0(:,:)
		real, allocatable	:: Sig1(:,:)
		integer		:: i,j,mmi,mmj,ci,cj
		allocate(Sig1(m%m,m%m))
		mmi=0
		do i=1,m%n
			mmj=0
			do j=1,m%n
				Sig1(m%clustid(j):m%clustid(j)+m%nclust(j)-1,m%clustid(i):m%clustid(i)+m%nclust(i)-1)=(Sig0(j,i)/(m%dxsi(i)*m%dxsi(j)))*matmul(reshape(m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1),[m%nclust(j),1]),reshape(m%x(m%clustid(i):m%clustid(i)+m%nclust(i)-1),[1,m%nclust(i)]))
				mmj=mmj+m%nclust(j)
			enddo
			mmi=mmi+m%nclust(i)
		enddo
	end function
	
	function geteSig(m,Sig0) result(Sig1)
		type(mtype)	:: m
		real		:: Sig0(:,:)
		real, allocatable	:: Sig1(:,:)
		real, parameter	:: rho=0.9
		real		:: Sige(capT,capT)
		integer		:: i,j
		do i=1,capT
			do j=1,capT
				Sige(j,i)=rho**abs(i-j)
			enddo
		enddo
		allocate(Sig1(m%m,m%m))
		do i=1,m%n
			do j=1,m%n
				Sig1(m%clustid(j):m%clustid(j)+m%nclust(j)-1,m%clustid(i):m%clustid(i)+m%nclust(i)-1)=Sig0(j,i)*Sige(1:m%nclust(j),1:m%nclust(i))
			enddo
		enddo
	end function
	
	function geteSig_het(m,Sig0) result(Sig1)
		type(mtype)	:: m
		real		:: Sig0(:,:)
		real, allocatable	:: Sig1(:,:)
		real, parameter	:: rho=0.9
		real		:: Sige(capT,capT)
		integer		:: i,j
		do i=1,capT
			do j=1,capT
				Sige(j,i)=rho**abs(i-j)
			enddo
		enddo
		allocate(Sig1(m%m,m%m))
		do i=1,m%n
			do j=1,m%n
				Sig1(m%clustid(j):m%clustid(j)+m%nclust(j)-1,m%clustid(i):m%clustid(i)+m%nclust(i)-1)=Sig0(j,i)*Sige(1:m%nclust(j),1:m%nclust(i))*		&
				matmul(reshape(sqrt(1+m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)**2),[m%nclust(j),1]),reshape(sqrt(1+m%x(m%clustid(i):m%clustid(i)+m%nclust(i)-1)**2),[1,m%nclust(i)]))
			enddo
		enddo
	end function

	subroutine setWx(m)
		type(mtype)	:: m
		integer		:: iq,j,q
		real		:: S(size(m%xZ,2),size(m%xZ,2))

		q=size(m%W,dim=2)-1

		if(allocated(m%Wx)) then
			if(.not. all(size(m%Wx)==[m%m,q+1])) then
				deallocate(m%Wx)
				allocate(m%Wx(m%m,q+1))
			endif
		else
			allocate(m%Wx(m%m,q+1))
		endif
		if(allocated(m%Weta)) then
			if(.not. all(size(m%Weta)==[m%n,q+1])) then
				deallocate(m%Weta)
				allocate(m%Weta(m%n,q+1))
			endif
		else
			allocate(m%Weta(m%n,q+1))
		endif

		m%Wx(:,1)=m%x
		do j=1,m%n
			do iq=1,q
				m%Wx(m%clustid(j):m%clustid(j)+m%nclust(j)-1,1+iq)=sqrt(real(m%n))*m%W(j,iq+1)*m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)
			enddo
		enddo
		S=matmul(transpose(m%xZ),m%xZ)
		call erset(0,0,0)
		call linds(S,S)
		if(iercd()>0) then
			m%n=-1
			return
		endif
		m%Wx(:,2:)=m%Wx(:,2:)-matmul(m%xZ,matmul(S,matmul(transpose(m%xZ),m%Wx(:,2:))))
		do j=1,m%n
			do iq=1,q+1
				m%Weta(j,iq)=sum(m%Wx(m%clustid(j):m%clustid(j)+m%nclust(j)-1,iq)*m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1))/m%dxsi(j)
			enddo
		enddo
	end subroutine

	function getWeigen(Sig,q) result(val)
		use evesf_int
		real	:: Sig(:,:)
		integer	:: q
		real	:: c,val(size(Sig,1),q+1)
		real	:: Sigm(size(Sig,1),size(Sig,1)),evs(q)
		integer	:: i,n
		n=size(Sig,1)
		do i=1,n
			Sigm(:,i)=Sig(:,i)-sum(Sig(:,i))/n
		enddo
		do i=1,n
			Sigm(i,:)=Sigm(i,:)-sum(Sigm(i,:))/n
		enddo
		call evesf(q,Sigm,.false.,evs,val(:,2:q+1))
		val(:,1)=1.0/sqrt(real(n))
	end function

	
	function getavcorr(Sig) result(val)
		real	:: Sig(:,:),val
		integer	:: i,j,n
		n=size(sig,1)
		val=0
		do i=1,n
			do j=1,i-1
				val=val+Sig(j,i)
			enddo
		enddo
		val=val/(n*(n-1)/2)
	end function
	
	function getavc_c(avc,p,s) result(c)
		use rootfinder
		real	:: avc,s(:,:),c
		integer	:: p
		integer	:: n,i,j,l
		real	:: c0,c1,tol,f,maxd
		real, allocatable	:: dns(:)
		integer(4)	:: rstat
		real	:: sx
		
		n=size(s,2)
		allocate(dns(n*(n-1)/2))
		l=1;maxd=0
		do i=1,n
			do j=1,i-1
				dns(l)=get_dist(s(:,i),s(:,j))
				maxd=max(maxd,dns(l))
				l=l+1
			enddo
		enddo
		dns=dns/maxd
		
		c0=5; c1=50
		do i=1,5
			if(getavcd(c0)>0) exit
			c0=.5*c0
		enddo
		if(i>5) then
			print *,"cannot find c0"
			stop
		endif
		do i=1,20
			if(getavcd(c1)<0) exit
			c1=2*c1
		enddo
		if(i>20) then
			print *,"cannot find c1 for avc=",avc
			stop
		endif
		rstat=0; tol=0.001
		call zero_rc(c0,c1,tol,c,rstat,f)
		do i=1,20
			f=getavcd(c)
			call zero_rc(c0,c1,tol,c,rstat,f)
			if(rstat==0) exit
		enddo
		if(i>20) then
			print *,"cannot locate c"
			stop
		endif
		
	contains
		function getavcd(c) result(val)
			real	:: c,val
			select case(p)
				case(-1)
					val=0
					do j=1,size(dns)
						val=val+materncov(c*dns(j))
					enddo
					val=val/size(dns)
				case(0)
					val=sum(exp(-c*dns))/size(dns)
				case(1)
					val=sum((1+c*dns)*exp(-c*dns))/size(dns)
				case(2)
					val=sum((1+c*dns+(c*dns)**2/3.0)*exp(-c*dns))/size(dns)
				case default
					val=sum(exp(-c*dns**2))/size(dns)
				end select
			val=val-avc	
		end function
	end function	

	
	subroutine setWHR(m)
		use evesf_int
		type(mtype)	:: m
		real	:: Xsi(m%m,m%n), evs(m%n), evcs(m%m,m%n),MxZ(m%m,m%m),val(m%m,m%n+1)
		integer	:: i,j,n,q
		
		call erset(0,0,0)
		MxZ=eye(m%m)-matmul(m%xZ,matmul(invertpd(matmul(transpose(m%xZ),m%xZ)),transpose(m%xZ)))
		if(iercd()>0) print *,"issue in MxZ generation",iercd()		
		Xsi=0
		do j=1,m%n
			Xsi(m%clustid(j):m%clustid(j)+m%nclust(j)-1,j)=m%x(m%clustid(j):m%clustid(j)+m%nclust(j)-1)
		enddo
		if(allocated(m%WHR)) deallocate(m%WHR)
		allocate(m%WHR(m%m,m%n+1))
		m%WHR(:,2:)=matmul(MxZ,Xsi)
		m%WHR(:,1)=m%x
	end subroutine
	
	function getKVSig(s,b) result(Sig)
		real	:: s(:,:),b,Sig(size(s,2),size(s,2))
		if(b==0) then
			Sig=eye(size(s,2))
		else
			Sig=max(1.0-distmat(s)/b,0.0)
		endif
	end function
	
	function robchol(mat) result(val)
		use evcsf_int
		real	:: mat(:,:),val(size(mat,1),size(mat,1))
		real	:: evs(size(mat,1)),evecs(size(mat,1),size(mat,1))
		integer	:: j		
		call evcsf(mat,evs,evecs)
		do j=1,size(mat,1)
			val(:,j)=sqrt(max(evs(j),0.0))*evecs(:,j)
		enddo
	end function	
		
	function getWKV(m,b) result(val)
		type(mtype)	:: m
		real		:: b
		real	:: val(m%m,m%n+1)
		if(b==0) then
			val=m%WHR
			return
		endif
		val(:,2:)=matmul(m%WHR(:,2:),robchol(getKVSig(m%s,b)))
		val(:,1)=m%x
	end function

end module
			
