Matrix models for quantifying competitive intransitivity from species abundance data
Werner Ulrich, Santiago Soliveres, Wojciech Kryszewski, Fernando T. Maestre, Nicholas J. Gotelli
The executable of this code is available at: http://www.keib.umk.pl/transitivity/

!Last change 15.12.2013
Module definitions
real*8 ran1
integer time1,time2,time3,time0,iseed,filen,length,nu,file
integer sp,si,spmin,spmax,simin,simax,varsp,varsi,iter,fill,run,replic,testpower,number,npairs,npairs1,npairs2,nsim
real up,down,r2,f,diag,diag1,diag2,int1,eps,epsor,spear,spearnew,degree,transmetric,transmetric1,transmetric2,prop,prop1,prop2,chlen,chlen1,chlen2
real transmean,transdown,transup,transprop,diagmean,diagdown,diagup,diagprop,propmean,propdown,propup,propprop,chlenmean,chlendown,chlenup,chlenprop 
real, allocatable:: mat(:,:),mat1(:,:),mat2(:,:),mat3(:,:),varmat(:,:),ptemp(:,:),ctemp(:,:),interaction(:,:),transition(:,:),rowa(:),cola(:),temp(:,:),trans(:),trans1(:),temp2(:,:),temp3(:,:),err(:,:)
real, allocatable:: row(:),col(:),transsim(:),propsim(:),diagsim(:),chlensim(:)
real*8, allocatable:: out(:,:),temp1(:,:),eigenvector(:,:),eigenvalue(:),yvar(:),bvar(:),evar(:)
character*10 time
character*10 rand,special,rep,dir,batch,mtype
character*65000 text
character*250 name1,name2,t,method,method1
character*50 tt(5000)
character*50, allocatable:: matspname(:),matsiname(:),varmatspname(:),varmatsiname(:)
end module definitions

program Transitivity
use definitions
use eigensolving
use mat_handling
use matrix_basics
use basic_stats
use mult_reg
use sorting
real*8 norm
call date_and_time(time=time)
read(time(1:2),'(i2)')time1
read(time(3:4),'(i2)')time2
read(time(5:6),'(i2)')time3
call copy
call random_seed
call random_number(ran1)
iseed=int(ran1*(10**9))
time0=time1*3600+time2*60+time3
special='no'
ijk=0
ijk1=0
nsim=100
number=100000
iter=200
method=' '
rep='n'
batch='n'
epsor=0.975
mtype=' '
if(special.eq.'yes')then
open(7,file='Special.txt',status='replace')
endif
202 ijk=ijk+1
ijk1=ijk1+1
length=65000
iter1=100000
filen=1
text=' '
irep=0
name1=' '
name2=' '
run=1
fill=0
testpower=0
if(rep.ne.'y')then
write(*,*)'Which method to use: environmental data (e),'
write(*,*)'time series (t), or abundance data only (n)(default)'
read(*,'(a)',advance='yes',eor=8)method
8 if(method.eq.' ')method='trial'
if(method.eq.'n')method='trial'
if(method.eq.'e')method='environment'
if(method.eq.'t')method='time'
if(method.eq.'trial'.or.method.eq.'time'.or.method.eq.'environment')then
write(*,*)'Number of random matrices to generate (default=100,000)?'
read(*,'(i6)',advance='yes',eor=92)number
92 if(number.le.0)number=100000
write(*,*)'Type of competition matrix: bivariate (b) or probability(p, default)?'
read(*,'(a)',advance='yes',eor=9)mtype
9 if(mtype.eq.' ')mtype='p'
endif
if(mtype.eq.'p')degree=0.00
if(mtype.eq.'b')degree=0.00
write(*,*)'Name of output file (with extension), default = Output.txt'
read(*,'(a)',advance='yes',eor=1)t
1 if(t.eq.' ')t='Output.txt'
open(3,file=t,status='replace')
t=' '
!write(3,'(a4,17x,a6,1x,23(a9,1x))')'File','Method','BenchM','Species','Sites','Metric1','Mean','DownCL','UpCL','Prop<1','Metric2','Mean','DownCL','UpCL','Prop<1','PropInTrS','Mean','DownCL','UpCL','Prop<1','TransChL','Mean','DownCL','UpCL','Prop<1'
write(3,'(a4,17x,a6,1x,8(a9,1x))')'File','Method','BenchM','Species','Sites','Metric1','Mean','DownCL','UpCL','Prop<1'
!if(method.eq.'trial')write(3,'(a4,17x,a6,1x,9(a9,1x))')'File','Method','BenchM','Species','Sites','RankCorr','Metric1','Mean','DownCL','UpCL','Prop<1'
write(*,*)'Name of matrix file (with extension), default = Matrix.txt'
read(*,'(a)',advance='yes',eor=2)t
2 if(t.eq.' ')t='Matrix.txt'
open(4,file=t,status='replace')
write(*,*)'Batch run: yes (y), or no (n, default)?'
read(*,'(a)',advance='yes',eor=3)batch
3 if(batch.eq.'y')then
write(*,*)'Name of batch file'
read(*,'(a)')text
open(2,file=text,status='old',err=3)
read(2,*)text
rep='y'
endif
endif   
if(special.eq.'yes')then
endif
if(batch.eq.'y')then
call reading (tt,nu,2)
name1=tt(1)
name2=tt(2)
if(name1.eq.' ')goto 1001
write(*,'(2(a20,1x))')name1,name2
endif
5 if(batch.ne.'y')then
write(*,*)'Name of species x sites file with extension.'
read(*,'(a)',advance='yes',eor=6)name1
6 if(name1.eq.' ')goto 5
if(method.ne.'environment')goto 7
write(*,*)'Name of variables x sites file with extension.'
read(*,'(a)',advance='yes',eor=7)name2
7 if(name2.eq.' ')name2='***'
endif
if(name1.ne.' ')then
open(filen,file=name1,status='old',err=202)
call matreading (mat,matspname,matsiname,sp,si,filen)
close(filen)
if(name2.ne.'***'.and.name2.ne.' ')then
!method='environment'
open(filen,file=name2,status='old',err=202)
call matreading (varmat,varmatspname,varmatsiname,varsp,varsi,filen)
close(filen)
endif
else
goto 5
endif
190 allocate(rowa(sp),cola(si),col(si),row(sp),trans(sp),trans1(sp),interaction(sp,sp),transition(sp,sp),transsim(nsim),diagsim(nsim),propsim(nsim),chlensim(nsim))
allocate(eigenvector(sp,1),eigenvalue(sp))
write(*,'(a,a12)')'Method: ',method
method1=' '
200 interaction=0
transition=0
rowa=0
cola=0
row=0
col=0
fill=0
abund=0
do 15 j=1,si
a=1000000000
do 14 i=1,sp
if(mat(i,j).gt.0)then
fill=fill+1
if(mat(i,j).lt.a)a=mat(i,j)
abund=abund+mat(i,j)
rowa(i)=rowa(i)+mat(i,j)
cola(j)=cola(j)+mat(i,j)
!row(i)=rowa(i)+1
!col(j)=col(j)+1
endif
14 continue
do 16 i=1,sp
16 mat(i,j)=mat(i,j)/a
15 continue
write(4,'(a)')'Species abundance matrix'
call filewrite(mat,matspname,matsiname,sp,si,4)
write(4,*)' '
write(4,'(a)')'Row total distribution'
call filewrite(rowa,matspname,matsiname,sp,1,4)
write(4,*)' '
if(method.eq.'trial'.or.method1.eq.'trial')then
dir='down'
call colstandard(mat,sp,si)
call colstandard(rowa,sp,1)
call sortrow(mat,rowa,matspname,sp,si,1,1,dir)
!number=max(100000,sp*si)
if(sp.le.2)then
write(*,*)'Too few species'
goto 1001 
endif
!a=0
!spear=0
!do 20 j=1,si-1
!do 21 j1=j+1,si
!do 23 i=1,sp
!trans(i)=mat(i,j)
!trans1(i)=mat(i,j1)
!23 continue
!spear=spear+spearman(trans,trans1,sp)
!a=a+1
!21 continue
!20 continue
!if(a.gt.0)spear=spear/a
if(mtype.eq.'p')call mtop(mat,transition,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
if(mtype.eq.'b')call mtoc(mat,interaction,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
if(mtype.eq.'p')call metric(transition,sp,n,npairs,transmetric,diag,prop,chlen)
if(mtype.eq.'b')call metricpairs(interaction,sp,n,npairs,transmetric,diag,prop,chlen)
if(mtype.eq.'b')call ctop(interaction,transition,sp,degree)
call powereigen(transition,sp,1,it,eigenvector,eigenvalue)
transmean=mean(transsim,nsim)
!diagmean=mean(diagsim,nsim)
!propmean=mean(propsim,nsim)
!chlenmean=mean(chlensim,nsim)
call conflimnew(transsim,nsim,transup,transdown,1.0,transprop)
!call conflimnew(diagsim,nsim,diagup,diagdown,1.0,diagprop)
!call conflimnew(propsim,nsim,propup,propdown,1.0,propprop)
!call conflimnew(chlensim,nsim,chlenup,chlendown,1.0,chlenprop)
!write(3,'(a20,1x,a6,1x,f9.3,1x,2(i9,1x),20(f9.3,1x))')name1,mtype,eps,sp,si,transmetric,transmean,transdown,transup,transprop,diag,diagmean,diagdown,diagup,diagprop,prop,propmean,propdown,propup,propprop,chlen,chlenmean,chlendown,chlenup,chlenprop
write(3,'(a20,1x,a6,1x,f9.3,1x,2(i9,1x),5(f9.3,1x))')name1,mtype,eps,sp,si,transmetric,transmean,transdown,transup,transprop
write(4,'(26a)')'Best fit transition matrix'
call filewrite(transition,matspname,matspname,sp,sp,4)
write(4,'(a4,17x,a7,14x,a19,1x,a20)')'File','Species','Obs. rel. abundance','Dominant eigenvector'
do 191 i=1,sp
191 write(4,'(a20,1x,a20,11x,f9.6,12x,f9.6)')name1,matspname(i),rowa(i),real(eigenvector(i,1))
write(4,*)' '
if(mtype.eq.'b')then
write(4,'(27a)')'Best fit competition matrix'
call filewrite(interaction,matspname,matspname,sp,sp,4)
endif
elseif(method.eq.'environment')then
if(varsi.ne.si)then
write(*,'(a,1x,i3,1x,i3)')'Site mismatch',varsi,si
goto 1001 
endif
allocate(out(si,varsp),yvar(si),bvar(varsp+1),evar(si))
call rowstandard(varmat,varsp,si) 
out=transpose(varmat)
do 100 i=1,sp
do 101 j=1,si
101 yvar(j)=alog(mat(i,j)+0.1)
call multreg(out,yvar,si,varsp,bvar,r2,f,evar)
do 102 j=1,si
102 mat(i,j)=exp(evar(j))
100 continue
method1='trial'
!write(4,'(a20)')'Environmental matrix'
!call filewrite(real(out),matspname,matsiname,si,varsp,4)
!write(4,*)' '
!write(4,'(a20)')'Adjusted abundance matrix'
!call filewrite(mat,matspname,matsiname,sp,si,4)
write(4,*)' '
deallocate(out,evar,yvar,bvar,varmat)
goto 200
elseif(method.eq.'time')then
dir='down'
call colstandard(mat,sp,si)
call colstandard(rowa,sp,1)
call sortrow(mat,rowa,matspname,sp,si,1,1,dir)
allocate(out(sp,sp),mat1(sp,si-1),mat2(sp,si-1),mat3(si-1,sp),temp(sp,sp))
if(si.le.sp)then
write(*,'(a,1x,i3,1x,i3)')'Site mismatch',sp,si
goto 1001 
endif
do 90 i=1,sp
do 91 j=1,si-1
mat1(i,j)=mat(i,j)
mat2(i,j)=mat(i,j+1)
91 continue
90 continue
!mat3=transpose (mat1)
!call dotproduct (mat1,mat3,out,sp,sp,si-1)
!call FINDInv(real(out),temp,sp,testpower)
!if(testpower.eq.-1)then
!write(3,'(a20,1x,f9.3,1x,4(i9,1x),7(f9.3,1x))')name1,eps,sp,si,0,n,0,0,0,0,0,0,0
!goto 120
!endif
!call dotproduct (mat2,mat3,out,sp,sp,si-1)
!call dotproduct (real(out),temp,out,sp,sp,sp)
!transition=real(out)
!call colstandard(transition,sp,sp) 
if(mtype.eq.'p')then
call ttop(mat1,mat2,transition,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
call powereigen(transition,sp,1,it,eigenvector,eigenvalue)
call metric(transition,sp,n,npairs,transmetric,diag,prop,chlen)
elseif(mtype.eq.'b')then
call ttop(mat1,mat2,transition,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
call powereigen(transition,sp,1,it,eigenvector,eigenvalue)
call ttoc(mat1,mat2,interaction,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
!call ptoc(transition,interaction,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
call metricpairs(interaction,sp,n,npairs,transmetric,diag,prop,chlen)
endif
transmean=mean(transsim,nsim)
!diagmean=mean(diagsim,nsim)
!propmean=mean(propsim,nsim)
!chlenmean=mean(chlensim,nsim)
call conflimnew(transsim,nsim,transup,transdown,1.0,transprop)
!call conflimnew(diagsim,nsim,diagup,diagdown,1.0,diagprop)
!call conflimnew(propsim,nsim,propup,propdown,1.0,propprop)
!call conflimnew(chlensim,nsim,chlenup,chlendown,1.0,chlenprop)
!write(3,'(a20,1x,a6,1x,f9.3,1x,2(i9,1x),20(f9.3,1x))')name1,mtype,eps,sp,si,transmetric,transmean,transdown,transup,transprop,diag,diagmean,diagdown,diagup,diagprop,prop,propmean,propdown,propup,propprop,chlen,chlenmean,chlendown,chlenup,chlenprop
write(3,'(a20,1x,a6,1x,f9.3,1x,2(i9,1x),5(f9.3,1x))')name1,mtype,eps,sp,si,transmetric,transmean,transdown,transup,transprop
write(4,*)'Best fit transition matrix'
call filewrite(transition,matspname,matspname,sp,sp,4)
write(4,'(a4,17x,a7,14x,a19,1x,a20)')'File','Species','Obs. rel. abundance','Dominant eigenvector'
do 193 i=1,sp
193 write(4,'(a20,1x,a20,11x,f9.6,12x,f9.6)')name1,matspname(i),rowa(i),real(eigenvector(i,1))
write(*,*)' '
if(mtype.eq.'b')then
write(4,*)'Best fit competition matrix'
call filewrite(interaction,matspname,matspname,sp,sp,4)
endif
120 deallocate(out,mat1,mat2,mat3,temp)
endif
deallocate(mat,matspname,matsiname,transition,interaction,transsim,diagsim,propsim,chlensim,row,rowa,col,cola,trans,trans1,eigenvalue,eigenvector)
name=' '
name1=' '
name2=' '
write(*,*)' '
if(special.eq.'yes')then 
endif
if(batch.eq.'y')goto 202
write(*,*)'Another run? (y/n), default = n'
read(*,'(a)',advance='yes',eor=201)rep
201 if(rep.eq.'y')then
goto 202
endif
goto 1001
1000 write(*,*)'Format error in data file'
1001 call date_and_time(time=time)
read(time(1:2),'(i2)')time1
read(time(3:4),'(i2)')time2
read(time(5:6),'(i2)')time3
time1=time1*3600+time2*60+time3-time0
time2=int(time1/3600)
time3=int(time1/60)-time2*60
time1=time1-time3*60-time2*3600
write(*,'(a,i3,a,i3,a,i3,a)')'Runtime of program: ',&
time2,'hh',time3,'min',time1,'sec'
write(3,'(a)')' '
write(3,'(a,i3,a,i3,a,i3,a)')'Runtime of program: ',&
time2,'hh',time3,'min',time1,'sec'
pause
stop
end

subroutine copy
use definitions
write(*,*)'*******************************************************'
write(*,*)'*                                                     *' 
write(*,*)'* Program Transitivity: Version 1; 15.12.2013         *' 
write(*,*)'*                                                     *' 
write(*,*)'* Copyright Dr. Werner Ulrich                         *'
write(*,*)'*                                                     *'
write(*,*)'* The author does not take responsibility for correct *'
write(*,*)'* program run or any damages caused by the program.   *'
write(*,*)'*                                                     *'
write(*,*)'*******************************************************'
write(*,*)
return 
end

real*8 function norm(iseed)
real*8 a
a=0
do while(a.lt.0.000000001)
a=ran(iseed)
enddo
norm=sqrt(-2.0*dlog(a))*cos(2.0*3.14159*ran(iseed))
!write(*,*) norm
end function

subroutine matrandpairs(mat,row,sp,si,irep,imax,iseed)
integer sp,si,iseed,irep,imax
real mat(sp,si),row(sp)
mat=1.0
do 1 i=1,sp
do 2 j=1,si
if(i.ne.j)then
a=ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)
a=(a-6.0)/12.0
mat(i,j)=row(i)*(1+a)
endif
2 continue
1 continue
do 3 i=1,sp
do 4 j=i+1,si
a=mat(i,j)+mat(j,i)
mat(i,j)=mat(i,j)/a
mat(j,i)=1.0-mat(i,j)
4 continue
3 continue
return
end 

subroutine matrand(mat,row,sp,si,irep,imax,iseed)
integer sp,si,iseed,irep,imax
real mat(sp,si),row(sp),col(si)
mat=0
col=0
do 1 j=1,si
do 2 i=1,sp
a=ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)+ran(iseed)
a=(a-6.0)/12.0
mat(i,j)=row(i)*(1+a)
!mat(i,j)=ran(iseed)
col(j)=col(j)+mat(i,j)
2 continue
1 continue
do 3 j=1,si
do 4 i=1,sp
mat(i,j)=mat(i,j)/col(j)
4 continue
3 continue
return
end 

subroutine metric(mat,sp,n,k,deg,diag,f,dist)
integer sp,n,n1
real mat(sp,sp),temp(sp),deg,diag,f,dist
character*10 dir
dir='down'
temp=0
deg=0
diag=0
f=0
n=0
dist=0
k=0
do 10 j=1,sp
do 11 i=1,sp-1
do 12 i1=i+1,sp
if(i.ne.j.and.i1.ne.j)then
k=k+1
if(mat(i,j).lt.mat(i1,j))then
n=n+1
temp(i1)=1.0
dist=dist+(i1-i-1)
endif
endif
12 continue
11 continue
do 13 i=j+1,sp
13 if(mat(j,j).lt.mat(i,i))diag=diag+1
10 continue
do 14 i=1,sp
14 if(temp(i).gt.0)f=f+1
if(k.gt.0)deg=1.0-n/real(k)
if(n.gt.0)dist=dist/real(n)
if(dist.eq.0)dist=sp-1.0
if(sp.gt.1)diag=1.0-2.0*diag/real(sp*(sp-1))
f=f/real(sp)
return
end

subroutine metricpairs(mat,sp,n,k,deg,diag,f,dist)
integer sp,n,k
real mat(sp,sp),temp(sp),deg,diag,f,dist
real*8 a
character*10 dir
temp=0
deg=0
f=0
dist=0
n=0
a=0
diag=0
do 3 i=1,sp
do 4 i1=i+1,sp
if(mat(i,i1).gt.mat(i1,i))n=n+1
if(mat(i,i1).lt.mat(i1,i))temp(i)=1.0
4 continue
3 continue
do 5 i=1,sp
5 if(temp(i).gt.0)f=f+1
if(sp.gt.1)then
f=f/real(sp)
k=sp*(sp-1)/2
deg=n/real(k)
endif
do 40 i=1,sp-1
do 41 i1=i+1,sp
do 42 j=sp,2,-1
do 43 j1=j-1,1,-1
a=a+1
if(mat(i,j).gt.mat(i,j1).and.mat(i,j).gt.mat(i1,j).and.mat(i1,j1).lt.mat(i,j1).and.mat(i1,j1).lt.mat(i1,j))diag=diag+1
43 continue
42 continue
41 continue
40 continue
!write(*,*)diag,a,diag/a
if(a.gt.0)diag=diag/a
return
end

subroutine conflimnew (field,i,up,down,a,p)
use basic_stats
integer i,i1,i2
real field(i),up,down,a,p
i1=nint(0.05*i)
if(i1.lt.1)i1=1
i2=nint(0.95*i)
if(i2.gt.i)i2=i
call sort_vector_down(field,i)
up=field(i1)
down=field(i2)
p=0.5
do 1 i1=1,i
if(i1.lt.0.5*i)then
if(field(i1).lt.a)then
p=max(real(i1)/i,1.0/i)
return
endif
else
if(field(i1).le.a)then
p=max(1.0-real(i1)/i,1.0/i)
return
endif
endif
1 continue
if(a.lt.field(i))p=1.0/i
return
end subroutine conflimnew

subroutine ctop(c,p,sp,degree)
use mat_handling
use basic_stats
use matrix_basics
integer sp
real c(sp,sp),p(sp,sp),col(sp),temp(sp),a,a1,degree
p=0
col=0
do 1 i=1,sp
a=1
do 2 j=1,sp
2 a=a*c(i,j)
p(i,i)=a
1 continue
do 3 i=1,sp
do 4 j=1,sp
if(j.ne.i)then
temp=0
do 5 j1=1,sp
5 if(j1.ne.i.and.j1.ne.j)temp(j1)=c(i,j1)
a=geomean(temp,sp)
a1=0
do 6 k=0,sp-2
6 a1=a1+a**k
p(j,i)=c(j,i)*a1/(real(sp-1))
endif
4 continue
3 continue
a1=0
do 10 j=1,sp
a=0
do 11 i=1,sp
11 a=a+p(i,j)
a1=a1+abs(1-a)
10 continue
a1=a1/sp
!write(4,'(a,f9.6,1x,i3,1x,f9.6)')'Error ',degree,sp,a1
call colstandard(p,sp,sp)
return
end subroutine ctop

subroutine ptoc(pmat,interaction,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
use matrix_basics
integer sp,si,irep,number,iseed,nsim
real eps,epsor,b,degree
real transsim(nsim),diagsim(nsim),propsim(nsim),chlensim(nsim),interaction(sp,sp),cmat(sp,sp),pmat(sp,sp),trans(sp),trans1(sp),temp(sp,sp),ctemp(sp,sp),rowa(sp) 
real*8 out(sp,sp)
eps=epsor
10 ktest=0
transsim=0
diagsim=0
propsim=0
chlensim=0
do 1 irep=1,number
!write(*,*)'Replicate: ',irep,'Maximum: ',number
call matrandpairs(ctemp,rowa,sp,sp,irep,number,iseed)
call ctop(ctemp,temp,sp,degree)
spear=0
b=eps
do 2 j=1,sp
do 3 i=1,sp
trans(i)=pmat(i,j)
trans1(i)=temp(i,j)
3 continue
spear=spear+spearman(trans,trans1,sp)
2 continue
spear=spear/real(si)
if(spear.gt.eps)then
if(spear.gt.b)then
b=spear
interaction=ctemp
endif
ktest=ktest+1
if(ktest.gt.nsim)goto 5
call metricpairs(ctemp,sp,n,npairs,transmetric,diag,prop,chlen)
transsim(ktest)=transmetric
diagsim(ktest)=diag
propsim(ktest)=prop
chlensim(ktest)=chlen
endif
1 continue
if(ktest.lt.nsim)then
write(*,*)'No convergence, reducing convergence limit'
eps=eps-0.05
goto 10
else
ktest=nsim
endif
5 return
end subroutine ptoc

subroutine mtop(mat,transition,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
use matrix_basics
integer sp,si,irep,number,iseed,nsim
real eps,epsor,b
real transsim(nsim),diagsim(nsim),propsim(nsim),chlensim(nsim),transition(sp,sp),mat(sp,si),trans(sp),trans1(sp),ptemp(sp,sp),rowa(sp) 
real*8 out(sp,si)
eps=epsor
10 ktest=0
transsim=0
diagsim=0
propsim=0
chlensim=0
do 1 irep=1,number
!write(*,*)'Replicate: ',irep,'Maximum: ',number
call matrand(ptemp,rowa,sp,sp,irep,number,iseed)
call dotproduct (ptemp,mat,out,sp,si,sp) 
spear=0
b=eps
do 2 j=1,si
do 3 i=1,sp
trans(i)=mat(i,j)
trans1(i)=real(out(i,j))
3 continue
spear=spear+spearman(trans,trans1,sp)
2 continue
spear=spear/real(si)
if(spear.gt.eps)then
if(spear.gt.b)then
b=spear
transition=ptemp
endif
ktest=ktest+1
if(ktest.gt.nsim)goto 5
call metric(ptemp,sp,n,npairs,transmetric,diag,prop,chlen)
transsim(ktest)=transmetric
diagsim(ktest)=diag
propsim(ktest)=prop
chlensim(ktest)=chlen
endif
1 continue
if(ktest.lt.nsim)then
write(*,*)'No convergence, reducing convergence limit'
eps=eps-0.05
goto 10
else
ktest=nsim
endif
5 return
end subroutine mtop

subroutine mtoc(mat,interaction,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
use matrix_basics
integer sp,si,irep,number,iseed,nsim
real eps,epsor,b,degree
real transsim(nsim),diagsim(nsim),propsim(nsim),chlensim(nsim),interaction(sp,sp),mat(sp,si),trans(sp),trans1(sp),ptemp(sp,sp),ctemp(sp,sp),rowa(sp) 
real*8 out(sp,si)
eps=epsor
10 ktest=0
transsim=0
diagsim=0
propsim=0
chlensim=0
do 1 irep=1,number
!write(*,*)'Replicate: ',irep,'Maximum: ',number
call matrandpairs(ctemp,rowa,sp,sp,irep,number,iseed)
call ctop(ctemp,ptemp,sp,degree)
call dotproduct (ptemp,mat,out,sp,si,sp) 
spear=0
b=eps
do 2 j=1,si
do 3 i=1,sp
trans(i)=mat(i,j)
trans1(i)=real(out(i,j))
3 continue
spear=spear+spearman(trans,trans1,sp)
2 continue
spear=spear/real(si)
if(spear.gt.eps)then
if(spear.gt.b)then
b=spear
interaction=ctemp
endif
ktest=ktest+1
if(ktest.gt.nsim)goto 5
call metricpairs(ctemp,sp,n,npairs,transmetric,diag,prop,chlen)
transsim(ktest)=transmetric
diagsim(ktest)=diag
propsim(ktest)=prop
chlensim(ktest)=chlen
endif
1 continue
if(ktest.lt.nsim)then
write(*,*)'No convergence, reducing convergence limit'
eps=eps-0.05
goto 10
else
ktest=nsim
endif
5 return
end subroutine mtoc

subroutine ttop(mat,mat1,transition,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
use matrix_basics
integer sp,si,irep,number,iseed,nsim
real eps,epsor,b
real transsim(nsim),diagsim(nsim),propsim(nsim),chlensim(nsim),transition(sp,sp),mat(sp,si-1),mat1(sp,si-1),mat2(sp,si-1),trans(sp),trans1(sp),ptemp(sp,sp),rowa(sp) 
eps=epsor
10 ktest=0
transsim=0
diagsim=0
propsim=0
chlensim=0
do 1 irep=1,number
!write(*,*)'Replicate: ',irep,'Maximum: ',number
call matrand(ptemp,rowa,sp,sp,irep,number,iseed)
call dotrealproduct (ptemp,mat,mat2,sp,si-1,sp) 
call colstandard(mat2,sp,si-1)
spear=0
b=eps
do 2 j=1,si-1
do 3 i=1,sp
trans(i)=mat1(i,j)
trans1(i)=mat2(i,j)
3 continue
spear=spear+spearman(trans,trans1,sp)
2 continue
spear=spear/real(si-1)
if(spear.gt.eps)then
if(spear.gt.b)then
b=spear
transition=ptemp
endif
ktest=ktest+1
if(ktest.gt.nsim)goto 5
call metric(ptemp,sp,n,npairs,transmetric,diag,prop,chlen)
transsim(ktest)=transmetric
diagsim(ktest)=diag
propsim(ktest)=prop
chlensim(ktest)=chlen
endif
1 continue
if(ktest.lt.nsim)then
write(*,*)'No convergence, reducing convergence limit'
eps=eps-0.05
goto 10
else
ktest=nsim
endif
5 return
end subroutine ttop

subroutine ttoc(mat,mat1,transition,rowa,transsim,diagsim,propsim,chlensim,epsor,eps,sp,si,irep,number,nsim,iseed)
use matrix_basics
integer sp,si,irep,number,iseed,nsim
real eps,epsor,b
real transsim(nsim),diagsim(nsim),propsim(nsim),chlensim(nsim),transition(sp,sp),mat(sp,si-1),mat1(sp,si-1),mat2(sp,si-1),trans(sp),trans1(sp),ptemp(sp,sp),ctemp(sp,sp),rowa(sp) 
eps=epsor
10 ktest=0
transsim=0
diagsim=0
propsim=0
chlensim=0
do 1 irep=1,number
!write(*,*)'Replicate: ',irep,'Maximum: ',number
call matrandpairs(ctemp,rowa,sp,sp,irep,number,iseed)
call ctop(ctemp,ptemp,sp,degree)
call dotrealproduct (ptemp,mat,mat2,sp,si-1,sp) 
call colstandard(mat2,sp,si-1)
spear=0
b=eps
do 2 j=1,si-1
do 3 i=1,sp
trans(i)=mat1(i,j)
trans1(i)=mat2(i,j)
3 continue
spear=spear+spearman(trans,trans1,sp)
2 continue
spear=spear/real(si-1)
if(spear.gt.eps)then
if(spear.gt.b)then
b=spear
transition=ctemp
endif
ktest=ktest+1
if(ktest.gt.nsim)goto 5
call metricpairs(ctemp,sp,n,npairs,transmetric,diag,prop,chlen)
transsim(ktest)=transmetric
diagsim(ktest)=diag
propsim(ktest)=prop
chlensim(ktest)=chlen
endif
1 continue
if(ktest.lt.nsim)then
write(*,*)'No convergence, reducing convergence limit'
eps=eps-0.05
goto 10
else
ktest=nsim
endif
5 return
end subroutine ttoc

module matrix_metrics
use basic_stats
use matrix_basics
use eigensolving
use sorting
contains

subroutine metrics (mat1,sp,si,metric,method,neigen)
integer sp,si,neigen
real embabs,absences,roweigen,coleigen,r2,tog,nodf,replace,clumping,discrepancy,cturn,csegr
real morisita,chao,cscore,wtog,wcs,wturn,wsegr
real mat1(sp,si),metric
character*50,null,method
if(method.eq.'roweigen')metric=roweigen(mat1,sp,si,neigen)
if(method.eq.'coleigen')metric=coleigen(mat1,sp,si,neigen)
if(method.eq.'absences')metric=absences(mat1,sp,si)
if(method.eq.'r2')metric=r2(mat1,sp,si)
if(method.eq.'tog')metric=tog(mat1,sp,si)
if(method.eq.'tog1')metric=tog1(mat1,sp,si)
if(method.eq.'embabs')metric=embabs(mat1,sp,si)
if(method.eq.'nodf')metric=nodf(mat1,sp,si)
if(method.eq.'replace')metric=replace(mat1,sp,si)
if(method.eq.'clumping')metric=clumping(mat1,sp,si)
if(method.eq.'discrepancy')metric=discrepancy(mat1,sp,si)
if(method.eq.'morisita')metric=morisita(mat1,sp,si)
if(method.eq.'chao')metric=chao(mat1,sp,si)
if(method.eq.'cscore')metric=cscore(mat1,sp,si)
if(method.eq.'cscorepair')metric=cscorepair(mat1,sp,si)
if(method.eq.'cscorevar')metric=cscorevar(mat1,sp,si)
if(method.eq.'cscoreskew')metric=cscoreskew(mat1,sp,si)
if(method.eq.'joccvar')metric=joccvar(mat1,sp,si)
if(method.eq.'cturn')metric=cturn(mat1,sp,si)
if(method.eq.'csegr')metric=csegr(mat1,sp,si)
if(method.eq.'wtog')metric=wtog(mat1,sp,si)
if(method.eq.'wcs')metric=wcs(mat1,sp,si)
if(method.eq.'test')metric=test(mat1,sp,si)
return
end subroutine metrics

real function nodf (mat1,sp,si)
integer sp,si
real mat1(sp,si),row(sp),col(si),a
row=0
col=0
do 1 i=1,sp
do 2 j=1,si
if(mat1(i,j).gt.0)then
row(i)=row(i)+1
col(j)=col(j)+1
endif
2 continue
1 continue
do 3 i=1,sp-1
do 4 i1=i+1,sp
a=0
do 5 j=1,si
if(mat1(i,j).gt.0.and.mat1(i1,j).gt.0)a=a+1
5 continue
if(row(i1).gt.0.and.row(i1).lt.row(i))then
a=a/row(i1)
nodf=nodf+a
endif
4 continue
3 continue
do 6 j=1,si-1
do 7 j1=j+1,si
a=0
do 8 i=1,sp
if(mat1(i,j).gt.0.and.mat1(i,j1).gt.0)a=a+1
8 continue
if(col(j1).gt.0.and.col(j1).lt.col(j))then
a=a/col(j1)
nodf=nodf+a
endif
7 continue
6 continue
nodf=2.0*nodf/(sp*(sp-1)+si*(si-1))
end function nodf

real function discrepancy (mat1,sp,si)
integer sp,si,row(sp),col(si)
real mat1(sp,si),br,br1,fill
row=0
col=0
fill=0
do 1 i=1,sp
do 2 j=1,si
if(mat1(i,j).gt.0)then
row(i)=row(i)+1
col(j)=col(j)+1
fill=fill+1
endif
2 continue
1 continue
br=0
br1=0
do 10 i=1,sp
do 11 j=1,row(i)
if(mat1(i,j).eq.0)br=br+1
11 continue
10 continue
do 20 j=1,si
do 21 i=1,col(j)
if(mat1(i,j).eq.0)br1=br1+1
21 continue
20 continue
if(br1.lt.br)br=br1
discrepancy=br/real(fill)
end function discrepancy

real function tog(mat,sp,si)
integer sp,si
real mat(sp,si)
tog=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
do 3 j=1,si-1
do 4 j1=j+1,si
if((mat(i,j).gt.0.and.mat(i1,j1).eq.0.and.mat(i1,j).gt.0.and.mat(i,j1).eq.0).or.&
(mat(i,j).eq.0.and.mat(i1,j1).gt.0.and.mat(i1,j).eq.0.and.mat(i,j1).gt.0))tog=tog+1
4 continue
3 continue
2 continue
1 continue
tog=tog/(sp*(sp-1)*si*(si-1)/4)
end function tog

real function cscore(mat,sp,si)
integer sp,si
real mat(sp,si)
cscore=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
do 3 j=1,si-1
do 4 j1=j+1,si
if((mat(i,j).gt.0.and.mat(i1,j1).gt.0.and.mat(i1,j).eq.0.and.mat(i,j1).eq.0).or.&
(mat(i,j).eq.0.and.mat(i1,j1).eq.0.and.mat(i1,j).gt.0.and.mat(i,j1).gt.0))cscore=cscore+1
4 continue
3 continue
2 continue
1 continue
cscore=cscore/(sp*(sp-1)*si*(si-1)/4)
end function cscore

real function cscorepair(mat,sp,si)
integer sp,si,count,row(sp)
real mat(sp,si)
cscorepair=0
row=0
do 10 i=1,sp
do 11 j=1,si
11 if(mat(i,j).gt.0)row(i)=row(i)+1
10 continue
do 1 i=1,sp-1
do 2 i1=i+1,sp
count=0
do 3 j=1,si
if(mat(i,j).eq.1.and.mat(i1,j).eq.1)count=count+1
3 continue
count=(row(i)-count)*(row(i1)-count)
cscorepair=cscorepair+count
2 continue
1 continue
cscorepair=cscorepair/(sp*(sp-1)/2)
end function cscorepair

real function cscorevar(mat,sp,si)
integer sp,si,count,row(sp)
real mat(sp,si),cs
cscorevar=0
cs=0
row=0
do 10 i=1,sp
do 11 j=1,si
11 if(mat(i,j).gt.0)row(i)=row(i)+1
10 continue
do 1 i=1,sp-1
do 2 i1=i+1,sp
count=0
do 3 j=1,si
if(mat(i,j).eq.1.and.mat(i1,j).eq.1)count=count+1
3 continue
count=(row(i)-count)*(row(i1)-count)
cs=cs+count
cscorevar=cscorevar+count*count
2 continue
1 continue
cs=cs/(sp*(sp-1)/2)
cscorevar=cscorevar/(sp*(sp-1)/2)-cs*cs
end function cscorevar

real function cscoreskew(mat,sp,si)
integer sp,si,count,row(sp)
real mat(sp,si),a,b
cscoreskew=0
row=0
do 10 i=1,sp
do 11 j=1,si
11 if(mat(i,j).gt.0)row(i)=row(i)+1
10 continue
do 1 i=1,sp-1
do 2 i1=i+1,sp
count=0
do 3 j=1,si
if(mat(i,j).eq.1.and.mat(i1,j).eq.1)count=count+1
3 continue
count=(row(i)-count)*(row(i1)-count)
cscoreskew=cscoreskew+count*count*count
2 continue
1 continue
a=cscorepair(mat,sp,si)
b=cscorevar(mat,sp,si)
cscoreskew=cscoreskew/(sp*(sp-1)/2)-3*a*b-a*a*a
if(b.gt.0)b=sqrt(b)
cscoreskew=cscoreskew/(b*b*b)
end function cscoreskew

real function joccvar(mat,sp,si)
integer sp,si,count
real mat(sp,si),row(sp*(sp-1)/2)
joccvar=0
row=0
k=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
k=k+1
count=0
do 3 j=1,si
if(mat(i,j).eq.1.and.mat(i1,j).eq.1)count=count+1
3 continue
row(k)=count
2 continue
1 continue
joccvar=variance(row,sp*(sp-1)/2)
joccvar=joccvar/(sp*(sp-1)*si*(si-1)/4)
end function joccvar

real function test(mat,sp,si)
integer sp,si
real mat(sp,si),k2,km,kn,count,row(sp)
test=0
row=0
k2=0
km=0
kn=0
do 10 i=1,sp
do 11 j=1,si
11 if(mat(i,j).gt.0)row(i)=row(i)+1
10 continue
do 1 i=1,sp-1
do 2 i1=i+1,sp
count=0
do 3 j=1,si
if(mat(i,j).eq.1.and.mat(i1,j).eq.1)count=count+1
3 continue
k2=k2+0.5*count*count
km=km+count*(row(i)+row(i1)-0.5)
kn=kn+row(i)*row(i1)
2 continue
1 continue
test=k2-km
test=test/(sp*(sp-1)*si*(si-1)/4)
end function test

real function clumping(mat,sp,si)
integer sp,si
real mat(sp,si)
clumping=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
do 3 j=1,si-1
do 4 j1=j+1,si
if(mat(i,j).gt.0.and.mat(i,j1).gt.0.and.mat(i1,j).gt.0.and.mat(i1,j1).gt.0)clumping=clumping+1
4 continue
3 continue
2 continue
1 continue
clumping=clumping/(sp*(sp-1)*si*(si-1)/4)
end function clumping

real function tog1(mat,sp,si)
integer sp,si
real mat(sp,si)
tog1=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
do 3 j=1,si-1
do 4 j1=j+1,si
if((mat(i,j).gt.0.and.mat(i1,j1).eq.0.and.mat(i1,j).eq.0.and.mat(i,j1).gt.0).or.&
(mat(i,j).eq.0.and.mat(i1,j1).gt.0.and.mat(i1,j).gt.0.and.mat(i,j1).eq.0))tog1=tog1+1
4 continue
3 continue
2 continue
1 continue
tog1=tog1/(sp*(sp-1)*si*(si-1)/4)
end function tog1

real function absences(mat,sp,si)
integer sp,si
real mat(sp,si)
absences=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
do 3 j=1,si-1
do 4 j1=j+1,si
if(mat(i,j).eq.0.and.mat(i,j1).eq.0.and.mat(i1,j).eq.0.and.mat(i1,j1).eq.0)absences=absences+1
4 continue
3 continue
2 continue
1 continue
absences=absences/(sp*(sp-1)*si*(si-1)/4)
end function absences

real function embabs(mat1,sp,si)
integer sp,si
real mat1(sp,si),row(sp),fill
row=0
fill=0
embabs=0
do 1 i=1,sp
do 2 j=1,si
if(mat1(i,j).gt.0)then
row(i)=row(i)+1
fill=fill+1
endif
2 continue
1 continue
do 10 i=1,sp
do 11 j1=1,si
if(mat1(i,j1).gt.0)goto 12
11 continue
12 do 13 j2=si,1,-1
if(mat1(i,j2).gt.0)goto 14
13 continue
goto 10
14 embabs=embabs+(j2-j1+1)-row(i)
10 continue
if(fill.gt.0)embabs=embabs/fill
if(embabs.lt.0)embabs=0
end function embabs

real function r2(mat1,sp,si)
integer sp,si
real mat1(sp,si)
real,allocatable:: s(:),s1(:)
k=0
do 1 j=1,si
do 2 i=1,sp
if(mat1(i,j).gt.0)then
k=k+1
endif
2 continue
1 continue
allocate(s(k),s1(k))
k=0
do 3 j=1,si
do 4 i=1,sp
if(mat1(i,j).gt.0)then
k=k+1
s(k)=i
s1(k)=j
endif
4 continue
3 continue
r2=pearson0(s,s1,k)
r2=r2*r2
end function r2

real function replace(mat1,sp,si)
integer sp,si
real mat1(sp,si)
replace=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
do 3 j=1,si
if((mat1(i,j).gt.0.and.mat1(i1,j).eq.0).or.(mat1(i,j).eq.0.and.mat1(i1,j).gt.0))replace=replace+1 
3 continue
2 continue
1 continue
a=sp*si*(sp-1)*(si-1)/4.0
replace=replace/a
end function replace

real function csegr(mat1,sp,si)
integer sp,si,fill
real mat1(sp,si)
fill=0
csegr=0
do 1 j=1,si
do 2 i=1,sp
if(mat1(i,j).gt.0)then
fill=fill+1
endif
2 continue
1 continue
if(sp.ge.si)then
jwide=fill/real(2*si)
else
iwide=fill/real(2*sp)
endif
do 10 i=1,sp-1
do 11 i1=i+1,sp
do 12 j=1,si-1
do 13 j1=j+1,si
if(sp.ge.si)then
j2=int(real(si)/real(sp)*(i-1))+1
j3=int(real(si)/real(sp)*(i1-1))+1
if((j.lt.j2-jwide.or.j.gt.j2+jwide).and.(j1.lt.j3-jwide.or.j1.gt.j3+jwide))then
if(mat1(i,j).gt.0.and.mat1(i1,j1).gt.0.and.mat1(i1,j).eq.0.and.mat1(i,j1).eq.0)csegr=csegr+1
endif
else
i2=int(real(sp)/real(si)*(j-1))+1
i3=int(real(sp)/real(si)*(j1-1))+1
if((i.lt.i2-iwide.or.i.gt.i2+iwide).and.(i1.lt.i3-iwide.or.i1.gt.i3+iwide))then
if(mat1(i,j).gt.0.and.mat1(i1,j1).gt.0.and.mat1(i1,j).eq.0.and.mat1(i,j1).eq.0)csegr=csegr+1
endif
endif
13 continue
12 continue
11 continue
10 continue
a=sp*si*(sp-1)*(si-1)/4.0
csegr=csegr/a
end function csegr

real function cturn(mat1,sp,si)
integer sp,si,fill
real mat1(sp,si)
fill=0
cturn=0
do 1 j=1,si
do 2 i=1,sp
if(mat1(i,j).gt.0)then
fill=fill+1
endif
2 continue
1 continue
if(sp.ge.si)then
jwide=fill/real(2*si)
else
iwide=fill/real(2*sp)
endif
do 10 i=1,sp-1
do 11 i1=i+1,sp
do 12 j=1,si-1
do 13 j1=j+1,si
if(sp.ge.si)then
j2=int(real(si)/real(sp)*(i-1))+1
j3=int(real(si)/real(sp)*(i1-1))+1
if((j.ge.j2-jwide.and.j.le.j2+jwide).and.(j1.ge.j3-jwide.and.j1.le.j3+jwide))then
if(mat1(i,j).gt.0.and.mat1(i1,j1).gt.0.and.mat1(i1,j).eq.0.and.mat1(i,j1).eq.0)cturn=cturn+1
endif
else
i2=int(real(sp)/real(si)*(j-1))+1
i3=int(real(sp)/real(si)*(j1-1))+1
if((i.ge.i2-iwide.and.i.le.i2+iwide).and.(i1.ge.i3-iwide.and.i1.le.i3+iwide))then
if(mat1(i,j).gt.0.and.mat1(i1,j1).gt.0.and.mat1(i1,j).eq.0.and.mat1(i,j1).eq.0)cturn=cturn+1
endif
endif
13 continue
12 continue
11 continue
10 continue
a=sp*si*(sp-1)*(si-1)/4.0
cturn=cturn/a
end function cturn

real function wcs(mat1,sp,si)
integer sp,si
real mat1(sp,si)
do 10 i=1,sp-1
do 11 i1=i+1,sp
do 12 j=1,si-1
do 13 j1=j+1,si
if((mat1(i,j).gt.mat1(i,j1).and.mat1(i,j).gt.mat1(i1,j).and.mat1(i1,j1).gt.mat1(i,j1).and.mat1(i1,j1).gt.mat1(i1,j)).or.&
(mat1(i,j).lt.mat1(i,j1).and.mat1(i,j).lt.mat1(i1,j).and.mat1(i1,j1).lt.mat1(i,j1).and.mat1(i1,j1).lt.mat1(i1,j)))wcs=wcs+1
13 continue
12 continue
11 continue
10 continue
a=sp*si*(sp-1)*(si-1)/4.0
wcs=wcs/a
end function wcs

real function wtog(mat1,sp,si)
integer sp,si
real mat1(sp,si)
do 10 i=1,sp-1
do 11 i1=i+1,sp
do 12 j=1,si-1
do 13 j1=j+1,si
if((mat1(i,j).gt.mat1(i,j1).and.mat1(i1,j).gt.mat1(i1,j1)).or.&
(mat1(i,j).lt.mat1(i,j1).and.mat1(i1,j).lt.mat1(i1,j1)))wtog=wtog+1
13 continue
12 continue
11 continue
10 continue
a=sp*si*(sp-1)*(si-1)/4.0
wtog=wtog/a
end function wtog

real function simpson (mat1,sp,si)
integer sp,si
real mat1(sp,si)
simpson=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
a=0
b=0
c=0
do 3 j=1,si
if(mat1(i,j).gt.0.and.mat1(i1,j).eq.0)a=a+1
if(mat1(i,j).eq.0.and.mat1(i1,j).gt.0)b=b+1
if(mat1(i,j).gt.0.and.mat1(i1,j).gt.0)c=c+1
3 continue
simpson=simpson+c/(a+b)
2 continue
1 continue
end function simpson

real function soerensen (mat1,sp,si)
integer sp,si
real mat1(sp,si)
soerensen=0
do 1 i=1,sp-1
do 2 i1=i+1,sp
a=0
b=0
c=0
do 3 j=1,si
if(mat1(i,j).gt.0.and.mat1(i1,j).eq.0)a=a+1
if(mat1(i,j).eq.0.and.mat1(i1,j).gt.0)b=b+1
if(mat1(i,j).gt.0.and.mat1(i1,j).gt.0)c=c+1
3 continue
soerensen=soerensen+2.0*c/(a+b+2.0*c)
2 continue
1 continue
end function soerensen

real function roweigen(mat1,sp,si,neigen)
integer sp,si,neigen
real mat1(sp,si),mat(sp,sp),mat2(si,sp)
real*8 eigenvector(sp,sp),eigenvalue(sp)
mat2=transpose(mat1)
!if(neigen.eq.0)call distance(mat2,mat,si,sp)
if(neigen.eq.0)call distsoerensen(mat2,mat,si,sp)
!call fieldwrite(mat,sp,sp,4)
if(neigen.eq.1)call distcorr(mat2,mat,si,sp)
call powereigen(mat,sp,1,itest,eigenvector,eigenvalue)
roweigen=real(eigenvalue(1))
end function roweigen

real function coleigen(mat1,sp,si,neigen)
integer sp,si,neigen
real mat1(sp,si),mat(si,si)
real*8 eigenvector(si,si),eigenvalue(si)
!if(neigen.eq.0)call distance(mat1,mat,sp,si)
if(neigen.eq.0)call distsoerensen(mat1,mat,sp,si)
if(neigen.eq.1)call distcorr(mat1,mat,sp,si)
call powereigen(mat,si,1,itest,eigenvector,eigenvalue)
coleigen=real(eigenvalue(1))
end function coleigen

real function morisita(mat1,sp,si)
integer sp,si,col(si)
real mat1(sp,si)
a=0
a1=0
col=0
do 60 i=1,sp
do 61 j=1,si
if(mat1(i,j).gt.0)then
col(j)=col(j)+1
a1=a1+1
goto 60
endif
61 continue
60 continue
do 62 i=1,sp
do 63 j=si,1,-1
if(mat1(i,j).gt.0)then
col(j)=col(j)+1
a1=a1+1
goto 62
endif
63 continue
62 continue
do 65 j=1,si
if(a1.gt.2)a=a+col(j)*(col(j)-1)/(a1*(a1-1))
65 continue
morisita=si*a
end function morisita

real function chao (mat1,sp,si)
integer sp,si,col(si)
real mat1(sp,si),mat(sp,si)
chao=0
mat=0
col=0
do 1 i=1,sp
do 2 j=1,si
if(mat1(i,j).gt.0)col(j)=col(j)+1
2 continue
1 continue
do 3 j=1,si
if(col(j).gt.0)then
do 4 i=1,sp
mat(i,j)=mat1(i,j)/col(j)
4 continue
endif
3 continue
a1=0
do 5 i=1,sp
a=0
do 6 j=1,si
a=a+mat(i,j)
6 continue
a=a*a
b=0
do 7 j=1,si
b=b+mat(i,j)*mat(i,j)
7 continue
chao=chao+a-b
a1=a1+b
5 continue
if(a1.gt.0.and.si.gt.1)then
chao=(chao/a1)/(si-1)
else
chao=0
endif
end function chao

subroutine roweigendist(mat1,sp,si,ev,neigen)
integer sp,si,neigen
real mat1(sp,si),mat(sp,sp),mat2(si,sp),ev(sp)
real*8 eigenvector(sp,sp),eigenvalue(sp)
character*10 dir
dir='down'
mat2=transpose(mat1)
if(neigen.eq.0)call distance(mat2,mat,si,sp)
if(neigen.eq.0)call distsoerensen(mat2,mat,si,sp)
if(neigen.eq.1)call distcorr(mat2,mat,si,sp)
call symeigen(mat,sp,eigenvector,eigenvalue,itest)
ev=real(eigenvalue)
call sort_vector(ev,sp,dir)
ev=ev
end subroutine roweigendist

end module matrix_metrics

module eigensolving
use matrix_basics
contains
!********************************************************
!* Eigenvalues and eigenvectors of a real square matrix *
!* by Rutishauser's method and inverse iteration method *
!* ---------------------------------------------------- *
!* Reference:                                           *
!*                                                      *
!*   "ALGEBRE Algorithmes et programmes en Pascal       *
!*    de Jean-Louis Jardrin - Dunod BO-PRE 1988."       *
!*    [BIBLI 10].                                       *
!*                                                      *
!*                    F90 Release By J-P Moreau, Paris. *
!* ---------------------------------------------------- *

  !***************************************************************
  !* Subroutine DECCRM determines the lower triangular matrix and*
  !* the upper triangukar matrix of Crout's decomposition of a   *
  !* given square real matrix, A.                                *
  !* ----------------------------------------------------------- *
  !* INPUTS:                                                     *
  !*         eps: required precision (double)                    *
  !*          n : size of matrix A (integer)                     *
  !*          A : input matrix (n x n)                           *
  !* OUTPUTS:                                                    *
  !*          it: flag, =0 if method does not apply              *
  !*                    =1 if method is ok.                      *
  !*           U: lower triangular matrix.                       *
  !*           V: upper triangular matrix.                       *
  !***************************************************************
  Subroutine DECCRM(eps, n, A, it, U, V)
    real*8 A(n,n), U(n,n), V(n,n)
    real*8 eps, s

    if (dabs(A(1,1)) < eps)  then
	  it=0
    else
      do i=1, n
	    U(i,1)=A(i,1)
      end do 
      V(1,1)=1.d0
      do j=2, n
	    V(1,j)=A(1,j)/U(1,1)
      end do
      it=1; k=2
      do while (it.ne.0.and.k<=n)
	    do i=1, n
	      if (i < k) then  
		    U(i,k)=0.d0
		  else
	        s=0.d0
	        do j=1, k-1
			  s = s + U(i,j)*V(j,k)
            end do
	        U(i,k) = A(i,k) - s
		  end if
        end do
	    if (dabs(U(k,k)) < eps) then 
          it=0
		else
	      do j=1, n
	        if (j < k) then 
			  V(k,j)=0.d0
	        else if (j==k) then
			  V(k,j)=1.d0
            else
	          s=0.d0
	          do i=1, k-1
			    s = s + U(k,i)*V(i,j)
              end do
	          V(k,j) = A(k,j)/U(k,k);
			end if
          end do
	      k = k + 1
		end if
      end do !while
    end if
	return
  End subroutine DECCRM

  !*********************************************************
  !* Calculate the eigenvalues of a real square matrix by  *
  !* Rutishauser's Method.                                 *
  !* ----------------------------------------------------- *
  !* INPUTS:                                               *
  !*        eps: absolute precision (double)               *
  !*        dta: relative precision (double)               *
  !*         m : maximum number of iterations (integer)    *
  !*         n : size of matrix A (integer)                *
  !*         A : input real square matrix (n x n)          *
  !* OUTPUTS:                                              *
  !*         it: flag, =-1 if convergence is not obtained  *
  !*                   =1 if convergence is ok.            *
  !*         R : contains in output the n eigenvalues of A *
  !*                                                       *         
  !*********************************************************
  Subroutine VAMR(eps, dta, m, n, A, it, R)
    real*8 A(n,n), R(n)
    real*8 eps, dta
    real*8 phi,s,t0
    real*8 U(n,n), V(n,n)

    t0=0.d0
    l=1
    do while (l<=m.and.it.ne.1)
      do i=1, n
	    R(i)=A(i,i)
      end do

      call DECCRM(eps, n, A, it, U, V)

      if (it==0) then
	    do i=1, n
		  A(i,i)=A(i,i) + 1.d0
        end do
	    t0 = t0 + 1.d0
      else
	    do i=1, n
		  do j=1, n
	        s=0.d0
	        do k=1, n
			  s = s + V(i,k) * U(k,j)
            end do
	        A(i,j) = s
		  end do
        end do
	    phi=0.d0
        do i=1, n
	      s= dabs(A(i,i)-R(i))
	      if (s > phi)  phi=s
		end do
	    if (phi < dta) then
	      do i=1, n
		    R(i) = A(i,i) - t0
          end do
		else
	      l = l + 1
	      it=-1
		end if
      end if
    end do !while
	return
  End subroutine VAMR

  !************************************************************
  !* Procedure IIM calculates a real eigenvalue and the asso- *
  !* ciated eigenvector of a real square matrix the inverse   *
  !* iteration method.                                        *
  !* -------------------------------------------------------- *
  !* INPUTS:                                                  *
  !*         eps : absolute precision (double)                *
  !*         dta : relative precision (double)                *
  !*          m  : maximum number of iterations (integer)     *
  !*          n  : size of matrix A                           *
  !*          A  : input real square matrix (n x n)           *
  !* OUTPUTS:                                                 *
  !*          it : flag, =-1 if convergence is not obtained   *
  !*                     =1 if convergence is ok.             *
  !*        Gamma: starting value for the eigenvalue as input *
  !*               approximation of the eigenvalue with preci-*
  !*               sion dta in output.                        *
  !*          X1 : contains in output the associated eigen-   *
  !*               vector.                                    *
  !*                                                          *
  !************************************************************
  Subroutine IIM(eps, dta, m, n, A, it, gamma, X1)
    real*8 A(n,n), X1(n)
    real*8 eps, dta, gamma
    real*8  p0,phi,s,t0
    real*8 W(n), X0(n)
    integer LP(n)

    do i=1, n
	  A(i,i) = A(i,i) - gamma
    end do
    do k=1, n-1
      p0=A(k,k); l0=k
      do i=k+1, n
		if (dabs(A(i,k)) > dabs(p0)) then
	      p0=A(i,k); l0=i
		end if
      end do 
      LP(k)=l0
      if (dabs(p0) < eps) then
	    p0=eps; A(l0,k)=eps
      end if
      if (l0.ne.k) then
	    do j=k, n
	      t0=A(k,j); A(k,j)=A(l0,j); A(l0,j)=t0
		end do
      end if
	  do i=k+1, n
	    A(i,k)=A(i,k)/p0
	    do j=k+1, n
	      A(i,j)=A(i,j)-A(i,k)*A(k,j)
        end do
      end do
	end do !k loop

    if (dabs(A(n,n)) < eps)  A(n,n)=eps
    do i=1, n
	  X0(i)=1.d0/sqrt(1.d0*i)
    end do 
    it=-1; l=1
    do while (it==-1.and.l<=m)
      do i=1, n
	    W(i)=X0(i)
      end do
      do k=1, n-1
        l0=LP(k)
        if (l0.ne.k) then
          t0=W(k); W(k)=W(l0); W(l0)=t0
        end if
        do i=k+1, n
		  W(i)=W(i)-A(i,k)*W(k)
        end do
      end do
      X1(n)=W(n)/A(n,n)
      do i=n-1, 1, -1
        s=0.d0
        do j=i+1, n
		  s = s + A(i,j)*X1(j)
        end do
	    X1(i)=(W(i)-s)/A(i,i)
      end do
      p0=0.d0
      do i=1, n
        if (dabs(X1(i)) > dabs(p0))  p0=X1(i)
      end do
      do i=1, n
        X1(i)=X1(i)/p0
      end do
      phi=0.d0
      do i=1, n
        s=dabs(X1(i)-X0(i))
        if (s > phi)  phi=s
      end do

      if (phi < dta) then
        gamma = gamma + 1.d0/p0
        it=1
      else
        do i=1, n
		  X0(i)=X1(i)
        end do 
        l = l + 1
      end if
	end do !while
	return
  End subroutine IIM !IIM  	  						    


  !******************************************************
  !* INPUTS:                                            *
  !* EPS : precision (Double)                           *
  !* D1  : precision d1 (Double)                        *
  !* D2  : precision d2 (Double)                        *
  !* M   : maximum number of iterations (integer)       *
  !* N   : order of matrix A (integer)                  *
  !* A   : input matrix to study (of MAT type)          *
  !* -------------------------------------------------- *
  !* OUTPUTS:                                           *
  !* IT  : -1 if no convergence is obtained (integer)   *
  !* R   : table of eigenvalues (of VEC type)           *
  !* VX  : table of eigenvectors (of MAT type)          *
  !******************************************************
  Subroutine eigen(n, A, it, R, VX)
    real*8 A(n,n), R(n), VX(n,n),eps,d1,d2
    real*8 X(n), A1(n,n)
    eps=1.d-10
    d1=1.d-8
    d2=1.d-8
    m=1000

    do i=1, n
      do j=1, n
        A1(i,j) = A(i,j)
      end do
    end do

    call VAMR(eps,d2,m,n,A1,it,R)

! restore A1 after VAMR
    do i=1, n
      do j=1, n
        A1(i,j) = A(i,j)
      end do
	end do   

    j=1
    do while (it==1.and.j<=n)

    call IIM(eps,d1,m,n,A1,it,R(j),X)
  
! restore A1 after IIM
      do i=1, n
        do k=1, n
          A1(i,k) = A(i,k)
        end do
        VX(i,j) = X(i)
      end do
      j = j + 1     
    end do !while
	
	return
  End subroutine eigen

Subroutine symeigen(mat,N,V,D,NROT)
integer N,NROT
real mat(N,N)
real*8  A(1:N,1:N),V(1:N,1:N),D(1:N)
real*8, pointer :: B(:), Z(:)
real*8  c,g,h,s,sm,t,tau,theta,tresh

allocate(B(1:N),stat=ialloc)
allocate(Z(1:N),stat=ialloc)
A=mat
  do ip=1, N    !initialize V to identity matrix
    do iq=1, N
      V(ip,iq)=0.d0 
    end do
      V(ip,ip)=1.d0
  end do  
  do ip=1, N
    B(ip)=A(ip,ip)
    D(ip)=B(ip)
    Z(ip)=0.d0    
  end do
  NROT=0
  do i=1, 50
    sm=0.d0
    do ip=1, N-1     !sum off-diagonal elements
      do iq=ip+1, N
        sm=sm+DABS(A(ip,iq))
      end do
    end do
    if(sm==0.d0) return  !normal return
    if(i.lt.4) then
      tresh=0.2d0*sm**2
    else
      tresh=0.d0
    end if
    do ip=1, N-1
      do iq=ip+1, N
        g=100.d0*DABS(A(ip,iq))
! after 4 sweeps, skip the rotation if the off-diagonal element is small
        if((i.gt.4).and.(DABS(D(ip))+g.eq.DABS(D(ip))) &
		.and.(DABS(D(iq))+g.eq.DABS(D(iq)))) then
		  A(ip,iq)=0.d0
        else if(DABS(A(ip,iq)).gt.tresh) then
	  h=D(iq)-D(ip)
	  if(DABS(h)+g.eq.DABS(h)) then
	    t=A(ip,iq)/h
          else
	    theta=0.5d0*h/A(ip,iq)  
            t=1.d0/(DABS(theta)+DSQRT(1.d0+theta**2))
	    if(theta.lt.0.d0) t=-t
          end if
	  c=1.d0/DSQRT(1.d0+t**2)
	  s=t*c
          tau=s/(1.d0+c)
	  h=t*A(ip,iq)
	  Z(ip)=Z(ip)-h
	  Z(iq)=Z(iq)+h
	  D(ip)=D(ip)-h
	  D(iq)=D(iq)+h
	  A(ip,iq)=0.d0
	  do j=1, ip-1
	    g=A(j,ip)
	    h=A(j,iq)
	    A(j,ip)=g-s*(h+g*tau)
	    A(j,iq)=h+s*(g-h*tau)
          end do
	  do j=ip+1, iq-1
	    g=A(ip,j)
	    h=A(j,iq)
	    A(ip,j)=g-s*(h+g*tau)
	    A(j,iq)=h+s*(g-h*tau)
          end do		      
	  do j=iq+1, N
	    g=A(ip,j)
	    h=A(iq,j)
	    A(ip,j)=g-s*(h+g*tau)
	    A(iq,j)=h+s*(g-h*tau)
          end do		  
	  do j=1, N
	    g=V(j,ip)
	    h=V(j,iq)
	    V(j,ip)=g-s*(h+g*tau)
	    V(j,iq)=h+s*(g-h*tau)
          end do		  
          NROT=NROT+1
        end if !if ((i.gt.4)...
      end do !main iq loop
    end do !main ip loop
    do ip=1, N
      B(ip)=B(ip)+Z(ip)
      D(ip)=B(ip)
      Z(ip)=0.d0
    end do
  end do !main i loop
  return
END subroutine symeigen   

subroutine powereigen(mat,n,k,testpower,eigenvector,eigenvalue)
use mat_handling
integer n,testpower,k
real mat(n,n)
real*8 mat1(n,n),eigenvector(n,k),eigenvalue(n),eigen(n,1),eigen2(1,n)
real*8 mat2(n,n),eigen1(n,1),dabsmaximum,a,b,eps
eps=0.0001
mat1=mat
eigenvalue=0
eigenvector=0
if(k.gt.n)k=n
do 10 j=1,k
a=1
b=0
testpower=0
eigen1=0
do 1 i=n,1,-1
eigen(i,1)=i
1 continue
do while(a.gt.eps)
testpower=testpower+1
if(testpower.gt.1000) goto 9
call dotproduct8(mat1,eigen,eigen1,n,1,n)
b=dabsmaximum(eigen1,n,1)
if(b.eq.0)goto 9
a=dabs((eigenvalue(j)-b)/b)
eigenvalue(j)=b
eigen=eigen1/b
end do
9 do 8 i=1,n
eigenvector(i,j)=eigen(i,1)
8 continue
a=0
do 7 i=1,n
eigen2(1,i)=eigen(i,1)
a=eigen(i,1)*eigen(i,1)+a
7 continue
if(a.eq.0)goto 10
call dotproduct8(eigen,eigen2,mat2,n,n,1)
mat2=b/a*mat2
mat1=mat1-mat2
10 continue
return
do 20 i=1,k
a=0
do 21 i1=1,n
a=a+eigenvector(i1,i)*eigenvector(i1,i)
21 continue
if(a.gt.0)then
do 22 i1=1,sp
22 eigenvector(i1,i)=eigenvector(i1,i)/a
endif
20 continue
end subroutine powereigen

subroutine dotproduct8 (mat1,mat2,out,sp1,sp2,si)
integer sp1,sp2,si
real*8 mat1(sp1,si),mat2(si,sp2)
real*8 out(sp1,sp2)
out=0
do 1 i=1,sp1
do 2 i1=1,sp2
do 3 j=1,si
out(i,i1)=mat1(i,j)*mat2(j,i1)+out(i,i1)
3 continue
2 continue
1 continue
return
end subroutine dotproduct8

real*8 function dabsmaximum(mat,sp,si)
integer sp,si
real*8 mat(sp,si)
dabsmaximum=-10**9
do 1 i=1,sp
do 2 j=1,si
if (dabs(mat(i,j)).gt.dabsmaximum)dabsmaximum=dabs(mat(i,j))
2 continue
1 continue
end function dabsmaximum

end module

module mult_reg
use matrix_basics
contains
subroutine multreg(x,y,sp,si,b,r2,f,e)
integer sp,si
real r2,f
real*8 x(sp,si),xaug(sp,si+1),xaugt(si+1,sp),xaugtxaug(si+1,si+1),xaugtxauginv(si+1,si+1),temp(si+1,sp),temp1(si,si),y(sp),b(si+1),e(sp)
real*8 rxx(si,si),r(si+1,si+1)
real a,a1
!call fieldwrite(x,sp,si,4)
!call fieldwrite(y,sp,1,4)
do 1 j=1,si+1
do 2 i=1,sp
if(j.eq.1)then
xaug(i,j)=1.0
else
xaug(i,j)=x(i,j-1)
endif
2 continue
1 continue
!call fieldwrite(xaug,sp,si+1,4)
xaugt=transpose(xaug)
!call fieldwrite(xaugt,si+1,sp,4)
call dotproduct_8(xaugt,xaug,xaugtxaug,si+1,si+1,sp)
!call fieldwrite(xaugtxaug,si+1,si+1,4)
call FINDInv_8(xaugtxaug,xaugtxauginv,si+1,ierror)
!call fieldwrite(xaugtxauginv,si+1,si+1,4)
call dotproduct_8(xaugtxauginv,xaugt,temp,si+1,sp,si+1)
!call fieldwrite(temp,si+1,sp,4)
call dotproduct_8(temp,y,b,si+1,1,sp)
!call fieldwrite(b,si+1,1,4)
do 5 i=1,sp
xaug(i,1)=y(i)
5 continue
call zcoltransform_8 (xaug,sp,si+1)
xaugt=transpose(xaug)
call dotproduct_8(xaugt,xaug,r,si+1,si+1,sp)
r=r/sp
!call fieldwrite(r,si+1,si+1,4)
do 6 j=2,si+1
do 7 j1=2,si+1
rxx(j-1,j1-1)=r(j,j1)
7 continue
6 continue
!call fieldwrite(rxx,si,si,4)
a=determinant(r,si+1)
a1=determinant(rxx,si)
if(a1.gt.0)then
r2=1-a/a1
else
r2=0
endif
if(r2.lt.1.and.r2.gt.0.and.sp.gt.2)then
f=sqrt((sp-2)*r2/(1.0-r2))
else
f=0
endif
call dotproduct_8(xaug,b,e,sp,1,si+1)
!call fieldwrite(y,sp,1,4)
!call fieldwrite(e,sp,1,4)
do 8 i=1,sp
e(i)=y(i)-e(i)
8 continue
return
end subroutine multreg

subroutine dotproduct_8(mat1,mat2,out,sp1,sp2,si)
integer sp1,sp2,si
real*8 mat1(sp1,si),mat2(si,sp2)
real*8 out(sp1,sp2)
out=0
do 1 i=1,sp1
do 2 i1=1,sp2
do 3 j=1,si
out(i,i1)=mat1(i,j)*mat2(j,i1)+out(i,i1)
3 continue
2 continue
1 continue
return
end subroutine dotproduct_8 

!Subroutine to find the inverse of a square matrix
!Author : Louisda16th a.k.a Ashwith J. Rego
!Reference : Algorithm has been well explained in:
!http://math.uww.edu/~mcfarlat/inverse.htm           
!http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html
SUBROUTINE FINDInv_8(matrix, inverse, n, errorflag)
	!IMPLICIT NONE
	!Declarations
	INTEGER n,i, i1,j, k, l
	INTEGER errorflag  !Return error status. -1 for error, 0 for normal
	REAL*8 matrix(n,n)  !Input matrix
	REAL*8 inverse(n,n) !Inverted matrix
    REAL*8 m
	REAL*8 augmatrix(n,2*n) !augmented matrix

	LOGICAL:: FLAG = .TRUE.
		
	!Augment input matrix with an identity matrix
	DO i = 1, n
		DO j = 1, 2*n
			IF (j <= n ) THEN
				augmatrix(i,j) = matrix(i,j)
			ELSE IF ((i+n) == j) THEN
				augmatrix(i,j) = 1
			Else
				augmatrix(i,j) = 0
			ENDIF
		END DO
	END DO
	
	!Reduce augmented matrix to upper traingular form
	DO k =1, n-1
		IF (augmatrix(k,k) == 0) THEN
			FLAG = .FALSE.
			DO i = k+1, n
				IF (augmatrix(i,k) /= 0) THEN
					DO j = 1,2*n
						augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j)
					END DO
					FLAG = .TRUE.
					EXIT
				ENDIF
				IF (FLAG .EQV. .FALSE.) THEN
					PRINT*, "Matrix is non - invertible"
					inverse = 0
					errorflag = -1
					return
				ENDIF
			END DO
		ENDIF
		DO j = k+1, n			
			m = augmatrix(j,k)/augmatrix(k,k)
			DO i = k, 2*n
				augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i)
			END DO
		END DO
	END DO
	
	!Test for invertibility
	DO i = 1, n
		IF (augmatrix(i,i) == 0) THEN
			PRINT*, "Matrix is non - invertible"
			inverse = 0
			errorflag = -1
			return
		ENDIF
	END DO
	
	!Make diagonal elements as 1
	DO i = 1 , n
		m = augmatrix(i,i)
		DO j = i , (2 * n)				
			   augmatrix(i,j) = (augmatrix(i,j) / m)
		END DO
	END DO
	
	!Reduced right side half of augmented matrix to identity matrix
	DO k = n-1, 1, -1
		DO i =1, k
		m = augmatrix(i,k+1)
			DO j = k, (2*n)
				augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m
			END DO
		END DO
	END DO				
	
	!store answer
	DO i =1, n
		DO j = 1, n
			inverse(i,j) = augmatrix(i,j+n)
		END DO
	END DO
test=0
do 1 i=1,n
do 2 i1=1,n
do 3 j=1,n
test=matrix(i,j)*inverse(j,i1)+test
3 continue
2 continue
1 continue
errorflag = nint(test)
END SUBROUTINE FINDinv_8

subroutine zcoltransform_8 (mat,sp,si)
integer sp,si
real row(sp),a,b
real*8 mat(sp,si)
do 1 j=1,si
do 2 i=1,sp
row(i)=mat(i,j)
2 continue
a=mean(row,sp)
b=variance(row,sp)
if(b.gt.0)b=sqrt(b)
do 3 i=1,sp
if(b.gt.0)then
mat(i,j)=(mat(i,j)-a)/b
else
mat(i,j)=0
endif
3 continue
1 continue
return
end subroutine zcoltransform_8 

end module

module matrix_basics
use basic_stats
use mat_handling
contains
subroutine dotproduct (mat1,mat2,out,sp1,sp2,si)
integer sp1,sp2,si
real mat1(sp1,si),mat2(si,sp2)
real*8 out(sp1,sp2)
out=0
do 1 i=1,sp1
do 2 i1=1,sp2
do 3 j=1,si
out(i,i1)=mat1(i,j)*mat2(j,i1)+out(i,i1)
3 continue
2 continue
1 continue
return
end subroutine dotproduct 

subroutine dotrealproduct (mat1,mat2,out,sp1,sp2,si)
integer sp1,sp2,si
real mat1(sp1,si),mat2(si,sp2),out(sp1,sp2)
out=0
do 1 i=1,sp1
do 2 i1=1,sp2
do 3 j=1,si
out(i,i1)=mat1(i,j)*mat2(j,i1)+out(i,i1)
3 continue
2 continue
1 continue
return
end subroutine dotrealproduct 

subroutine trendcoltransform (mat,sp,si)
integer sp,si
real row(sp),a
real mat(sp,si)
do 1 j=1,si
do 2 i=1,sp
row(i)=mat(i,j)
2 continue
a=mean(row,sp)
do 3 i=1,sp
mat(i,j)=(mat(i,j)-a)
3 continue
1 continue
return
end subroutine trendcoltransform 

subroutine trendrowtransform (mat,sp,si)
integer sp,si
real col(si),a
real mat(sp,si)
do 1 i=1,sp
do 2 j=1,si
col(j)=mat(i,j)
2 continue
a=mean(col,si)
do 3 j=1,si
mat(i,j)=(mat(i,j)-a)
3 continue
1 continue
return
end subroutine trendrowtransform 

subroutine zcoltransform (mat,sp,si)
integer sp,si
real row(sp),a,b
real mat(sp,si)
do 1 j=1,si
do 2 i=1,sp
row(i)=mat(i,j)
2 continue
a=mean(row,sp)
b=variance(row,sp)
if(b.gt.0)b=sqrt(b)
do 3 i=1,sp
if(b.gt.0)then
mat(i,j)=(mat(i,j)-a)/b
else
mat(i,j)=0
endif
3 continue
1 continue
return
end subroutine zcoltransform 

subroutine zrowtransform (mat,sp,si)
integer sp,si
real col(si),a,b
real mat(sp,si)
do 1 i=1,sp
do 2 j=1,si
col(j)=mat(i,j)
2 continue
a=mean(col,si)
b=variance(col,si)
if(b.gt.0)b=sqrt(b)
do 3 j=1,si
if(b.gt.0)then
mat(i,j)=(mat(i,j)-a)/b
else
mat(i,j)=0
endif
3 continue
1 continue
return
end subroutine zrowtransform 

subroutine proptransform (mat,sp,si)
integer sp,si
real mat(sp,si),a
a=0
do 1 j=1,si
do 2 i=1,sp
a=a+mat(i,j)
2 continue
1 continue
if(a.gt.0)mat=mat/a
return
end subroutine proptransform 

subroutine patransform (mat,sp,si)
integer sp,si
real mat(sp,si)
do 1 i=1,sp
do 2 j=1,si
if(mat(i,j).gt.0)mat(i,j)=1.0
2 continue
1 continue
return
end subroutine patransform 

subroutine indtransform (mat,sp,si)
integer sp,si
real mat(sp,si),a
a=10.0**9
do 1 j=1,si
do 2 i=1,sp
if(mat(i,j).gt.0.and.mat(i,j).lt.a)a=mat(i,j)
2 continue
1 continue
if(a.gt.0.and.a.lt.1)mat=mat/a
return
end subroutine indtransform 

real function maximum(mat,sp,si)
integer sp,si
real mat(sp,si)
maximum=-10**9
do 1 i=1,sp
do 2 j=1,si
if (mat(i,j).gt.maximum)maximum=mat(i,j)
2 continue
1 continue
end function maximum

real function minimum(mat,sp,si)
integer sp,si
real mat(sp,si)
minimum=10**9
do 1 i=1,sp
do 2 j=1,si
if (mat(i,j).lt.minimum)minimum=mat(i,j)
2 continue
1 continue
end function minimum

real function absmaximum(mat,sp,si)
integer sp,si
real mat(sp,si)
absmaximum=-10**9
do 1 i=1,sp
do 2 j=1,si
if (abs(mat(i,j)).gt.absmaximum)absmaximum=abs(mat(i,j))
2 continue
1 continue
end function absmaximum

subroutine distcorr(mat,distmat,n,n1)
integer n,n1
real mat(n,n1),mat1(n,n1),distmat(n1,n1)
mat1=mat
call zcoltransform(mat1,n,n1)
do 1 j=1,n1
do 2 j1=j,n1
a=0
do 3 i=1,n
a=a+mat1(i,j)*mat1(i,j1)
3 continue
a=a/real(n)
distmat(j,j1)=a
distmat(j1,j)=a
2 continue
1 continue
return
!distmat=1.0-distmat
end subroutine distcorr

subroutine disteuclid(mat,distmat,n,n1)
integer n,n1
real mat(n,n1)
real distmat(n,n),mat1(n,n1)
mat1=mat
call zcoltransform(mat1,n,n1)
!call fieldwrite(mat1,n,n1,4)
do 1 i=1,n
do 2 i1=i,n
a=0
do 3 j=1,n1
b=(mat1(i,j)-mat1(i1,j))
a=a+b*b
3 continue
if(a.gt.0)a=sqrt(a)
distmat(i,i1)=a
distmat(i1,i)=a
2 continue
1 continue
!call fieldwrite(distmat,n,n,4)
return
end subroutine disteuclid

subroutine distance(mat,distmat,n,n1)
integer n,n1
real mat(n,n1)
real distmat(n1,n1),mat1(n1,n)
!call fieldwrite(mat,n,n1,4)
mat1=transpose(mat)
call dotrealproduct(mat1,mat,distmat,n1,n1,n)
!call fieldwrite(distmat,n,n,4)
return
end subroutine distance

subroutine distsoerensen (mat1,distmat,sp,si)
integer sp,si
real mat1(sp,si),distmat(si,si)
soerensen=0
do 1 j=1,si
do 2 j1=j,si
a=0
b=0
c=0
do 3 i=1,sp
if(mat1(i,j).gt.0.and.mat1(i,j1).eq.0)a=a+1
if(mat1(i,j).eq.0.and.mat1(i,j1).gt.0)b=b+1
if(mat1(i,j).gt.0.and.mat1(i,j1).gt.0)c=c+1
3 continue
distmat(j,j1)=2.0*c/(a+b+2.0*c)
distmat(j1,j)=distmat(j,j1)
2 continue
1 continue
!call fieldwrite(distmat,si,si,4)
!distmat=1-distmat
!call fieldwrite(distmat,si,si,4)
end subroutine distsoerensen

subroutine covmat(mat,distmat,n,n1)
integer n,n1
real mat(n,n1),mat1(n,n1),distmat(n1,n1)
mat1=mat
call trendcoltransform(mat1,n,n1)
!call fieldwrite(mat1,n,n1,4)
do 1 i=1,n1
do 2 i1=i,n1
b=0
do 3 j=1,n
b=(mat1(j,i)*mat1(j,i1))+b
3 continue
distmat(i,i1)=b/n
distmat(i1,i)=distmat(i,i1)
2 continue
1 continue
return
end subroutine covmat

subroutine chi2distance(mat,sp,si)
integer sp,si
real mat(sp,si),row(sp),col(si)
row=0
col=0
a=0
do 1 i=1,sp
do 2 j=1,si
a=a+mat(i,j)
row(i)=row(i)+mat(i,j)
col(j)=col(j)+mat(i,j)
2 continue
1 continue
mat=mat/a
do 3 i=1,sp
do 4 j=1,si
if(row(i)*col(j).gt.0)then
mat(i,j)=(mat(i,j)*a-row(i)*col(j))/(a*sqrt(row(i)*col(j)))
else
mat(i,j)=0
endif
4 continue
3 continue
return
end subroutine chi2distance

subroutine rowstandard(mat,sp,si)
integer sp,si
real mat(sp,si)
do 1 i=1,sp
a=0
do 2 j=1,si
2 a=mat(i,j)+a
if(a.gt.0)then
do 3 j=1,si
3 mat(i,j)=mat(i,j)/a
endif
1 continue
return
end subroutine rowstandard

subroutine colstandard(mat,sp,si)
integer sp,si
real mat(sp,si)
do 1 j=1,si
a=0
do 2 i=1,sp
2 a=mat(i,j)+a
if(a.gt.0)then
do 3 i=1,sp
3 mat(i,j)=mat(i,j)/a
endif
1 continue
return
end subroutine colstandard

!Subroutine to find the inverse of a square matrix
!Author : Louisda16th a.k.a Ashwith J. Rego
!Reference : Algorithm has been well explained in:
!http://math.uww.edu/~mcfarlat/inverse.htm           
!http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html
SUBROUTINE FINDInv(mat, matout, n, errorflag)
	!IMPLICIT NONE
	!Declarations
	INTEGER n,i, i1,j, k, l
	INTEGER errorflag  !Return error status. -1 for error, 0 for normal
	real mat(n,n),matout(n,n)
	REAL*8 matrix(n,n)  !Input matrix
	REAL*8 inverse(n,n) !Inverted matrix
    REAL*8 m
	REAL*8 augmatrix(n,2*n) !augmented matrix

	LOGICAL:: FLAG = .TRUE.
		
	!Augment input matrix with an identity matrix
	matrix=mat
	DO i = 1, n
		DO j = 1, 2*n
			IF (j <= n ) THEN
				augmatrix(i,j) = matrix(i,j)
			ELSE IF ((i+n) == j) THEN
				augmatrix(i,j) = 1
			Else
				augmatrix(i,j) = 0
			ENDIF
		END DO
	END DO
	
	!Reduce augmented matrix to upper traingular form
	DO k =1, n-1
		IF (augmatrix(k,k) == 0) THEN
			FLAG = .FALSE.
			DO i = k+1, n
				IF (augmatrix(i,k) /= 0) THEN
					DO j = 1,2*n
						augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j)
					END DO
					FLAG = .TRUE.
					EXIT
				ENDIF
				IF (FLAG .EQV. .FALSE.) THEN
					PRINT*, "Matrix is non - invertible"
					inverse = 0
					errorflag = -1
					return
				ENDIF
			END DO
		ENDIF
		DO j = k+1, n			
			m = augmatrix(j,k)/augmatrix(k,k)
			DO i = k, 2*n
				augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i)
			END DO
		END DO
	END DO
	
	!Test for invertibility
	DO i = 1, n
		IF (augmatrix(i,i) == 0) THEN
			PRINT*, "Matrix is non - invertible"
			inverse = 0
			errorflag = -1
			return
		ENDIF
	END DO
	
	!Make diagonal elements as 1
	DO i = 1 , n
		m = augmatrix(i,i)
		DO j = i , (2 * n)				
			   augmatrix(i,j) = (augmatrix(i,j) / m)
		END DO
	END DO
	
	!Reduced right side half of augmented matrix to identity matrix
	DO k = n-1, 1, -1
		DO i =1, k
		m = augmatrix(i,k+1)
			DO j = k, (2*n)
				augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m
			END DO
		END DO
	END DO				
	
	!store answer
	DO i =1, n
		DO j = 1, n
			inverse(i,j) = augmatrix(i,j+n)
		END DO
	END DO
test=0
do 1 i=1,n
do 2 i1=1,n
do 3 j=1,n
test=matrix(i,j)*inverse(j,i1)+test
3 continue
2 continue
1 continue
errorflag = nint(test)
matout=real(inverse)
END SUBROUTINE FINDinv

!Function to find the determinant of a square matrix
!Author : Louisda16th a.k.a Ashwith J. Rego
!Description: The subroutine is based on two key points:
!1] A determinant is unaltered when row operations are performed: Hence, using this principle,
!row operations (column operations would work as well) are used
!to convert the matrix into upper traingular form
!2]The determinant of a triangular matrix is obtained by finding the product of the diagonal elements
!
REAL FUNCTION determinant(matrix, n)
	IMPLICIT NONE
	REAL*8, DIMENSION(n,n) :: matrix
	INTEGER, INTENT(IN) :: n
	REAL :: m, temp
	INTEGER :: i, j, k, l
	LOGICAL :: DetExists = .TRUE.
	l = 1
	!Convert to upper triangular form
	DO k = 1, n-1
		IF (matrix(k,k) == 0) THEN
			DetExists = .FALSE.
			DO i = k+1, n
				IF (matrix(i,k) /= 0) THEN
					DO j = 1, n
						temp = matrix(i,j)
						matrix(i,j)= matrix(k,j)
						matrix(k,j) = temp
					END DO
					DetExists = .TRUE.
					l=-l
					EXIT
				ENDIF
			END DO
			IF (DetExists .EQV. .FALSE.) THEN
				determinant = 0
				return
			END IF
		ENDIF
		DO j = k+1, n
			m = matrix(j,k)/matrix(k,k)
			DO i = k+1, n
				matrix(j,i) = matrix(j,i) - m*matrix(k,i)
			END DO
		END DO
	END DO
		!Calculate determinant by finding product of diagonal elements
	determinant = l
	DO i = 1, n
		determinant = determinant * matrix(i,i)
	END DO
END FUNCTION determinant

end module matrix_basics

module mat_handling
contains
subroutine matreading (mat,spname,siname,sp,si,filen)
integer sp,si,nu,filen
real,allocatable:: mat(:,:)
character*50,allocatable:: spname(:),siname(:)
character*65000 text
character*50 tt(5000)
text=' '
read(filen,'(a)')text
call reading (tt,nu,filen)
si=nu-1
allocate (siname(si))
rewind(filen)
call reading (tt,nu,filen)
do 1 j=2,si+1
read(tt(j),'(a50)',err=2,end=2)siname(j-1)
1 continue
2 sp=0
do while (sp.ge.0)
read(filen,'(a)',err=3,end=3)text
if(text.eq.' ')goto 3
sp=sp+1
enddo
3 allocate (mat(sp,si),spname(sp))
mat=0
rewind(filen)
read(filen,'(a)',err=10,end=10)text
tt=' '
nu=0
do 4 i=1,sp
call reading (tt,nu,filen)
if(tt(1).eq.' ')goto 7
read(tt(1),'(a50)')spname(i)
do 5 j=2,nu
read(tt(j),*)mat(i,j-1)
5 continue
!write(*,*) spname(i)
4 continue
7 return
10 write(*,*)'Error'
end subroutine matreading 

subroutine textreading (mat,spname,siname,sp,si,filen)
integer sp,si,nu,filen
character*50,allocatable:: mat(:,:),spname(:),siname(:)
character*65000 text
character*50 tt(5000)
text=' '
call reading (tt,nu,filen)
si=nu-1
allocate (siname(si))
rewind(filen)
call reading (tt,nu,filen)
do 1 j=2,nu
read(tt(j),'(a50)',err=2,end=2)siname(j-1)
1 continue
2 sp=0
do while (sp.ge.0)
read(filen,'(a)',err=3,end=3)text
if(text.eq.' ')goto 3
sp=sp+1
enddo
3 allocate (mat(sp,si),spname(sp))
mat=' '
rewind(filen)
read(filen,'(a)',err=10,end=10)text
tt=' '
nu=0
!write(*,'(<si>(a5,1x))')(siname(j),j=1,si)
do 4 i=1,sp
call reading (tt,nu,filen)
if(tt(1).eq.' ')goto 7
read(tt(1),'(a50)')spname(i)
do 5 j=2,nu
read(tt(j),'(a50)')mat(i,j-1)
5 continue
4 continue
7 return
10 write(*,*)'Error'
end subroutine textreading 

subroutine matwrite(mat,spname,siname,species,sites,name,filen,file)
integer species,sites,file,filen
real mat(species,sites)
character*50 spname(species),siname(sites),name,text
if(file.gt.0)then
open(14,file='temp.txt',status='replace')
write(14,*)file
backspace (14)
read(14,*)text
backspace (14)
name='TraitMat'
i1=len_trim(name)
name=name(1:i1)//text
endif
i1=len_trim(name)
name=name(1:i1)//'.txt'
open(filen,file=name,status='replace')
write(filen,'(a1,20x,<sites>(a12,1x))')'S',(siname(j),j=1,sites)
!call spaces (filen)
do 8 i=1,species
if(spname(i).ne.' ')write(filen,'(a20,1x,<sites>(f12.3,1x))')spname(i),(mat(i,j),j=1,sites)
!call spaces (filen)
8 continue
close(filen)
close(14)
return
end subroutine matwrite

subroutine filewrite(mat,spname,siname,species,sites,filen)
integer species,sites,filen
real mat(species,sites)
character*50 spname(species),siname(sites)
write(filen,'(a1,20x,<sites>(a14,1x))')'S',(siname(j),j=1,sites)
!call spaces (filen)
do 8 i=1,species
write(filen,'(a20,1x,<sites>(f14.6,1x))')spname(i),(mat(i,j),j=1,sites)
!call spaces (filen)
8 continue
write(filen,*)' '
return
end subroutine filewrite

subroutine fieldwrite(mat,species,sites,filen)
integer species,sites,filen
real mat(species,sites)
character*50 spname(species),siname(sites)
siname='nd'
spname='nd'
write(filen,'(a1,20x,<sites>(a14,1x))')'S',(siname(j),j=1,sites)
!call spaces (filen)
do 8 i=1,species
write(filen,'(a20,1x,<sites>(f14.6,1x))')spname(i),(mat(i,j),j=1,sites)
!call spaces (filen)
8 continue
write(filen,*)' '
return
end subroutine fieldwrite

subroutine reading (tt,nu,filen)
integer nu,length,filen
character*65000 text
character*50 tt(5000)
length=65000
text=' '
tt=' '
read(filen,'(a)',err=4,end=4) text
nu=1
i=1
do while (i.le.length-1)
if(text(i:i).ne.' ')then
do 1 j=i,length
if(text(j:j).eq.' ') then
if(nu.gt.5000) goto 4 
tt(nu)=text(i:j-1)
i=j
nu=nu+1
if(i.gt.length-2.and.nu.lt.2)goto 4
goto 3
endif
1 continue
endif
i=i+1
if(i.gt.length-2.and.nu.lt.2)goto 4
3 enddo
if(nu.lt.2)goto 4
nu=nu-1
4 return
end subroutine reading 

subroutine rowcolelimination(mator,spnameor,sinameor,mat,spname,siname,sp,si)
integer sp,si,row(sp),col(si)
real mator(sp,si)
character*50 spnameor(sp),sinameor(si)
real, allocatable:: mat(:,:),temp(:,:)
character*50, allocatable:: spname(:),siname(:)
row=0
col=0
do 1 i=1,sp
do 2 j=1,si
if(mator(i,j).ne.0)then
row(i)=row(i)+1
col(j)=col(j)+1
endif
2 continue
1 continue 
k=0
do 3 i=1,sp
3 if(row(i).gt.0)k=k+1
k1=0
do 4 j=1,si
4 if(col(j).gt.0)k1=k1+1
allocate (mat(k,k1),temp(k,si),spname(k),siname(k1))
n=0
n1=0
do 5 i=1,sp
if(row(i).gt.0)then
n=n+1
do 6 j=1,si
6 temp(n,j)=mator(i,j)
spname(n)=spnameor(i)
endif
5 continue
n=0
do 7 j=1,si
if(col(j).gt.0)then
n=n+1
do 8 i=1,k
8 mat(i,n)=temp(i,j)
siname(n)=sinameor(j)
endif
7 continue
sp=k
si=k1
return
end subroutine rowcolelimination

end module mat_handling

Module basic_stats
contains

real function mean(a,k)
real m,a(k)
m=0
if(k.gt.0)then
do 1 i=1,k
m=m+a(i)
1 continue
mean=m/k
else
mean=0
endif
end function mean

real function geomean(a,k)
real m,a(k)
m=0
k1=k
geomean=0
if(k.gt.0)then
do 1 i=1,k
if(a(i).gt.0)then
m=m+alog(a(i))
else
k1=k1-1
if(k1.lt.1)return
endif
1 continue
geomean=exp(m/k1)
endif
end function geomean

real function mean0(a,k)
real m,a(k)
m=0
if(k.ge.1)then
j=0
do 1 i=1,k
if(a(i).ne.0)then
m=m+a(i)
j=j+1
endif
1 continue
if(j.gt.0)then
mean0=m/j
else
mean0=0
endif
else
mean0=0
endif
end function mean0

real function variance(a,k)
real mean,square,a(k)
mean=0
square=0
if(k.gt.1)then
do 1 i=1,k
mean=mean+a(i)
square=square+a(i)*a(i)
1 continue
variance=square/k-(mean*mean)/(k*k)
if(variance.lt.0)variance=0
else
variance=0
endif
end function variance

real function standarddeviation(a,k)
real mean,square,a(k)
mean=0
square=0
if(k.gt.1)then
do 1 i=1,k
mean=mean+a(i)
square=square+a(i)*a(i)
1 continue
standarddeviation=square/k-(mean*mean)/(k*k)
if(standarddeviation.le.0)standarddeviation=0
if(standarddeviation.gt.0)standarddeviation=sqrt(standarddeviation)
else
standarddeviation=0
endif
end function standarddeviation

real function variance0(a,k)
real mean,square,a(k)
mean=0
square=0
if(k.gt.1)then
j=0
do 1 i=1,k
if(a(i).ne.0)then
mean=mean+a(i)
square=square+a(i)*a(i)
j=j+1
endif
1 continue
if(j.gt.0)then
variance0=square/j-(mean*mean)/(j*j)
else
variance0=0
endif
else
variance0=0
endif
if(variance0.lt.0)variance0=0
end function variance0

real function standarddeviation0(a,k)
real mean,square,a(k)
mean=0
square=0
if(k.gt.1)then
j=0
do 1 i=1,k
if(a(i).ne.0)then
mean=mean+a(i)
square=square+a(i)*a(i)
j=j+1
endif
1 continue
if(j.gt.0)then
standarddeviation0=square/j-(mean*mean)/(j*j)
if(standarddeviation0.le.0)standarddeviation0=0
if(standarddeviation0.gt.0)standarddeviation0=sqrt(standarddeviation0)
else
standarddeviation0=0
endif
else
standarddeviation0=0
endif
if(standarddeviation0.lt.0)standarddeviation0=0
end function standarddeviation0

real function skewness(a,k)
real b,c,skew,a(k)
!real variance,mean
skew=0
if(k.gt.2)then
b=variance(a,k)
c=mean(a,k)
if(b.gt.0)then
b=sqrt(b)
do 1 i=1,k
skew=skew+(a(i)-c)*(a(i)-c)*(a(i)-c)
1 continue
skewness=k/((k-1.0)*(k-2.0))*skew/(b*b*b)
if(skewness.gt.999.0)skewness=999.0
else
skewness=0
endif
else
skewness=0
endif
end function skewness

real function skewness0(a,k)
real b,c,skew,a(k)
!real mean0,variance0
skew=0
k1=0
if(k.gt.2)then
b=variance0(a,k)
c=mean0(a,k)
if(b.gt.0)then
b=sqrt(b)
do 1 i=1,k
if(a(i).ne.0)then
skew=skew+(a(i)-c)*(a(i)-c)*(a(i)-c)
k1=k1+1
endif
1 continue
if(k1.gt.2)then
skewness0=k1/((k1-1.0)*(k1-2.0))*skew/(b*b*b)
else
skewness0=0
endif
if(skewness0.gt.999.0)skewness0=999.0
else
skewness0=0
endif
else
skewness0=0
endif
end function skewness0

subroutine conflim (field,i,up,down,a,p)
integer i,i1,i2
real field(i),up,down,a,p
i1=nint(0.025*i)
if(i1.lt.1)i1=1
i2=nint(0.975*i)
if(i2.gt.i)i2=i
call sort_vector_down(field,i)
up=field(i1)
down=field(i2)
p=0.5
do 1 i1=1,i
if(i1.lt.0.5*i)then
if(field(i1).lt.a)then
p=max(real(i1)/i,1.0/i)
return
endif
else
if(field(i1).le.a)then
p=max(1.0-real(i1)/i,1.0/i)
return
endif
endif
1 continue
if(a.lt.field(i))p=1.0/i
return
end subroutine conflim

subroutine conflim0 (field,i,up,down,a,p)
integer i,i1,i2
real field(i),up,down,a,p
call sort_vector_down(field,i)
do 2 i3=1,i
2 if(i3.gt.0)goto 3
3 i3=i3-1
i1=nint(0.025*(i-i3))+i3
if(i1.lt.1)i1=1
i2=nint(0.975*(i-i3))+i3
if(i2.gt.i)i2=i
up=field(i1)
down=field(i2)
p=0.5
do 1 i1=1,i
if(i1.lt.0.5*i)then
if(field(i1).lt.a)then
p=max(real(i1)/i,1.0/i)
return
endif
else
if(field(i1).le.a)then
p=max(1.0-real(i1)/i,1.0/i)
return
endif
endif
1 continue
if(a.lt.field(i))p=1.0/i
return
end subroutine conflim0

real function error(a)
real a,b,c,pi
pi=3.14159
c=a*0.70711
b=0.140012*c*c
error=c*c*(4.0/pi+b)/(1.0+b)
error=(1.0-exp(-error))**0.5
end function error

real function spearman(a,b,k)
real a(k),b(k)
!real pearson
call rank (a,k,1)
call rank (b,k,1)
spearman=pearson(a,b,k)
end function spearman

real function pearson(a,b,k)
real a(k),b(k),a1,a2,b1,b2,c
!real mean,variance
pearson=0
a1=mean(a,k)
a2=mean(b,k)
b1=variance(a,k)
b2=variance(b,k)
if(b1.gt.0.and.b2.gt.0)then
b1=sqrt(b1)
b2=sqrt(b2)
do 1 i=1,k
pearson=pearson+(a(i)-a1)*(b(i)-a2)
1 continue
pearson=pearson/(k*b1*b2)
else
pearson=0
endif
end function pearson


real function pearson0(a,b,k)
real a(k),b(k),a1,a2,b1,b2,c
!real mean0,variance0
pearson0=0
a1=mean0(a,k)
a2=mean0(b,k)
b1=variance0(a,k)
b2=variance0(b,k)
if(b1.gt.0.and.b2.gt.0)then
b1=sqrt(b1)
b2=sqrt(b2)
do 1 i=1,k
if(a(i).ne.0.and.b(i).ne.0)pearson0=pearson0+(a(i)-a1)*(b(i)-a2)
1 continue
pearson0=pearson0/(k*b1*b2)
else
pearson0=0
endif
end function pearson0

subroutine rank (out,sp,si)
integer sp,si
real out(sp,si),row(sp)
do 10 j=1,si
do 1 i=1,sp
row(i)=out(i,j)
1 continue
k=1
do while(k.le.sp)
a=10000000
do 2 i=1,sp
if(row(i).lt.a)a=row(i)
2 continue
i1=0
do 3 i=1,sp
if(row(i).eq.a)then
out(i,j)=k
row(i)=100000000
i1=i1+1
endif
3 continue
do 4 i2=1,sp
if(out(i2,j).eq.k)out(i2,j)=k+0.5*(i1-1)
4 continue
k=k+i1
enddo
10 continue
return
end subroutine rank 

subroutine seriation (field,rowscore,colscore,sp,si,iseed)
integer sp,si,iseed
real field(sp,si),rowscore(sp),colscore(si),rowscoreor(sp),colscoreor(si),a,b,sum,m,s
!real mean,variance
do 10 i=1,sp
rowscore(i)=ran(iseed)
10 continue
colscore=0
rowscoreor=0
colscoreor=0
a=1
b=1
k=1
do while (a+b.gt.0.0001)
a=0
b=0
do 2 j=1,si
m=0
sum=0
do 3 i=1,sp
m=m+field(i,j)*rowscore(i)
sum=sum+field(i,j)
!write(4,*)rowscore(i),m
3 continue
if(sum.ne.0)colscore(j)=m/sum
2 continue
m=mean(colscore,si)
s=variance(colscore,si)
do 1 j=1,si
if(s.gt.0)then
colscore(j)=(colscore(j)-m)/sqrt(s)
else
colscore(j)=0
endif
a=a+abs(colscore(j)-colscoreor(j))
1 continue
rowscore=0
do 4 i=1,sp
m=0
sum=0
do 5 j=1,si
m=m+field(i,j)*colscore(j)
sum=sum+field(i,j)
5 continue
if(sum.ne.0)rowscore(i)=m/sum
4 continue
m=mean(rowscore,sp)
s=variance(rowscore,sp)
do 6 i=1,sp
if(s.gt.0)then
rowscore(i)=(rowscore(i)-m)/sqrt(s)
else
rowscore(i)=0
endif
b=b+abs(rowscore(i)-rowscoreor(i))
6 continue
rowscoreor=rowscore
colscoreor=colscore
k=k+1
if(k.gt.1000)then
write(*,*)'Ordination did not converge after 1000 steps.' 
return
endif
enddo
return
end subroutine seriation

subroutine sort_vector_down(field,i)
integer i
real field(i),an
do 2 j1=1,i-1
do 1 i1=j1+1,i
an=field(j1)
if(field(i1).gt.an)then
field(j1)=field(i1)
field(i1)=an
endif
1 continue
2 continue
return
end subroutine sort_vector_down

integer function counter0(a,k)
real a(k)
counter0=0
do 1 i=1,k
1 if(a(i).ne.0)counter0=counter0+1
end function counter0

integer function countitem(b,k)
character*50 b(k)
countitem=0
do 1 i=1,k-1
do 2 i1=i+1,k
2 if(b(i1).eq.b(i))b(i1)=' '
1 continue
do 3 i=1,k
3 if(b(i).ne.' ')countitem=countitem+1
end function countitem

real function chi2fit(a,b,n)
integer n
real a(n),b(n)
chi2fit=0
do 1 i=1,n
1 chi2fit=((a(i)-b(i))*(a(i)-b(i)))+chi2fit
end function chi2fit

real function relerror(a,b,m,n,n1)
integer m,n,n1
real a(m,n),b(m,n1)
relerror=0
k=0
do 1 j=1,n
do 2 j1=1,n1
do 3 i=1,m
k=k+1
3 if(a(i,j).ne.0.and.b(i,j1).ne.0)relerror=abs(a(i,j)-b(i,j1))/abs(a(i,j)+b(i,j1))+relerror
2 continue
1 continue
if(k.gt.0)relerror=relerror/k
end function relerror

subroutine probdist (in,sp,ipos,iseed)
integer ipos,sp
real an,in(sp),probf(sp)
probf=0
probf(1)=in(1)
do 1 i=2,sp
probf(i)=in(i)+probf(i-1)
1 continue
an=ran(iseed)*probf(sp)
do 2 ipos=1,sp
if(an.lt.probf(ipos))return
2 continue
ipos=0
return
end subroutine probdist

subroutine mantel(mat,mat1,sp,si,rmat)
integer sp,si
real mat(sp,si),mat1(sp,si),a(sp*si),b(sp*si),rmat
k=0
do 1 i=1,sp
do 2 j=1,si
k=k+1
a(k)=mat(i,j)
b(k)=mat1(i,j)
2 continue
1 continue
!rmat=pearson(a,b,sp*si)
rmat=spearman(a,b,sp*si)
end subroutine mantel

subroutine manteloffdiag(mat,mat1,sp,rmat)
integer sp,si
real mat(sp,sp),mat1(sp,sp),a(sp*sp-sp),b(sp*sp-sp),rmat
k=0
do 1 i=1,sp
do 2 j=1,sp
if(i.ne.j)then
k=k+1
a(k)=mat(i,j)
b(k)=mat1(i,j)
endif
2 continue
1 continue
!rmat=pearson(a,b,sp*sp-sp)
rmat=spearman(a,b,sp*sp-sp)
end subroutine manteloffdiag

end module


module sorting
contains

subroutine sortrow(field,rowscore,spname,sp,si,si1,k,dir)
integer sp,si,si1
real field(sp,si),rowscore(sp,si1),an(si),an1(si1)
character*10 dir
character*50,spname(sp),b
do 2 j1=1,sp-1
do 1 i1=j1+1,sp
if(dir.eq.'up')then
if(rowscore(i1,k).lt.rowscore(j1,k))then
do 5 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
5 continue
do 6 k1=1,si1
an1(k1)=rowscore(j1,k1)
rowscore(j1,k1)=rowscore(i1,k1)
rowscore(i1,k1)=an1(k1)
6 continue
b=spname(j1)
spname(j1)=spname(i1)
spname(i1)=b
endif
elseif(dir.eq.'down')then
if(rowscore(i1,k).gt.rowscore(j1,k))then
do 15 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
15 continue
do 16 k1=1,si1
an1(k1)=rowscore(j1,k1)
rowscore(j1,k1)=rowscore(i1,k1)
rowscore(i1,k1)=an1(k1)
16 continue
b=spname(j1)
spname(j1)=spname(i1)
spname(i1)=b
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sortrow

subroutine sortcol(field,colscore,siname,sp,si,si1,k,dir)
integer sp,si,si1
real field(sp,si),colscore(si,si1),an(sp),an1(si1)
character*10 dir
character*50,siname(si),b
do 2 j1=1,si-1
do 1 i1=j1+1,si
if(dir.eq.'up')then
if(colscore(i1,k).lt.colscore(j1,k))then
do 5 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
5 continue
do 6 k1=1,si1
an1(k1)=colscore(j1,k1)
colscore(j1,k1)=colscore(i1,k1)
colscore(i1,k1)=an1(k1)
6 continue
b=siname(j1)
siname(j1)=siname(i1)
siname(i1)=b
endif
elseif(dir.eq.'down')then
if(colscore(i1,k).gt.colscore(j1,k))then
do 15 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
15 continue
do 16 k1=1,si1
an1(k1)=colscore(j1,k1)
colscore(j1,k1)=colscore(i1,k1)
colscore(i1,k1)=an1(k1)
16 continue
b=siname(j1)
siname(j1)=siname(i1)
siname(i1)=b
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sortcol

subroutine sortmatrow(field,rowscore,sp,si,si1,k,dir)
integer sp,si,si1
real field(sp,si),rowscore(sp,si1),an(si),an1(si1)
character*10 dir
do 2 j1=1,sp-1
do 1 i1=j1+1,sp
if(dir.eq.'up')then
if(rowscore(i1,k).lt.rowscore(j1,k))then
do 5 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
5 continue
do 6 k1=1,si1
an1(k1)=rowscore(j1,k1)
rowscore(j1,k1)=rowscore(i1,k1)
rowscore(i1,k1)=an1(k1)
6 continue
endif
elseif(dir.eq.'down')then
if(rowscore(i1,k).gt.rowscore(j1,k))then
do 15 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
15 continue
do 16 k1=1,si1
an1(k1)=rowscore(j1,k1)
rowscore(j1,k1)=rowscore(i1,k1)
rowscore(i1,k1)=an1(k1)
16 continue
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sortmatrow

subroutine sortmatcol(field,colscore,sp,si,si1,k,dir)
integer sp,si,si1
real field(sp,si),colscore(si,si1),an(sp),an1(si1)
character*10 dir
do 2 j1=1,si-1
do 1 i1=j1+1,si
if(dir.eq.'up')then
if(colscore(i1,k).lt.colscore(j1,k))then
do 5 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
5 continue
do 6 k1=1,si1
an1(k1)=colscore(j1,k1)
colscore(j1,k1)=colscore(i1,k1)
colscore(i1,k1)=an1(k1)
6 continue
endif
elseif(dir.eq.'down')then
if(colscore(i1,k).gt.colscore(j1,k))then
do 15 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
15 continue
do 16 k1=1,si1
an1(k1)=colscore(j1,k1)
colscore(j1,k1)=colscore(i1,k1)
colscore(i1,k1)=an1(k1)
16 continue
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sortmatcol

subroutine sort_vector(field,i,dir)
integer i
real field(i),an
character*10 dir
do 2 j1=1,i-1
do 1 i1=j1+1,i
an=field(j1)
if(dir.eq.'up')then
if(field(i1).lt.an)then
field(j1)=field(i1)
field(i1)=an
endif
elseif(dir.eq.'down')then
if(field(i1).gt.an)then
field(j1)=field(i1)
field(i1)=an
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sort_vector

subroutine sort_textvector(field,i,dir)
integer i
character*50 field(i),an
character*10 dir
do 2 j1=1,i-1
do 1 i1=j1+1,i
an=field(j1)
if(dir.eq.'up')then
if(field(i1).lt.an)then
field(j1)=field(i1)
field(i1)=an
endif
elseif(dir.eq.'down')then
if(field(i1).gt.an)then
field(j1)=field(i1)
field(i1)=an
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sort_textvector

subroutine sortrowname(field,spname,sp,si,dir)
integer sp,si
real field(sp,si),an(si)
character*10 dir
character*50,spname(sp),b
do 2 j1=1,sp-1
do 1 i1=j1+1,sp
if(dir.eq.'up')then
if(spname(i1).lt.spname(j1))then
do 5 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
5 continue
b=spname(j1)
spname(j1)=spname(i1)
spname(i1)=b
endif
elseif(dir.eq.'down')then
if(spname(i1).gt.spname(j1))then
do 15 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
15 continue
b=spname(j1)
spname(j1)=spname(i1)
spname(i1)=b
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sortrowname

subroutine sortcolname(field,siname,sp,si,dir)
integer sp,si
real field(sp,si),an(sp)
character*10 dir
character*50,siname(si),b
do 2 j1=1,si-1
do 1 i1=j1+1,si
if(dir.eq.'up')then
if(siname(i1).lt.siname(j1))then
do 5 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
5 continue
b=siname(j1)
siname(j1)=siname(i1)
siname(i1)=b
endif
elseif(dir.eq.'down')then
if(siname(i1).gt.siname(j1))then
do 15 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
15 continue
b=siname(j1)
siname(j1)=siname(i1)
siname(i1)=b
endif
else
write(*,*)'No sorting direction'
return
endif
1 continue
2 continue
return
end subroutine sortcolname

subroutine sortmatspabund(field,row,sp,si,si1)
integer sp,si,si1
real row(sp,si1),a,c
real field(sp,si),an(si)
do 2 j1=1,sp-1
do 1 i1=j1+1,sp
if(row(i1,1).eq.row(j1,1).and.row(i1,2).gt.row(j1,2))then
do 5 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
5 continue
a=row(j1,1)
c=row(j1,2)
row(j1,1)=row(i1,1)
row(j1,2)=row(i1,2)
row(i1,1)=a
row(i1,2)=c
endif
1 continue
2 continue
return
end subroutine sortmatspabund

subroutine sortmatsiabund(field,col,sp,si,si1)
integer sp,si,si1
real field(sp,si),an(sp)
real col(si,si1),a,c
do 2 j1=1,si-1
do 1 i1=j1+1,si
if(col(i1,1).eq.col(j1,1).and.col(i1,2).gt.col(j1,2))then
do 5 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
5 continue
a=col(j1,1)
c=col(j1,2)
col(j1,1)=col(i1,1)
col(j1,2)=col(i1,2)
col(i1,1)=a
col(i1,2)=c
endif
1 continue
2 continue
return
end subroutine sortmatsiabund

subroutine sortspabund(field,row,spname,sp,si,si1)
integer sp,si,si1
real row(sp,si1),a,c
real field(sp,si),an(si)
character*50,spname(sp),b
do 2 j1=1,sp-1
do 1 i1=j1+1,sp
if(row(i1,1).eq.row(j1,1).and.row(i1,2).gt.row(j1,2))then
do 5 k1=1,si
an(k1)=field(j1,k1)
field(j1,k1)=field(i1,k1)
field(i1,k1)=an(k1)
5 continue
a=row(j1,1)
c=row(j1,2)
b=spname(j1)
row(j1,1)=row(i1,1)
row(j1,2)=row(i1,2)
spname(j1)=spname(i1)
row(i1,1)=a
row(i1,2)=c
spname(i1)=b
endif
1 continue
2 continue
return
end subroutine sortspabund

subroutine sortsiabund(field,col,siname,sp,si,si1)
integer sp,si,si1
real field(sp,si),an(sp)
real col(si,si1),a,c
character*50,siname(si),b
do 2 j1=1,si-1
do 1 i1=j1+1,si
if(col(i1,1).eq.col(j1,1).and.col(i1,2).gt.col(j1,2))then
do 5 k1=1,sp
an(k1)=field(k1,j1)
field(k1,j1)=field(k1,i1)
field(k1,i1)=an(k1)
5 continue
a=col(j1,1)
c=col(j1,2)
b=siname(j1)
col(j1,1)=col(i1,1)
col(j1,2)=col(i1,2)
siname(j1)=siname(i1)
col(i1,1)=a
col(i1,2)=c
siname(i1)=b
endif
1 continue
2 continue
return
end subroutine sortsiabund

end module sorting