1 c===========================================================================
3 c This file is part of TISEAN
5 c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
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.
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.
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
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===========================================================================
28 dimension d(npmax,npmax), iu(npmax), ifix(npmax)
29 character*72 file, fout, filex
32 call whatido("clustering a dissimilarity matrix",iverb)
35 call stcan('X',filex,' ')
36 isout=igetout(fout,iverb)
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")
46 1 read(iunit,*,end=999) i,j,dij
50 999 if(iv_io(iverb).eq.1) write(0,'(a,i)') "matrix size ", np
55 if(d(i,j).ne.-1e20) then
62 30 if(d(i,j).eq.-1e20) d(i,j)=dmean/nd
66 open(10,file=filex,status='old',err=998)
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
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)
77 50 write(iunit,*) iu(n), (costi(np,iu,d,n,ic,iflag),ic=1,ncl)
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")
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")
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)
100 dimension d(npmax,npmax), iu(*), ifix(*)
104 ntotmax=(np-nfix)*ntotmaxf
105 nsuccmax=(np-nfix)*nsuccmaxf
109 call ranconf(np,iu,ncl,ifix)
110 e=cost(np,iu,d,ncl,iflag)
113 t=tfac*sqrt(se2/nt0-(se/nt0)**2)
117 1 call cconf(np,iu,ncl,nch,iuold,ifix)
118 ec=cost(np,iu,d,ncl,iflag)
120 if(ec.lt.e.or.(rand(0.0).lt.exp(-(ec-e)/t))) then
126 if(ntot.eq.ntotmax .or. nsucc.eq.nsuccmax) then
127 if(nsucc.eq.0) return
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
138 function cost(np,iu,d,ncl,iflag)
139 parameter(npmax=1000)
140 dimension d(npmax,npmax), iu(*), ictab(npmax)
146 if(iu(n).ne.ic) goto 20
156 10 if(nic.gt.0) cost=cost+cc/(1+(1-iflag)*(nic-1))
159 function costi(np,iu,d,nn,ic,iflag)
160 parameter(npmax=1000)
161 dimension d(npmax,npmax), iu(*), ictab(npmax)
166 if(iu(n).ne.ic) goto 20
173 30 cc=cc+d(nn,j)+d(j,nn)
174 if(nic.gt.0) costi=0.5*cc/(1+(1-iflag)*(nic-1))
177 subroutine ranconf(np,iu,ncl,ifix)
179 dimension iu(*), ifix(*)
183 10 if(ifix(n).eq.0) iu(n)=min(int(rand(0.0)*ncl)+1,ncl)
186 subroutine cconf(np,iu,ncl,nch,iuold,ifix)
188 dimension iu(*), ifix(*)
190 1 nch=min(int(rand(0.0)*np)+1,np)
191 if(ifix(nch).ne.0) goto 1
193 iu(nch)=iuold+int(rand(0.0)*(ncl-1))+1
194 if(iu(nch).gt.ncl) iu(nch)=iu(nch)-ncl