ccc 3nd version: TO USE, TO USE. Compute affinity, EXCLUDING case where 
ccc monkey met only one other monkey during the whole run. 
ccc Affinity matrix lines: TIME.
ccc Comparison with random series.
      program preferencial
      implicit none
      real*8 x(300,300),y(300,300)
      real*8 avaff(300),dx,dy,err,varaff,varaffj,avaffjl
      real*8 varaffrun,rndvaraffrun,rndvaraff,rndvaraffj
      real*8 rndavaffjl
      real*4 rkiss,r
      integer ntime,nmonos,t0,i,j,id(300),k(300,300),meet(300,300)
      integer affinity(300,300,300),l,n,prod,na,nd,nbrun,nrun
      integer narb(300),duration(300,300),nb,nbnosolos
      integer alone(300),nvecinos(300),runok
      integer rndaffinity(300,300,300),nbmeet(300),ndum(300)
      integer rndnbmeet,kisset
      character*40 inputdat,outputdat,namerun(0:100)
#     include "namerun.f"
ccccccccc PARAMETERS: (t0: transient time)
      ntime=50
      nmonos=200
      t0=1
      nbrun=50
      outputdat='affinity.time.rnd'
ccccccccccccc
      i=kisset(1,1,1)
      r=rkiss()
      varaffrun=0.
      rndvaraffrun=0.
      runok=0
      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)
      do 2 j=1,nmonos
       alone(j)=1
       nvecinos(j)=0
       do 12 l=1,nmonos
        i=0
        meet(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
13      continue 
        nvecinos(j)=nvecinos(j)+meet(j,l)       
12     continue
       if(nvecinos(j).ge.2) alone(j)=0
cc compute average affinity of j only among the l's it met (nvecinos.ge.2):
       avaff(j)=0.
       nbmeet(j)=0
       if(alone(j).eq.1) goto 16
       nb=0
       do 14 l=1,nmonos
        if(meet(j,l).eq.0) goto 14
        do 15 n=1,ntime
         nbmeet(j)=nbmeet(j)+affinity(j,l,n)
         avaff(j)=avaff(j)+real(affinity(j,l,n))
         nb=nb+1
15      continue
14     continue   
       if(nb.gt.0) avaff(j)=avaff(j)/real(nb)
16     continue
2     continue
ccccccccccccccccccc
cccc create a random affinity matrix with same avaff(j) and 1's in first raw:
      do 20 j=1,nmonos
       if(nvecinos(j).lt.2) goto 20
       do 21 l=1,nvecinos(j)
        n=1
        rndaffinity(j,l,n)=1
        ndum(l)=2
        do 22 n=2,ntime       
         rndaffinity(j,l,n)=0
22      continue
21     continue
       rndnbmeet=nvecinos(j)
       if(rndnbmeet.eq.nbmeet(j)) goto 24
       do 23 n=1,10*nbmeet(j)
        r=rkiss()
        l=int(r*real(nvecinos(j)))+1   
        if(ndum(l).le.ntime) then     
         rndaffinity(j,l,ndum(l))=1
         ndum(l)=ndum(l)+1
         rndnbmeet=rndnbmeet+1
         if(rndnbmeet.eq.nbmeet(j)) goto 24
        endif
23     continue
24     continue
20    continue
cccccccccccccccccc
ccccc compute average affinities and variance:
      varaff=0.
      rndvaraff=0.
      nbnosolos=0
      do 3 j=1,nmonos
       if(alone(j).eq.1) goto 30
       nbnosolos=nbnosolos+1
       varaffj=0.
       nb=0
       do 4 l=1,nmonos
        if(meet(j,l).eq.0) goto 40
        avaffjl=0.
        do 5 n=1,ntime
         avaffjl=avaffjl+real(affinity(j,l,n))
5       continue
        avaffjl=avaffjl/real(ntime)
        varaffj=varaffj+(avaffjl-avaff(j))**2/(avaff(j))**2
        nb=nb+1
40      continue
4      continue
       varaffj=(varaffj/real(nb))**.5
cc       write(2,*) 'mono #',j, ':  variance in affinity=',varaffj
       varaff=varaff+varaffj
cc compute variance of rndaffinity
       rndvaraffj=0.
       do 41 l=1,nvecinos(j)
        rndavaffjl=0.
        do 42 n=1,ntime
         rndavaffjl=rndavaffjl+real(rndaffinity(j,l,n))
42      continue
        rndavaffjl=rndavaffjl/real(ntime)
        rndvaraffj=rndvaraffj
     &       +(rndavaffjl-avaff(j))**2/(avaff(j))**2
41     continue
       rndvaraffj=(rndvaraffj/real(nvecinos(j)))**.5
       rndvaraff=rndvaraff+rndvaraffj
cccccccccccccccccccccccccccc 
30     continue
3     continue
      if(nbnosolos.gt.0) then
       varaff=varaff/real(nbnosolos)
       rndvaraff=rndvaraff/real(nbnosolos)
       runok=runok+1
       write(2,*) 'variance in affinity=', varaff, 
     & '; random=', rndvaraff
      endif
      varaffrun=varaffrun+varaff
      rndvaraffrun=rndvaraffrun+rndvaraff
100   continue
      varaffrun=varaffrun/real(runok)
      rndvaraffrun=rndvaraffrun/real(runok)
      write(2,*) 'variance averaged over runs=', varaffrun,
     &  '; random=', rndvaraffrun
      close(2)        
      end


      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
      

