ccc TO USE, TO USE: Analysis of the network of meetings between 
ccc individuals (strong links)  + size of largest cluster.
ccc If monkey j y l meet at time t, afinity(j,l,t)=1
ccc connect(j,l) = sum over t of affinity(j,l,t)
ccc Comparison with average degree of weak+strong links network.
      program preferencial
      implicit none
      common/dim/lx
      common/var/s(300)
      common/bond/statusbond(300,300),nbneigh(300),
     & ineigh(300,300)           
      real*8 x(300,300),y(300,300)
      real*8 avaff(300),dx,dy,err,varaff,varaffj,avaffjl
      real*8 degavrun,vecavrun,clustavrun,sizemaxavrun,clusternoavrun
      real*8 sizemaxstavrun,clusternostavrun
      real*8 rndavaffjl,avdeg,degj,avvec,clust,pdeg(0:300)
      real*4 rkiss,r
      integer ntime,nmonos,lx,t0,i,j,id(300),k(300,300),meet(300,300)
      integer affinity(300,300,300),l,n,prod,na,nd,nbrun,nrun,runok
      integer connect(300,300),rndconnect(300,300),nb,nbnosolos
      integer newconnect(300,300),newdegree(300),deg,degmax,ntot
      integer alone(300),nvecinos(300),degree(300),nlink,kc,kb
      integer rndaffinity(300,300,300),nbmeet(300),ndum(300),s,sizemax
      integer rndnbmeet,kisset
      integer clusterno,clusternost,sizemaxst
      integer statusbond,nbneigh,ineigh
      character*40 inputdat,outputdat,namerun(0:100)
#     include "namerun.f"
cccccccccccccc PARAMETERS: (t0: transient time)
      ntime=100
      nmonos=100
      t0=1
      nbrun=50
      outputdat='affinity.network.perco'
ccccccccccccccc
      i=kisset(1,1,1)
      r=rkiss()
      lx=nmonos
      runok=0
      degavrun=0.
      vecavrun=0.
      clustavrun=0.
      sizemaxavrun=0.
      sizemaxstavrun=0.
      clusternoavrun=0.
      clusternostavrun=0.
      degmax=0
      ntot=0
      do 101 deg=0,nmonos
       pdeg(deg)=0.
101   continue
      do 100 nrun=0,nbrun-1
      inputdat=namerun(nrun)
      open(unit=1,file=inputdat,
     &    status='old',form='formatted')
      do 1 i=1,ntime
       read(1,*) (id(j),x(j,i),y(j,i),k(j,i), j=1,nmonos)
1     continue
      close(1)
cccccccccccccccccccc
      open(unit=2,file=outputdat,
     &    status='new',form='formatted')
cccc determine affinity matrix (1 if coincide at time t), and connection 
ccc  (integrated)
      do 2 j=1,nmonos
       alone(j)=1
       nvecinos(j)=0
       degree(j)=0
       do 12 l=1,nmonos
        i=0
        meet(j,l)=0
        connect(j,l)=0
        do 13 n=1,ntime
         affinity(j,l,n)=0
         dx=abs(x(j,n)-x(l,n))
         dy=abs(y(j,n)-y(l,n))
         err=0.0000001
         if(dx.le.err.AND.dy.le.err) then
          affinity(j,l,n)=1
          meet(j,l)=1
         endif
         if(j.eq.l) then 
          affinity(j,l,n)=0
          meet(j,l)=0
         endif
         connect(j,l)=connect(j,l)+affinity(j,l,n)
13      continue 
        nvecinos(j)=nvecinos(j)+meet(j,l)
        degree(j)=degree(j)+connect(j,l)       
12     continue
       if(nvecinos(j).ge.1) alone(j)=0
2     continue
cc compute link number and average degree:
      nlink=0
      avdeg=0.
      do 6 j=1,nmonos
       avdeg=avdeg+real(degree(j))
       nlink=nlink+degree(j)
6     continue
      avdeg=avdeg/real(nmonos)
      nlink=nlink/2  
ccccccccccccccccc
cc compute perco. prop. (weak+strong links) of meet(i,j) -> largest cluster
      do 66 j=1,nmonos
       s(j)=1
       nbneigh(j)=0
       do 67 l=1,nmonos
        if(l.eq.j) goto 67
        if(meet(j,l).eq.1) then
         nbneigh(j)=nbneigh(j)+1
         ineigh(j,nbneigh(j))=l
         statusbond(j,nbneigh(j))=1
        endif
67     continue
66    continue
      call clustercounting(sizemax,clusterno)
cc      write(*,*) sizemax, clusterno
ccccccccccccccccc
cc make a new connection matrix based on p{ER}(k_c)=0.05. If 
cc connect(j,l) > kc then newconnect(j,l)=1, 0 otherwise.
      do 23 j=1,nmonos
       call kcrit(degree(j),nvecinos(j),kc)
       do 24 l=1,nmonos
        if(connect(j,l).ge.kc) then
         newconnect(j,l)=1
        else
         newconnect(j,l)=0
        endif
24     continue
23    continue
cccccccccccc
cc compute perco. prop. of newconnect (strong links) -> largest cluster
      do 166 j=1,nmonos
       s(j)=1
       nbneigh(j)=0
       do 167 l=1,nmonos
        if(l.eq.j) goto 167
        if(newconnect(j,l).eq.1) then
         nbneigh(j)=nbneigh(j)+1
         ineigh(j,nbneigh(j))=l
         statusbond(j,nbneigh(j))=1
        endif
167    continue
166   continue
      call clustercounting(sizemaxst,clusternost)
cc degre distribution of the new network: 
      do 25 j=1,nmonos
       newdegree(j)=0
       do 26 l=1,nmonos
        if(newconnect(j,l).eq.1) newdegree(j)=newdegree(j)+1
26     continue
25    continue
      do 27 j=1,nmonos
       pdeg(newdegree(j))=pdeg(newdegree(j))+1.
       ntot=ntot+1
       if(newdegree(j).gt.degmax) degmax=newdegree(j)
27    continue     
cccccccccccccccccc
cc compute new average degree, compared with original neighbours, + clustering:
      avdeg=0.
      avvec=0.
      do 7 j=1,nmonos
       degj=0.
       do 8 l=1,nmonos 
        degj=degj+real(newconnect(j,l))
8      continue 
       avdeg=avdeg+degj
       avvec=avvec+nvecinos(j)
7     continue
      avdeg=avdeg/real(nmonos)
      avvec=avvec/real(nmonos)
      call clustering(newconnect,nmonos,clust,nb)
cc      write(*,*) 'mean degree=',avdeg,'  original nb. neighbors=', avvec
cc      write(*,*) 'clustering=', clust, ' nb node  deg>1 =', nb
ccccc general averages:
      degavrun=degavrun+avdeg
      vecavrun=vecavrun+avvec
      sizemaxavrun=sizemaxavrun+sizemax
      sizemaxstavrun=sizemaxstavrun+sizemaxst
      clusternoavrun=clusternoavrun+clusterno
      clusternostavrun=clusternostavrun+clusternost
      if(clust.le.1.01) then
       runok=runok+1
       clustavrun=clustavrun+clust
      endif
100   continue
      degavrun=degavrun/real(nbrun)
      vecavrun=vecavrun/real(nbrun)
      sizemaxavrun=sizemaxavrun/real(nbrun)
      sizemaxstavrun=sizemaxstavrun/real(nbrun)
      clusternoavrun=clusternoavrun/real(nbrun)
      clusternostavrun=clusternostavrun/real(nbrun)
      clustavrun=clustavrun/real(runok)
      write(2,*) 'mean degree of strong links (over runs)=', degavrun,
     &  '   mean degree of strong+weak links (over runs)', vecavrun
      write(2,*) 'mean clustering coef. (strong links)  
     &  over runs=', clustavrun
      write(2,*)'mean size of largest cluster in network(weak+strong)=',
     & sizemaxavrun
      write(2,*)'mean number of clusters in network(weak+strong)=',
     & clusternoavrun
      write(2,*)'mean size of largest cluster in network(strong)=',
     & sizemaxstavrun
      write(2,*)'mean number of clusters in network(strong)=',
     & clusternostavrun
      write(2,*) 'degree distribution of strong links [k,p(k)]:'
      do 102 deg=0,degmax
       write(2,*) deg, pdeg(deg)/real(ntot)
102   continue 
      close(2)        
      end
 
      include "subrout.clust.net.f"
      
      function rkiss()
	real*8 rm
c	data rm / 4294967295.0d0 /
	data rm / 2147483648.0d0 /
	rkiss=dfloat(abs(kiss())+1)/rm
      return
      end



      function kiss()
c      data i,j,k/123456789,362436069,521288629/
c      data i,j,k/607504,6057433,609760129/
      common/rnd/i,j,k
      m(k,n)=ieor(k,ishft(k,n))
        i=69069*i+23606797
        j=m(m(j,17),-15)
        k=m(iand(m(k,18),2147483647),-13)
        kiss=i+j+k
      return
      end
      
      function kisset(ii,jj,kk)
      common/rnd/i,j,k
c	print*,' Current seeds:',i,j,k
      i=ii
      j=jj
      k=kk
      kisset=1
      return
      end
      
      subroutine kcrit(degree,vecinos,kc)
cc compute prob that a node with given degree and # of neighbors
cc has k connections to a given neighbor, if all connections are random.
cc p(kc)=0.05 as a definition.
      implicit none
      real*8 p(0:40000),npair,prod
      integer k,degree,i,vecinos,kc
      kc=1
      npair=real(vecinos)
      if(degree.eq.0) then
       kc=1
       goto 3
      endif
      p(0)=(1.-1./npair)**degree
      do 1 k=1,degree
       prod=1
       do 2 i=0,k-1
        prod=prod*real((degree-i))/real(i+1)
2      continue
       p(k)=prod
     &    *((1./npair)**k)*((1.-1./npair)**(degree-k))
       if(p(k-1).ge.0.05.AND.p(k).le.0.05) then
         kc=k
         goto 3
       endif
1     continue
3     continue
      end
       
      subroutine clustering(connect,nmonos,clust,nb)
      implicit none
      integer connect(300,300),nvecinos(300),neigh(300,300)
      integer nmonos,j,l,nb,cj,n1,n2,l1,l2
      real*8 clust
      do 1 j=1,nmonos
       nvecinos(j)=0
       do 2 l=1,nmonos
cc        if(connect(j,l).ge.1.OR.connect(l,j).ge.1) then
        if(connect(j,l).ge.1) then
         nvecinos(j)=nvecinos(j)+1
         neigh(j,nvecinos(j))=l
        endif
2      continue
1     continue
      nb=0
      clust=0.
      do 3 j=1,nmonos
       if(nvecinos(j).lt.2) goto 3
       nb=nb+1
       cj=0
       do 4 n1=1,nvecinos(j)
        l1=neigh(j,n1)
        do 5 n2=1,nvecinos(j)
         if(n1.eq.n2) goto 5
         l2=neigh(j,n2)
         if(connect(l1,l2).ge.1) cj=cj+1
5       continue
4      continue
       cj=cj
       clust=clust
     &   +real(cj)/real(nvecinos(j)*(nvecinos(j)-1))
3     continue
      if(nb.eq.0) then
       clust=2
      else
       clust=clust/real(nb)
      endif
      end        


       
 
