binaries/src/disembl/Tisean_3.0.1/source_f/cluster.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c   clustering a dissimilarity matrix
  23 c   see Schreiber and Schmitz, Phys. Rev. Lett. 79 (1997) 1475
  24 c   author T. Schreiber (1998)
  25 c===========================================================================
  26
  27       parameter(npmax=1000)
  28       dimension d(npmax,npmax), iu(npmax), ifix(npmax)
  29       character*72 file, fout, filex
  30       data iverb/3/
  31
  32       call whatido("clustering a dissimilarity matrix",iverb)
  33       ncl=imust("#")
  34       iflag=lopt("=",1)
  35       call stcan('X',filex,' ')
  36       isout=igetout(fout,iverb)
  37
  38       call nthstring(1,file)
  39       call infile(file,iunit,iverb)
  40       if(file.eq."-") file="stdin"
  41       if(isout.eq.1) call addsuff(fout,file,"_clust")
  42       do 10 i=1,npmax
  43          do 10 j=1,npmax
  44  10         d(i,j)=-1e20
  45       np=0
  46  1    read(iunit,*,end=999) i,j,dij
  47       d(i,j)=dij
  48       np=max(i,j,np)
  49       goto 1
  50  999  if(iv_io(iverb).eq.1) write(0,'(a,i)') "matrix size ", np
  51       dmean=0
  52       nd=0
  53       do 20 i=1,np
  54          do 20 j=1,np
  55             if(d(i,j).ne.-1e20) then
  56                nd=nd+1
  57                dmean=dmean+d(i,j)
  58             endif
  59  20         continue
  60       do 30 i=1,np
  61          do 30 j=1,np
  62  30         if(d(i,j).eq.-1e20) d(i,j)=dmean/nd
  63       do 40 i=1,np
  64  40      ifix(i)=0
  65       if(filex.ne." ") then
  66          open(10,file=filex,status='old',err=998)
  67          nfix=0
  68  2       read(10,*,end=998,err=2) i, iff
  69          if(i.lt.1.or.i.gt.np.or.iff.gt.ncl.or.iff.lt.1) goto 1
  70          ifix(i)=iff
  71          nfix=nfix+1
  72       endif
  73  998  if(nfix.eq.np) stop "all fixed."
  74       call clustering(np,d,npmax,ncl,nfix,ifix,iu,iverb,iflag)
  75       call outfile(fout,iunit,iverb)
  76       do 50 n=1,np
  77  50      write(iunit,*) iu(n), (costi(np,iu,d,n,ic,iflag),ic=1,ncl)
  78       end
  79
  80       subroutine usage()
  81 c usage message
  82
  83       call whatineed("-## [-= -X xfile] file")
  84       call popt("#","number of clusters")
  85       call popt("=","if set, bias towards similar size clusters")
  86       call popt("X","list of indices with fixed cluster assignments")
  87       call pout("file_clust")
  88       call pall()
  89       call ptext("Verbosity levels (add what you want):")
  90       call ptext("          1 = input/output" )
  91       call ptext("          2 = state of clustering")
  92       call ptext("          8 = temperature / cost at cooling")
  93       stop
  94       end
  95
  96       subroutine clustering(np,d,npmax,ncl,nfix,ifix,iu,iverb,iflag)
  97       parameter(nt0=20,tfac=10.,tstep=0.99,ntotmaxf=20,nsuccmaxf=2)
  98       external rand
  99       character*1 c
 100       dimension d(npmax,npmax), iu(*), ifix(*)
 101       equivalence (c,ic)
 102       data c/'A'/
 103
 104       ntotmax=(np-nfix)*ntotmaxf
 105       nsuccmax=(np-nfix)*nsuccmaxf
 106       se=0.
 107       se2=0.
 108       do 10 nt=1,nt0
 109          call ranconf(np,iu,ncl,ifix)
 110          e=cost(np,iu,d,ncl,iflag)
 111          se=se+e
 112  10      se2=se2+e**2
 113       t=tfac*sqrt(se2/nt0-(se/nt0)**2)
 114
 115       ntot=0
 116       nsucc=0
 117  1    call cconf(np,iu,ncl,nch,iuold,ifix)
 118       ec=cost(np,iu,d,ncl,iflag)
 119       ntot=ntot+1
 120       if(ec.lt.e.or.(rand(0.0).lt.exp(-(ec-e)/t))) then
 121          e=ec
 122          nsucc=nsucc+1
 123       else
 124          iu(nch)=iuold
 125       endif
 126       if(ntot.eq.ntotmax .or. nsucc.eq.nsuccmax) then
 127          if(nsucc.eq.0) return
 128          ntot=0
 129          nsucc=0
 130          if(iv_clust(iverb).eq.1) write(istderr(),'(80a1)')
 131      .      (ic+iu(n)-1,n=1,np)
 132          if(iv_cool(iverb).eq.1) write(istderr(),*) t, e
 133          t=t*tstep
 134       endif
 135       goto 1
 136       end
 137
 138       function cost(np,iu,d,ncl,iflag)
 139       parameter(npmax=1000)
 140       dimension d(npmax,npmax), iu(*), ictab(npmax)
 141
 142       cost=0
 143       do 10 ic=1,ncl
 144          nic=0
 145          do 20 n=1,np
 146             if(iu(n).ne.ic) goto 20
 147             nic=nic+1
 148             ictab(nic)=n
 149  20         continue
 150          cc=0
 151          do 30 ii=1,nic
 152             i=ictab(ii)
 153             do 30 jj=1,nic
 154                j=ictab(jj)
 155  30            cc=cc+d(i,j)
 156  10      if(nic.gt.0) cost=cost+cc/(1+(1-iflag)*(nic-1))
 157       end
 158
 159       function costi(np,iu,d,nn,ic,iflag)
 160       parameter(npmax=1000)
 161       dimension d(npmax,npmax), iu(*), ictab(npmax)
 162
 163       costi=0
 164       nic=0
 165       do 20 n=1,np
 166          if(iu(n).ne.ic) goto 20
 167          nic=nic+1
 168          ictab(nic)=n
 169  20      continue
 170       cc=0
 171       do 30 jj=1,nic
 172          j=ictab(jj)
 173  30      cc=cc+d(nn,j)+d(j,nn)
 174       if(nic.gt.0) costi=0.5*cc/(1+(1-iflag)*(nic-1))
 175       end
 176
 177       subroutine ranconf(np,iu,ncl,ifix)
 178       external rand
 179       dimension iu(*), ifix(*)
 180
 181       do 10 n=1,np
 182          iu(n)=ifix(n)
 183  10      if(ifix(n).eq.0) iu(n)=min(int(rand(0.0)*ncl)+1,ncl)
 184       end
 185
 186       subroutine cconf(np,iu,ncl,nch,iuold,ifix)
 187       external rand
 188       dimension iu(*), ifix(*)
 189
 190  1    nch=min(int(rand(0.0)*np)+1,np)
 191       if(ifix(nch).ne.0) goto 1
 192       iuold=iu(nch)
 193       iu(nch)=iuold+int(rand(0.0)*(ncl-1))+1
 194       if(iu(nch).gt.ncl) iu(nch)=iu(nch)-ncl
 195       end