--- /dev/null
+c===========================================================================
+c
+c This file is part of TISEAN
+c
+c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+c
+c TISEAN is free software; you can redistribute it and/or modify
+c it under the terms of the GNU General Public License as published by
+c the Free Software Foundation; either version 2 of the License, or
+c (at your option) any later version.
+c
+c TISEAN is distributed in the hope that it will be useful,
+c but WITHOUT ANY WARRANTY; without even the implied warranty of
+c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+c GNU General Public License for more details.
+c
+c You should have received a copy of the GNU General Public License
+c along with TISEAN; if not, write to the Free Software
+c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+c
+c===========================================================================
+c clustering a dissimilarity matrix
+c see Schreiber and Schmitz, Phys. Rev. Lett. 79 (1997) 1475
+c author T. Schreiber (1998)
+c===========================================================================
+
+ parameter(npmax=1000)
+ dimension d(npmax,npmax), iu(npmax), ifix(npmax)
+ character*72 file, fout, filex
+ data iverb/3/
+
+ call whatido("clustering a dissimilarity matrix",iverb)
+ ncl=imust("#")
+ iflag=lopt("=",1)
+ call stcan('X',filex,' ')
+ isout=igetout(fout,iverb)
+
+ call nthstring(1,file)
+ call infile(file,iunit,iverb)
+ if(file.eq."-") file="stdin"
+ if(isout.eq.1) call addsuff(fout,file,"_clust")
+ do 10 i=1,npmax
+ do 10 j=1,npmax
+ 10 d(i,j)=-1e20
+ np=0
+ 1 read(iunit,*,end=999) i,j,dij
+ d(i,j)=dij
+ np=max(i,j,np)
+ goto 1
+ 999 if(iv_io(iverb).eq.1) write(0,'(a,i)') "matrix size ", np
+ dmean=0
+ nd=0
+ do 20 i=1,np
+ do 20 j=1,np
+ if(d(i,j).ne.-1e20) then
+ nd=nd+1
+ dmean=dmean+d(i,j)
+ endif
+ 20 continue
+ do 30 i=1,np
+ do 30 j=1,np
+ 30 if(d(i,j).eq.-1e20) d(i,j)=dmean/nd
+ do 40 i=1,np
+ 40 ifix(i)=0
+ if(filex.ne." ") then
+ open(10,file=filex,status='old',err=998)
+ nfix=0
+ 2 read(10,*,end=998,err=2) i, iff
+ if(i.lt.1.or.i.gt.np.or.iff.gt.ncl.or.iff.lt.1) goto 1
+ ifix(i)=iff
+ nfix=nfix+1
+ endif
+ 998 if(nfix.eq.np) stop "all fixed."
+ call clustering(np,d,npmax,ncl,nfix,ifix,iu,iverb,iflag)
+ call outfile(fout,iunit,iverb)
+ do 50 n=1,np
+ 50 write(iunit,*) iu(n), (costi(np,iu,d,n,ic,iflag),ic=1,ncl)
+ end
+
+ subroutine usage()
+c usage message
+
+ call whatineed("-## [-= -X xfile] file")
+ call popt("#","number of clusters")
+ call popt("=","if set, bias towards similar size clusters")
+ call popt("X","list of indices with fixed cluster assignments")
+ call pout("file_clust")
+ call pall()
+ call ptext("Verbosity levels (add what you want):")
+ call ptext(" 1 = input/output" )
+ call ptext(" 2 = state of clustering")
+ call ptext(" 8 = temperature / cost at cooling")
+ stop
+ end
+
+ subroutine clustering(np,d,npmax,ncl,nfix,ifix,iu,iverb,iflag)
+ parameter(nt0=20,tfac=10.,tstep=0.99,ntotmaxf=20,nsuccmaxf=2)
+ external rand
+ character*1 c
+ dimension d(npmax,npmax), iu(*), ifix(*)
+ equivalence (c,ic)
+ data c/'A'/
+
+ ntotmax=(np-nfix)*ntotmaxf
+ nsuccmax=(np-nfix)*nsuccmaxf
+ se=0.
+ se2=0.
+ do 10 nt=1,nt0
+ call ranconf(np,iu,ncl,ifix)
+ e=cost(np,iu,d,ncl,iflag)
+ se=se+e
+ 10 se2=se2+e**2
+ t=tfac*sqrt(se2/nt0-(se/nt0)**2)
+
+ ntot=0
+ nsucc=0
+ 1 call cconf(np,iu,ncl,nch,iuold,ifix)
+ ec=cost(np,iu,d,ncl,iflag)
+ ntot=ntot+1
+ if(ec.lt.e.or.(rand(0.0).lt.exp(-(ec-e)/t))) then
+ e=ec
+ nsucc=nsucc+1
+ else
+ iu(nch)=iuold
+ endif
+ if(ntot.eq.ntotmax .or. nsucc.eq.nsuccmax) then
+ if(nsucc.eq.0) return
+ ntot=0
+ nsucc=0
+ if(iv_clust(iverb).eq.1) write(istderr(),'(80a1)')
+ . (ic+iu(n)-1,n=1,np)
+ if(iv_cool(iverb).eq.1) write(istderr(),*) t, e
+ t=t*tstep
+ endif
+ goto 1
+ end
+
+ function cost(np,iu,d,ncl,iflag)
+ parameter(npmax=1000)
+ dimension d(npmax,npmax), iu(*), ictab(npmax)
+
+ cost=0
+ do 10 ic=1,ncl
+ nic=0
+ do 20 n=1,np
+ if(iu(n).ne.ic) goto 20
+ nic=nic+1
+ ictab(nic)=n
+ 20 continue
+ cc=0
+ do 30 ii=1,nic
+ i=ictab(ii)
+ do 30 jj=1,nic
+ j=ictab(jj)
+ 30 cc=cc+d(i,j)
+ 10 if(nic.gt.0) cost=cost+cc/(1+(1-iflag)*(nic-1))
+ end
+
+ function costi(np,iu,d,nn,ic,iflag)
+ parameter(npmax=1000)
+ dimension d(npmax,npmax), iu(*), ictab(npmax)
+
+ costi=0
+ nic=0
+ do 20 n=1,np
+ if(iu(n).ne.ic) goto 20
+ nic=nic+1
+ ictab(nic)=n
+ 20 continue
+ cc=0
+ do 30 jj=1,nic
+ j=ictab(jj)
+ 30 cc=cc+d(nn,j)+d(j,nn)
+ if(nic.gt.0) costi=0.5*cc/(1+(1-iflag)*(nic-1))
+ end
+
+ subroutine ranconf(np,iu,ncl,ifix)
+ external rand
+ dimension iu(*), ifix(*)
+
+ do 10 n=1,np
+ iu(n)=ifix(n)
+ 10 if(ifix(n).eq.0) iu(n)=min(int(rand(0.0)*ncl)+1,ncl)
+ end
+
+ subroutine cconf(np,iu,ncl,nch,iuold,ifix)
+ external rand
+ dimension iu(*), ifix(*)
+
+ 1 nch=min(int(rand(0.0)*np)+1,np)
+ if(ifix(nch).ne.0) goto 1
+ iuold=iu(nch)
+ iu(nch)=iuold+int(rand(0.0)*(ncl-1))+1
+ if(iu(nch).gt.ncl) iu(nch)=iu(nch)-ncl
+ end
git://source.jalview.org / jabaws.git / blobdiff