a1a089e430498d663bc48cfeebdc5eb4c7d1fd6c
[unres.git] / source / cluster / wham / src-M / hc.f
1 C***********************  Contents  ****************************************
2 C* Sample driver program, VAX-11 Fortran; **********************************
3 C* HC: O(n^2) time, O(n^2) space hierarchical clustering, Fortran 77 *******
4 C* HCASS: determine cluster-memberships, Fortran 77. *********************** 
5 C* HCDEN: draw upper part of dendrogram, VAX-11 Fortran. *******************
6 C* Sample data set: last 36 lines. *****************************************
7 C***************************************************************************
8 C      REAL DATA(18,16),CRIT(18),MEMBR(18)
9 C      REAL CRITVAL(9)
10 C      INTEGER IA(18),IB(18)
11 C      INTEGER ICLASS(18,9),HVALS(9)
12 C      INTEGER IORDER(9),HEIGHT(9)
13 C      DIMENSION NN(18),DISNN(18)
14 C      REAL D(153)
15 C      LOGICAL FLAG(18)
16 C IN ABOVE, 18=N, 16=M, 9=LEV, 153=N(N-1)/2.
17 C
18 C
19 C       OPEN(UNIT=21,STATUS='OLD',FILE='SPECTR.DAT')
20 C
21 C
22 C      N = 18
23 C      M = 16
24 C      DO I=1,N
25 C        READ(21,100)(DATA(I,J),J=1,M)        
26 C      ENDDO
27 C 100  FORMAT(8F7.1)
28 C
29 C
30 C      LEN = (N*(N-1))/2
31 C      IOPT=1
32 C      CALL HC(N,M,LEN,IOPT,DATA,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,D)
33 C
34 C
35 C      LEV = 9
36 C      CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
37 C
38 C
39 C      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
40 C
41 C
42 C      END
43 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
44 C                                                            C
45 C  HIERARCHICAL CLUSTERING using (user-specified) criterion. C
46 C                                                            C
47 C  Parameters:                                               C
48 C                                                            C
49 Cremoved  DATA(N,M)         input data matrix,               C
50 C  DISS(LEN)         dissimilarities in lower half diagonal  C
51 C                    storage; LEN = N.N-1/2,                 C
52 C  IOPT              clustering criterion to be used,        C
53 C  IA, IB, CRIT      history of agglomerations; dimensions   C
54 C                    N, first N-1 locations only used,       C
55 C  MEMBR, NN, DISNN  vectors of length N, used to store      C 
56 C                    cluster cardinalities, current nearest  C
57 C                    neighbour, and the dissimilarity assoc. C
58 C                    with the latter.                        C
59 C  FLAG              boolean indicator of agglomerable obj./ C
60 C                    clusters.                               C
61 C                                                            C
62 C  F. Murtagh, ESA/ESO/STECF, Garching, February 1986.       C
63 C                                                            C
64 C------------------------------------------------------------C
65       SUBROUTINE HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,
66      X                FLAG,DISS)
67       REAL MEMBR(N)
68       REAL DISS(LEN)
69       INTEGER IA(N),IB(N)
70       REAL CRIT(N)
71       DIMENSION NN(N),DISNN(N)
72       LOGICAL FLAG(N)
73       REAL INF
74       DATA INF/1.E+20/
75 C
76 C  Initializations
77 C
78       DO I=1,N
79          MEMBR(I)=1.
80          FLAG(I)=.TRUE.
81       ENDDO
82       NCL=N
83 C
84 C  Construct dissimilarity matrix
85 C
86       DO I=1,N-1
87          DO J=I+1,N
88             IND=IOFFSET(N,I,J)
89 cinput            DISS(IND)=0.
90 cinput            DO K=1,M
91 cinput               DISS(IND)=DISS(IND)+(DATA(I,K)-DATA(J,K))**2
92 cinput            ENDDO
93             IF (IOPT.EQ.1) DISS(IND)=DISS(IND)/2.
94 C           (Above is done for the case of the min. var. method
95 C            where merging criteria are defined in terms of variances
96 C            rather than distances.)
97           ENDDO
98        ENDDO
99 C
100 C  Carry out an agglomeration - first create list of NNs
101 C
102       DO I=1,N-1
103          DMIN=INF
104          DO J=I+1,N
105             IND=IOFFSET(N,I,J)
106             IF (DISS(IND).GE.DMIN) GOTO 500
107                DMIN=DISS(IND)
108                JM=J
109   500    CONTINUE
110          ENDDO
111          NN(I)=JM
112          DISNN(I)=DMIN
113       ENDDO
114 C
115   400 CONTINUE
116 C     Next, determine least diss. using list of NNs
117       DMIN=INF
118       DO I=1,N-1
119          IF (.NOT.FLAG(I)) GOTO 600
120          IF (DISNN(I).GE.DMIN) GOTO 600
121             DMIN=DISNN(I)
122             IM=I
123             JM=NN(I)
124   600    CONTINUE
125       ENDDO
126       NCL=NCL-1
127 C
128 C  This allows an agglomeration to be carried out.
129 C
130       I2=MIN0(IM,JM)
131       J2=MAX0(IM,JM)
132       IA(N-NCL)=I2
133       IB(N-NCL)=J2
134       CRIT(N-NCL)=DMIN
135 C
136 C  Update dissimilarities from new cluster.
137 C
138       FLAG(J2)=.FALSE.
139       DMIN=INF
140       DO K=1,N
141          IF (.NOT.FLAG(K)) GOTO 800
142          IF (K.EQ.I2) GOTO 800
143          X=MEMBR(I2)+MEMBR(J2)+MEMBR(K)
144          IF (I2.LT.K) THEN
145                            IND1=IOFFSET(N,I2,K)
146                       ELSE
147                            IND1=IOFFSET(N,K,I2)
148          ENDIF
149          IF (J2.LT.K) THEN
150                            IND2=IOFFSET(N,J2,K)
151                       ELSE
152                            IND2=IOFFSET(N,K,J2)
153          ENDIF
154          IND3=IOFFSET(N,I2,J2)
155          XX=DISS(IND3)
156 C
157 C  WARD'S MINIMUM VARIANCE METHOD - IOPT=1.
158 C
159          IF (IOPT.EQ.1) THEN
160             DISS(IND1)=(MEMBR(I2)+MEMBR(K))*DISS(IND1)+
161      X                 (MEMBR(J2)+MEMBR(K))*DISS(IND2)-
162      X                 MEMBR(K)*XX
163             DISS(IND1)=DISS(IND1)/X
164          ENDIF
165 C
166 C  SINGLE LINK METHOD - IOPT=2.
167 C
168          IF (IOPT.EQ.2) THEN
169             DISS(IND1)=MIN(DISS(IND1),DISS(IND2))
170          ENDIF
171 C
172 C  COMPLETE LINK METHOD - IOPT=3.
173 C
174          IF (IOPT.EQ.3) THEN
175             DISS(IND1)=MAX(DISS(IND1),DISS(IND2))
176          ENDIF
177 C
178 C  AVERAGE LINK (OR GROUP AVERAGE) METHOD - IOPT=4.
179 C
180          IF (IOPT.EQ.4) THEN
181             DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2))/
182      X                 (MEMBR(I2)+MEMBR(J2))
183          ENDIF
184 C
185 C  MCQUITTY'S METHOD - IOPT=5.
186 C
187          IF (IOPT.EQ.5) THEN
188             DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2)
189          ENDIF
190 C
191 C  MEDIAN (GOWER'S) METHOD - IOPT=6.
192 C
193          IF (IOPT.EQ.6) THEN
194             DISS(IND1)=0.5*DISS(IND1)+0.5*DISS(IND2)-0.25*XX
195          ENDIF
196 C
197 C  CENTROID METHOD - IOPT=7.
198 C
199          IF (IOPT.EQ.7) THEN
200             DISS(IND1)=(MEMBR(I2)*DISS(IND1)+MEMBR(J2)*DISS(IND2)-
201      X          MEMBR(I2)*MEMBR(J2)*XX/(MEMBR(I2)+MEMBR(J2)))/
202      X          (MEMBR(I2)+MEMBR(J2))
203             ENDIF
204 C
205          IF (I2.GT.K) GOTO 800
206          IF (DISS(IND1).GE.DMIN) GOTO 800
207             DMIN=DISS(IND1)
208             JJ=K
209   800    CONTINUE
210       ENDDO
211       MEMBR(I2)=MEMBR(I2)+MEMBR(J2)
212       DISNN(I2)=DMIN
213       NN(I2)=JJ
214 C
215 C  Update list of NNs insofar as this is required.
216 C
217       DO I=1,N-1
218          IF (.NOT.FLAG(I)) GOTO 900
219          IF (NN(I).EQ.I2) GOTO 850
220          IF (NN(I).EQ.J2) GOTO 850
221          GOTO 900
222   850    CONTINUE
223 C        (Redetermine NN of I:)
224          DMIN=INF
225          DO J=I+1,N
226             IND=IOFFSET(N,I,J)
227             IF (.NOT.FLAG(J)) GOTO 870
228             IF (I.EQ.J) GOTO 870
229             IF (DISS(IND).GE.DMIN) GOTO 870
230                DMIN=DISS(IND)
231                JJ=J
232   870       CONTINUE
233          ENDDO
234          NN(I)=JJ
235          DISNN(I)=DMIN
236   900    CONTINUE
237       ENDDO
238 C
239 C  Repeat previous steps until N-1 agglomerations carried out.
240 C
241       IF (NCL.GT.1) GOTO 400
242 C
243 C
244       RETURN
245       END
246 C
247 C
248       FUNCTION IOFFSET(N,I,J)
249 C  Map row I and column J of upper half diagonal symmetric matrix 
250 C  onto vector.
251       IOFFSET=J+(I-1)*N-(I*(I+1))/2
252       RETURN
253       END
254 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
255 C                                                               C
256 C  Given a HIERARCHIC CLUSTERING, described as a sequence of    C
257 C  agglomerations, derive the assignments into clusters for the C
258 C  top LEV-1 levels of the hierarchy.                           C
259 C  Prepare also the required data for representing the          C
260 C  dendrogram of this top part of the hierarchy.                C
261 C                                                               C
262 C  Parameters:                                                  C
263 C                                                               C
264 C  IA, IB, CRIT: vectors of dimension N defining the agglomer-  C
265 C                 ations.                                       C
266 C  LEV:          number of clusters in largest partition.       C
267 C  HVALS:        vector of dim. LEV, used internally only.      C
268 C  ICLASS:       array of cluster assignments; dim. N by LEV.   C
269 C  IORDER, CRITVAL, HEIGHT: vectors describing the dendrogram,  C
270 C                all of dim. LEV.                               C
271 C                                                               C
272 C  F. Murtagh, ESA/ESO/STECF, Garching, February 1986.          C
273 C                                                               C
274 C HISTORY                                                       C
275 C                                                               C
276 C Bounds bug fix, Oct. 1990, F. Murtagh.                        C
277 C Inserted line "IF (LOC.GT.LEV) GOTO 58" on line 48.  This was C
278 C occassioned by incorrect termination of this loop when I      C
279 C reached its (lower) extremity, i.e. N-LEV.  Without the       C
280 C /CHECK=(BOUNDS) option on VAX/VMS compilation, this inserted  C
281 C statement was not necessary.                                  C
282 C---------------------------------------------------------------C
283       SUBROUTINE HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,
284      X        CRITVAL,HEIGHT)
285       include 'sizesclu.dat'
286       include 'COMMON.IOUNITS'
287       integer ICLASS(maxconf,maxconf-1)
288       INTEGER IA(N),IB(N),HVALS(LEV),IORDER(LEV),
289      X        HEIGHT(LEV)
290       REAL CRIT(N),CRITVAL(LEV)
291 C
292 C  Pick out the clusters which the N objects belong to,
293 C  at levels N-2, N-3, ... N-LEV+1 of the hierarchy.
294 C  The clusters are identified by the lowest seq. no. of
295 C  their members.
296 C  There are 2, 3, ... LEV clusters, respectively, for the
297 C  above levels of the hierarchy.
298 C
299       HVALS(1)=1
300       HVALS(2)=IB(N-1)
301       LOC=3
302       DO 59 I=N-2,N-LEV,-1
303          DO 52 J=1,LOC-1
304             IF (IA(I).EQ.HVALS(J)) GOTO 54
305   52     CONTINUE
306          HVALS(LOC)=IA(I)
307          LOC=LOC+1
308   54     CONTINUE
309          DO 56 J=1,LOC-1
310             IF (IB(I).EQ.HVALS(J)) GOTO 58
311   56     CONTINUE
312          IF (LOC.GT.LEV) GOTO 58
313          HVALS(LOC)=IB(I)
314          LOC=LOC+1
315   58     CONTINUE
316   59  CONTINUE
317 C
318       DO 400 LEVEL=N-LEV,N-2
319          DO 200 I=1,N
320             ICL=I
321             DO 100 ILEV=1,LEVEL
322   100       IF (IB(ILEV).EQ.ICL) ICL=IA(ILEV)
323             NCL=N-LEVEL
324             ICLASS(I,NCL-1)=ICL
325   200    CONTINUE
326   400  CONTINUE
327 C
328       DO 120 I=1,N
329       DO 120 J=1,LEV-1
330       DO 110 K=2,LEV
331       IF (ICLASS(I,J).NE.HVALS(K)) GOTO 110
332          ICLASS(I,J)=K
333          GOTO 120
334   110 CONTINUE
335   120 CONTINUE
336 C
337       WRITE (iout,450) (j,j=2,LEV)
338   450 FORMAT(4X,' SEQ NOS',8(i2,'CL'),10000(i3,'CL'))
339       WRITE (iout,470) (' ---',j=2,LEV)
340   470 FORMAT(4X,' -------',10000a4)
341       DO 500 I=1,N
342       WRITE (iout,600) I,(ICLASS(I,J),J=1,LEV-1) 
343   600 FORMAT(I11,8I4,10000i5)                    
344   500 CONTINUE
345 C
346 C  Determine an ordering of the LEV clusters (at level LEV-1)
347 C  for later representation of the dendrogram.
348 C  These are stored in IORDER.
349 C  Determine the associated ordering of the criterion values
350 C  for the vertical lines in the dendrogram.
351 C  The ordinal values of these criterion values may be used in
352 C  preference, and these are stored in HEIGHT.
353 C  Finally, note that the LEV clusters are renamed so that they
354 C  have seq. nos. 1 to LEV.
355 C
356       IORDER(1)=IA(N-1)
357       IORDER(2)=IB(N-1)
358       CRITVAL(1)=0.0
359       CRITVAL(2)=CRIT(N-1)
360       HEIGHT(1)=LEV
361       HEIGHT(2)=LEV-1
362       LOC=2
363       DO 700 I=N-2,N-LEV+1,-1
364          DO 650 J=1,LOC
365             IF (IA(I).EQ.IORDER(J)) THEN
366 C              Shift rightwards and insert IB(I) beside IORDER(J):
367                DO 630 K=LOC+1,J+1,-1
368                   IORDER(K)=IORDER(K-1)
369                   CRITVAL(K)=CRITVAL(K-1)
370                   HEIGHT(K)=HEIGHT(K-1)
371   630          CONTINUE
372                IORDER(J+1)=IB(I)
373                 CRITVAL(J+1)=CRIT(I)
374                 HEIGHT(J+1)=I-(N-LEV)
375                LOC=LOC+1
376             ENDIF
377   650   CONTINUE
378   700 CONTINUE
379       DO 705 I=1,LEV
380          DO 703 J=1,LEV
381             IF (HVALS(I).EQ.IORDER(J)) THEN
382                IORDER(J)=I
383                GOTO 705
384             ENDIF
385   703    CONTINUE
386   705 CONTINUE
387 C
388       RETURN
389       END
390 C+++++++++++++++++++++++++++++++++++++++++++++++++C
391 C                                                 C
392 C  Construct a DENDROGRAM of the top 8 levels of  C
393 C  a HIERARCHIC CLUSTERING.                       C
394 C                                                 C
395 C  Parameters:                                    C
396 C                                                 C
397 C  IORDER, HEIGHT, CRITVAL: vectors of length LEV C
398 C          defining the dendrogram.               C
399 C          These are: the ordering of objects     C
400 C          along the bottom of the dendrogram     C
401 C          (IORDER); the height of the vertical   C
402 C          above each object, in ordinal values   C
403 C          (HEIGHT); and in real values (CRITVAL).C
404 C                                                 C
405 C  NOTE: these vectors MUST have been set up with C
406 C        LEV = 9 in the prior call to routine     C
407 C        HCASS.
408 C                                                 C
409 C  F. Murtagh, ESA/ESO/STECF, Garching, Feb. 1986.C
410 C                                                 C 
411 C-------------------------------------------------C
412       SUBROUTINE HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
413       include 'COMMON.IOUNITS'
414       CHARACTER*80 LINE
415       INTEGER IORDER(LEV),HEIGHT(LEV)
416       REAL CRITVAL(LEV)
417 c      INTEGER OUT(3*LEV,3*LEV)
418 c      INTEGER UP,ACROSS,BLANK
419       CHARACTER*1 OUT(3*LEV,3*LEV)
420       CHARACTER*1 UP,ACROSS,BLANK
421       DATA UP,ACROSS,BLANK/'|','-',' '/
422 C
423 C
424       DO I=1,3*LEV
425         DO J=1,3*LEV
426           OUT(I,J)=BLANK
427         ENDDO
428       ENDDO
429 C
430 C
431       DO I=3,3*LEV,3
432          I2=I/3
433 C
434          J2=3*LEV+1-3*HEIGHT(I2)
435          DO J=3*LEV,J2,-1
436             OUT(J,I)=UP
437          ENDDO
438 C
439          DO K=I,3,-1
440             I3=INT((K+2)/3)
441             IF ( (3*LEV+1-HEIGHT(I3)*3).LT.J2) GOTO 100
442             OUT(J2,K)=ACROSS
443          ENDDO
444   100    CONTINUE
445 C
446       ENDDO
447 C
448 C
449       IC=3
450       DO I=1,3*LEV
451       IF (I.EQ.IC+1) THEN
452                    IDUM=IC/3
453                    IDUM=LEV-IDUM
454                    DO L=1,LEV
455                       IF (HEIGHT(L).EQ.IDUM) GOTO 190
456                    ENDDO
457   190              IDUM=L
458                    WRITE(iout,200) CRITVAL(IDUM),(OUT(I,J),J=1,3*LEV)
459                    IC=IC+3
460                    ELSE
461                    LINE = ' '
462                    WRITE(iout,210) (OUT(I,J),J=1,3*LEV)
463       ENDIF
464   200 FORMAT(1H ,8X,F12.2,4X,27000A1)
465   210 FORMAT(1H ,24X,27000A1)
466       ENDDO
467       WRITE(iout,250)
468       WRITE(iout,220)(IORDER(J),J=1,LEV)
469       WRITE(iout,250)
470   220 FORMAT(1H ,24X,9000I3)
471       WRITE(iout,230) LEV
472   230 FORMAT(1H ,13X,'CRITERION        CLUSTERS 1 TO ',i3)
473       WRITE(iout,240) LEV-1
474   240 FORMAT(1H ,13X,'VALUES.      (TOP ',i3,' LEVELS OF HIERARCHY).')
475   250 FORMAT(/)
476 C
477 C
478       RETURN
479       END