C C*********************************************************************** C C Generalized Agglomerative Hierarchical Clustering Algorithms C Programmed By Glenn W. Milligan C Prepared March 1980 - The Ohio State University C Revised By Richard Cheng, October 1993 C C*********************************************************************** C C Complier Specific Information: C C The present software code makes calls to several subroutines or C functions that may not be supported by other FORTRAN compilers. C First, these complier specific routines can be safely disabled C with no adverse effect to the actual clustering computations. C All routines are directed at improving the user interface. C For example, consider CALL BEEP. This sounds a warning beep C when the user should be aware of a data entry error or problem. C The second is the INQUIRE() function that tests for the existence C of the specified input data file. C The PAUSE() function simply waits for the user to hit any C key. The function IARGC() counts the number of command line C arguments present when HCINFLU.EXE was started. C Screen graphics functions include: C SET_VIDEO_MODE() C GET_VIDEO_MODE() C GRAPHICS_MODE() C CLEAR() C HOME() C SET_PALETTE() C TEXT_MODE C CLEAR_TEXT C Again, these functions can be disabled when adopting the C program for a different complier. C C*********************************************************************** C C INCLUDE 'GREX.FH' PARAMETER(MAXPTS=200,MAXVAR=30) C C MAXPTS = Maximum number of data points C MAXVAR = Maximum number of variables C COMMON E(MAXPTS,MAXPTS),BETA,BETALW REAL LEVEL(MAXPTS) INTEGER OBJPOS(MAXPTS*2),PAUSE,ERR,OMODE,IARGV C C PAUSE, ERR, OMODE, IARGV are display graphics variables C DIMENSION TITLE(20),IP(MAXPTS),D(MAXPTS,MAXPTS) DIMENSION LK(MAXPTS),LKI(MAXPTS),LKJ(MAXPTS),XX(MAXPTS) DIMENSION I9(10),P(MAXPTS*2),PX(MAXPTS),GS(MAXPTS) DIMENSION GM(MAXPTS),PCORR(MAXPTS),DATA(MAXPTS,MAXPTS) DIMENSION R(MAXPTS,MAXPTS),RCORR(MAXPTS,MAXPTS),SP(MAXPTS) DIMENSION XR(MAXPTS),DR(MAXPTS),IQ(MAXPTS) DIMENSION LLSET(MAXPTS),LLSTOR(MAXPTS),U(MAXPTS*MAXVAR) LOGICAL W(MAXPTS),WD(MAXPTS,MAXPTS),TEST CHARACTER*1 LABEL(6,MAXPTS) CHARACTER*20 FILIN,OUTFIL CHARACTER*40 FORM CHARACTER*60 MC(10) CHARACTER PAL(17) DATA XXXX / 'X'/,BLANK / ' '/,PERIOD / '.'/ DATA MC /'Single Linkage Method','Complete Linkage Method', * 'Group Average Method','Weighted Average Method', * 'Centroid Method','Median Method', * 'Ward''s Minimum Variance Method', * 'Lance & William''s Beta-Flexible Method', * 'Belbin, Faith, & Milligan''s Beta-Flexible Method', * 'User Specified Method'/ DATA PAL /'4','3','3','3','3','3','3','3','3','3','3','3', * '3','3','3','3','4'/ C C CALL ICIDIS (2,'&','&','=',' ') C C USER INPUT SECTION C IARGV = IARGC() IF (IARGV.EQ.0) THEN OMODE = GET_VIDEO_MODE(I,J) ERR = GRAPHICS_MODE (16) ERR = CLEAR() ERR = HOME() ERR = SET_PALETTE(PAL) ELSE ERR = GRAPHICS_MODE (6) ERR = CLEAR() ERR = TEXT_MODE() ENDIF WRITE (*,10) 10 FORMAT (/,13X, * 'ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»', */,13X,'º HCINFLU º', */,13X,'º º', */,13X,'º Generalized Agglomerative Hierarchical Clustering º', */,13X,'º Algorithms With Internal Influence Measurement º', */,13X,'º º', */,13X,'º Developed by Glenn W. Milligan º', */,13X,'º Revised by Richard Cheng º', */,13X,'º Ohio State University º', */,13X,'º Version 2.0 -- Oct. 1993 º', */,13X,'ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ', *////,13X,' In the event of a data entry error,' */,13X,' program execution can be terminated by' */,13X,' entering -99 for any numerical prompt.' *////,13X,' Press enter to continue...') IKEY = PAUSE() IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF WRITE (*,*) 'Enter a one line title for your analysis:' READ (*,20) (TITLE(I),I=1,20) 20 FORMAT (20A4) MODIFY = 0 MODBETA = 0 MODBFM = 0 30 WRITE (*,*) 'How many elements in the data set?' 40 READ (*,50,ERR=60) N 50 FORMAT (I5) IF (N.EQ.-99) GO TO 2970 GO TO 80 60 CALL BEEP WRITE (*,70) 70 FORMAT (' Integer value expected, please re-enter:') GO TO 40 80 IF (N.GT.MAXPTS) THEN CALL BEEP WRITE (*,90) 90 FORMAT (' You have exceeded the maximum number of objects ', * 'allowed.',/,' Please re-enter:') GO TO 40 ENDIF IF (N.LT.2) THEN CALL BEEP IF (N.LT.0) THEN WRITE (*,100) 100 FORMAT (' You cannot have a negative number of objects!',/, * ' Please re-enter:') ELSE WRITE (*,110) 110 FORMAT (' You cannot cluster less than two objects!',/, * ' Please re-enter:') ENDIF GO TO 40 ENDIF IF (MODIFY.EQ.1.AND.ICLUS2.GT.(N-1)) GO TO 650 IF (MODIFY.EQ.1) GO TO 1070 120 WRITE (*,*) 'How many variables in the data set?' 130 READ (*,50,ERR=140) NV IF (NV.EQ.-99) GO TO 2970 GO TO 150 140 CALL BEEP WRITE (*,70) GO TO 130 150 IF (NV.GT.MAXVAR) THEN CALL BEEP WRITE (*,160) 160 FORMAT (' You have exceeded the maximum number of variables', * ' allowed.',/,' Please re-enter:') GO TO 130 ENDIF IF (NV.LT.1) THEN CALL BEEP IF (NV.LT.0) THEN WRITE (*,170) 170 FORMAT (' You cannot have a negative number of variables!',/ * ,' Please re-enter:') ELSE WRITE (*,180) 180 FORMAT (' You cannot cluster with less than one ', * 'variable!',/, * ' If you are using your own proximity matrix,',/, * ' enter the number of variables used to compute',/, * ' the proximities. If not known or not applicable,',/, * ' enter 1',//,' If you are using your own proximity ', * 'matrix',/, * ' and the number if variables exceeds 20, enter 20',//, * ' The correct number of variables is needed for the ', * 'test',/, * ' of clustering tendency based on the point-biserial',/, * ' correlation coefficient.',//,' Please re-enter:') ENDIF GO TO 130 ENDIF IF (MODIFY.EQ.1) GO TO 1070 190 WRITE (*,*) 'Select the type of proximity measure:' WRITE (*,*) ' 3 = Variables input data:' WRITE (*,*) ' User supplied proximity measure' WRITE (*,*) ' (Requires source code modification)' WRITE (*,*) ' 2 = Variables input data:' WRITE (*,*) ' Program to compute Pearson correlations' WRITE (*,*) ' 1 = Variables input data:' WRITE (*,*) ' Program to compute Euclidean distances' WRITE (*,*) ' 0 = User supplied proximity matrix:' WRITE (*,*) ' Lower-half triangular matrix' 200 READ (*,50,ERR=210) MTYPE IF (MTYPE.EQ.-99) GO TO 2970 GO TO 230 210 CALL BEEP WRITE (*,220) 220 FORMAT (' An integer between 0 and 3 is expected, please ', * 're-enter:') GO TO 200 230 IF (MTYPE.LT.0.OR.MTYPE.GT.3) THEN CALL BEEP WRITE (*,220) GO TO 200 ENDIF IF (MTYPE.EQ.2.AND.NV.LE.2) THEN CALL BEEP WRITE (*,240) 240 FORMAT (' This option requires 3 or more variables,', * ' select a different entry.') GO TO 200 ENDIF IF (MTYPE.EQ.2) IS = -1 250 IF (MTYPE.EQ.1) THEN IS = 1 WRITE (*,*) 'Standardize distances by range?' WRITE (*,*) ' 1 = Yes' WRITE (*,*) ' 0 = No (Default)' 260 READ (*,50,ERR=270) ISTD IF (ISTD.EQ.-99) GO TO 2970 GO TO 280 270 CALL BEEP WRITE (*,70) GO TO 260 280 IF (ISTD.NE.0.AND.ISTD.NE.1) THEN CALL BEEP WRITE (*,*) 'Entry error, please re-enter:' GO TO 260 ENDIF ENDIF MTYPE1 = MTYPE IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF IF (MTYPE.EQ.3.OR.MTYPE.EQ.0) THEN WRITE (*,*) 'Enter nature of the proximity measure:' WRITE (*,*) ' 1 = Dissimilarity Measure' WRITE (*,*) * ' (Smaller values represent greater similarity)' WRITE (*,*) ' -1 = Similarity Measure' WRITE (*,*) * ' (Larger values represent greater similarity)' 290 READ (*,50,ERR=300) IS IF (IS.EQ.-99) GO TO 2970 GO TO 310 300 CALL BEEP WRITE (*,70) GO TO 290 310 IF (IS.NE.1.AND.IS.NE.-1) THEN CALL BEEP WRITE (*,*) 'Entry error, please re-enter:' GO TO 290 ENDIF ENDIF IF (MTYPE.EQ.0) THEN WRITE (*,*) 'Note: Cluster centroids cannot be computed' WRITE (*,*) ' when a user supplied proximity matrix' WRITE (*,*) ' is used.' WRITE (*,*) ' ' WRITE (*,*) 'Press enter to continue...' WRITE (*,*) IKEY = PAUSE() ENDIF IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF IF (MODIFY.EQ.1) GO TO 1070 DO 320 I = 1, 10 I9(I) = 0 320 CONTINUE 330 ISKIP = 0 WRITE (*,*) 'Clustering method codes:' WRITE (*,*) ' 1-Single Linkage Method' WRITE (*,*) ' 2-Complete Linkage Method' WRITE (*,*) ' 3-Group Average Method' WRITE (*,*) ' 4-Weighted Average Method' WRITE (*,*) ' 5-Centroid Method' WRITE (*,*) ' 6-Median Method' WRITE (*,*) ' 7-Ward''s Minimum Variance Method' WRITE (*,*) ' 8-Lance & William''s Beta-Flexible Method' WRITE (*,340) 340 FORMAT (4X,' 9-Belbin, Faith, & Milligan''s Beta-Flexible', * ' Method') WRITE (*,350) 350 FORMAT (4X,'10-User Specified Method (Requires Coding', * ' in the Source Program)') WRITE (*,*) 'Enter code for the first selected method:' 360 IF (ISKIP.EQ.0) THEN 370 READ (*,50,ERR=380) ITMP IF (ITMP.EQ.-99) GO TO 2970 GO TO 400 380 CALL BEEP WRITE (*,390) 390 FORMAT (' An integer between 1 and 10 is expected, please ', * 're-enter:') GO TO 360 400 IF (ITMP.GT.10.OR.ITMP.LE.0) THEN CALL BEEP WRITE (*,410) 410 FORMAT (' That method code is invalid, please ','re-enter:') GO TO 370 ELSE I9(1) = ITMP I = 2 ENDIF ENDIF 420 WRITE (*,430) 430 FORMAT (' Enter next method code or press enter to end selection:' * ) 440 READ (*,50,ERR=450) ITMP IF (ITMP.EQ.-99) GO TO 2970 GO TO 460 450 CALL BEEP WRITE (*,70) GO TO 360 460 IF (ITMP.GT.10.OR.ITMP.LT.0) THEN CALL BEEP ISKIP = 1 WRITE (*,470) 470 FORMAT (' That method code is invalid, please re-enter:') GO TO 440 ELSE I9(I) = ITMP ENDIF IF (I9(I).EQ.0) GO TO 480 I = I+1 GO TO 420 IF (MODIFY.EQ.1) GO TO 1070 480 BETALW = 0.0 DO 570 I = 1, 10 IF (I9(I).EQ.8) THEN WRITE (*,490) 490 FORMAT (' Enter beta value for the Lance & William''s', * ' method',/,' (Suggested value: -0.25):') 500 READ (*,510,ERR=520) BETALW 510 FORMAT (G10.0) IF (BETALW.EQ.-99.) GO TO 2970 GO TO 530 520 BETALW = 10. 530 IF (BETALW.LE.-1.0.OR.BETALW.GE.1.0) THEN CALL BEEP WRITE (*,540) 540 FORMAT (' The beta value must be between -1.0 and 1.0.', * ' Please re-enter:') GO TO 500 ENDIF IF (BETALW.EQ.0.0) THEN CALL BEEP WRITE (*,550) 550 FORMAT (' Do you really mean 0.0 as the value for', * ' beta?',/,' 0 = Yes - use 0.0',/, * ' 1 = No - request new value',/, * ' 2 = No - use suggested value of -0.25') 560 READ (*,50) ITEST IF (ITEST.EQ.-99) GO TO 2970 IF (ITEST.NE.0.AND.ITEST.NE.1.AND.ITEST.NE.2) THEN CALL BEEP WRITE (*,*) 'Entry error, please re-enter:' GO TO 560 ENDIF IF (ITEST.EQ.1) GO TO 480 IF (ITEST.EQ.2) BETALW = -.25 ENDIF ENDIF 570 CONTINUE IF (MODBETA.EQ.1) GO TO 1070 580 BETA = 0.0 DO 640 I = 1, 10 IF (I9(I).EQ.9) THEN WRITE (*,590) 590 FORMAT (' Enter beta value for the Belbin, Faith, & ', * 'Milligan''s method',/,' (Suggested value: -0.10):') 600 READ (*,510,ERR=610) BETA IF (BETA.EQ.-99.) GO TO 2970 GO TO 620 610 BETA = 10. 620 IF (BETA.LE.-1.0.OR.BETA.GE.1.0) THEN CALL BEEP WRITE (*,*) * 'The beta value must fall between -1.0 to +1.0.' WRITE (*,*) 'Please re-enter:' GO TO 600 ENDIF IF (BETA.EQ.0.0) THEN CALL BEEP WRITE (*,550) 630 READ (*,50) ITEST IF (ITEST.EQ.-99) GO TO 2970 IF (ITEST.NE.0.AND.ITEST.NE.1.AND.ITEST.NE.2) THEN CALL BEEP WRITE (*,*) 'Entry error - please re-enter:' GO TO 630 ENDIF IF (ITEST.EQ.1) GO TO 580 IF (ITEST.EQ.2) BETA = -.10 ENDIF ENDIF 640 CONTINUE IF (MODBFM.EQ.1) GO TO 1070 IF (MODIFY.EQ.1) GO TO 1070 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF 650 WRITE (*,*) 'Enter the lower bound for the number of clusters' WRITE (*,*) 'for which detailed analyses are to be reported.' 660 READ (*,50,ERR=670) ICLUS1 IF (ICLUS1.EQ.-99) GO TO 2970 GO TO 680 670 CALL BEEP WRITE (*,70) GO TO 660 680 IF (ICLUS1.LT.2) THEN CALL BEEP WRITE (*,690) 690 FORMAT (' Minimum level is 2 clusters',/,' Please re-enter:') GO TO 660 ENDIF IF (ICLUS1.GT.(N-1)) THEN CALL BEEP N1 = N-1 WRITE (*,700) N1 700 FORMAT (' Maximum level is ',I3,' clusters',/, * ' Please re-enter:') GO TO 660 ENDIF WRITE (*,*) 'Enter the upper bound for the number of clusters' WRITE (*,*) 'for which detailed analyses are to be reported.' 710 READ (*,50,ERR=720) ICLUS2 IF (ICLUS2.EQ.-99) GO TO 2970 GO TO 730 720 CALL BEEP WRITE (*,70) GO TO 710 730 IF (ICLUS2.LT.2) THEN CALL BEEP WRITE (*,690) GO TO 710 ENDIF IF (ICLUS2.GT.(N-1)) THEN CALL BEEP N1 = N-1 WRITE (*,700) N1 GO TO 710 ENDIF C C MAKE SURE ICLUS2 IS THE UPPER BOUND OF RANGE - OTHERWISE SWITCH THEM C IF (ICLUS1.GT.ICLUS2) THEN ITEMP = ICLUS1 ICLUS1 = ICLUS2 ICLUS2 = ITEMP ENDIF IF (ICLUS1.EQ.ICLUS2) THEN WRITE (*,740) 740 FORMAT (' Since the upper and lower bounds match, the cluster', * /,' centroids will be written to file CENT.DAT using format' * ,/,' (30G15.5)',/) ENDIF IF (MODIFY.EQ.1) GO TO 1070 750 WRITE (*,*) 'How many columns does your output device support?' WRITE (*,*) '(Minimum = 75, Default = 80)' 760 READ (*,50,ERR=770) NCOL IF (NCOL.EQ.-99) GO TO 2970 GO TO 780 770 CALL BEEP WRITE (*,70) GO TO 760 780 IF (NCOL.EQ.0) NCOL = 80 IF (NCOL.LT.75) THEN CALL BEEP WRITE (*,*) 'The minimum number of output columns is 75' WRITE (*,*) 'Please re-enter:' GO TO 760 ENDIF IF (MODIFY.EQ.1) GO TO 1070 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF 790 WRITE (*,*) 'Test for the existence of cluster structure?' WRITE (*,*) ' 1 = Yes' WRITE (*,*) ' 0 = No (Default)' 800 READ (*,50,ERR=810) IR IF (IR.EQ.-99) GO TO 2970 GO TO 820 810 CALL BEEP WRITE (*,920) GO TO 800 820 IF (IR.NE.1.AND.IR.NE.0) THEN CALL BEEP WRITE (*,*) 'Entry error, please re-enter:' GO TO 800 ENDIF IF (IR.EQ.1) THEN CALL BEEP WRITE (*,830) 830 FORMAT (/,1X,'Warning:',/, * ' It may take a considerable amount of time ',/, * ' to conduct the test of clustering tendency.',/) ENDIF IF (IR.EQ.1.AND.(MTYPE.EQ.0.OR.MTYPE.EQ.3)) THEN WRITE (*,840) 840 FORMAT (' The program is configured to conduct the',/, * ' test based on Euclidean distances or',/, * ' Pearson correlations only.',/, * ' If your proximity matrix was computed from',/, * ' a different measure, then the test may not',/, * ' be valid. You need to modify the software',/, * ' to have the subroutine compute the measure',/, * ' that you used for the original data.',/) WRITE (*,*) 'Select the proximity measure to be used:' WRITE (*,*) ' 1 = Simulation based on distances' WRITE (*,*) ' 2 = Simulation based on correlations' 850 READ (*,50,ERR=860) MTYPE1 IF (MTYPE1.EQ.-99) GO TO 2970 GO TO 880 860 CALL BEEP WRITE (*,870) 870 FORMAT (' Please select 1 or 2:') GO TO 850 880 IF (MTYPE1.LT.1.OR.MTYPE1.GT.2) THEN CALL BEEP WRITE (*,870) GO TO 850 ENDIF IF (MTYPE1.EQ.2.AND.NV.LE.2) THEN CALL BEEP WRITE (*,240) GO TO 850 ENDIF WRITE (*,*) ' ' WRITE (*,*) 'Press enter to continue...' IKEY = PAUSE() IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF ENDIF IF (MODIFY.EQ.1) GO TO 1070 890 WRITE (*,*) 'Internal influence measure for each data point?' WRITE (*,*) ' 1 = Yes' WRITE (*,*) ' 0 = No (Default)' 900 READ (*,50,ERR=910) INF IF (INF.EQ.-99) GO TO 2970 GO TO 930 910 CALL BEEP WRITE (*,920) 920 FORMAT (' Please type 1 or press enter:') GO TO 900 930 IF (INF.NE.1.AND.INF.NE.0) THEN CALL BEEP WRITE (*,*) 'Entry error, please re-enter:' GO TO 900 ENDIF IF (INF.EQ.1.AND.N.GT.75) THEN CALL BEEP WRITE (*,940) 940 FORMAT (/,1X,'Warning:',/, * ' It may take a considerable amount of time ',/, * ' to compute the internal influence measures.',/) WRITE (*,*) ' ' WRITE (*,*) 'Press enter to continue...' WRITE (*,*) IKEY = PAUSE() ENDIF IF (MODIFY.EQ.1) GO TO 1070 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF 950 WRITE (*,*) 'Tree output option:' WRITE (*,*) ' 0 = Skyline plot (Default)' WRITE (*,*) ' 1 = Icicle plot' WRITE (*,*) ' 2 = No tree output' 960 READ (*,50,ERR=970) IPLOT IF (IPLOT.EQ.-99) GO TO 2970 GO TO 980 970 CALL BEEP WRITE (*,70) GO TO 960 980 IF (IPLOT.NE.0.AND.IPLOT.NE.1.AND.IPLOT.NE.2) THEN CALL BEEP WRITE (*,*) 'Entry error, please re-enter:' GO TO 960 ENDIF IF (IPLOT.NE.1) GO TO 990 C c- c- Read in parameters c- nlab: max number of char.s in any label, c- but if 0, labels absent, c- or if -1, use same labels as with preceding set of data. c- klev: level options: -1, 0, 1, 2, 3 (first arg to "icilev"). c- 0=use log spacing; -1=same, and return display levels; c- 1=user provides levels; 2,3=one level for each join c- Levels (sec. arg to "icilev") are present only if klev=1 c- nlev: number of levels, but can use 0 and let program decide c- (used if klev=-1,0,1) (third arg to "icilev") c- kobjpo: if 1, user gives obj. positions; if 0, user doesn't; c- if -1, return the object positions used c- and choice of submethods: any subset of 1,2,3, in any order c- (1=single-link, 2=complete-link, 3=average submethod) c- C c write (*,*)' ENTER THE MAXIMUM # OF CHARACTERS IN ANY LABEL:' c WRITE (*,*) ' N = NUMBER OF CHARACTERS (6 MAXIMUM)' c write (*,*) ' 0 = LABELS ARE ABSENT' c write (*,*) ' -1 = USE THE SAME LABELS FROM THE' c write (*,*) ' PRECEDING DATA SET' c read (*,*) nlab C NLAB = 2 IF (N.GE.100) NLAB = 3 C c write (*,*) 'SPACING LEVEL OPTIONS:' c write (*,*) ' 2 = ONE LEVEL FOR EACH MERGER' c write (*,*) ' (RECOMMENDED)' c write (*,*) ' 1 = USER PROVIDED LEVELS' c write (*,*) ' 0 = LOG SPACING' c write (*,*) ' -1 = LOG SPACING & LIST LEVELS' c read (*,*) klev C KLEV = 2 C c write (*,*) 'NUMBER OF LEVELS OPTION:' c write (*,*) ' 0 = PROGRAM DETERMINES NUMBER OF LEVELS' c write (*,*) ' (RECOMMENDED)' c write (*,*) ' N > 0 USER SPECIFIES A NUMBER GREATER THAN 0' c read (*,*) nlev C NLEV = 0 C c write (*,*) 'OBJECT POSITIONS:' c write (*,*) ' 0 = PROGRAM SELECTS OBJECT POSITIONS ' c write (*,*) ' (RECOMMENDED)' c write (*,*) ' 1 = USER SPECIFIED OBJECT POSITIONS' c read (*,*) kobjpo c read (*,111), nlab, klev, nlev, kobjpo c111 format(7i3) C c- c- Read new labels, or use old labels, or use no labels c if(.not.(nlab.gt.0))goto 23015 c nlabls=nlab c write (*,*) 'ENTER NAME OF FILE CONTAINING LABELS' c read (*,2) infil2 c open (7,file=infil2,status='unknown') c write (*,*) "ENTER THE FORMAT FOR READING THE LABELS" c read (*,2) (form1(i),i=1,20) c read (7,form1) ((label(i,j),i=1,nlab),j=1,n) c write (6,*), c & "LABELS WERE READ IN ACCORDING TO THE FOLLOWING FORMAT" c write (6,23) form1 c23 format(1x,20a4) C C c if(.not.(nlab.lt.0))goto 23017 c nlab=nlabls C C c if(.not.(nlab.eq.0))goto 23019 c nlabls=0 C c- If levels are present, read them in C C c if(.not.(klev.eq.1))goto 23021 c read (*,2) form1 c read (*,form1) (level(i),i=1,nlev) c write (*,*) c & "LEVELS WERE READ IN ACCORDING TO THE FOLLOWING FORMAT" c write (*,23) form1 c c- If the user is supplying object positions, read them in C C c if(.not.(kobjpo.eq.1))goto 23023 c read (*,2) form1 c read (*,form1) (objpos(i),i=1,n) c write (*,*) c &"OBJECT POSITIONS WERE READ IN ACCORDING TO THE FOLLOWING FORMAT" c write (*,23) form1 C 990 CONTINUE IF (MODIFY.EQ.1) GO TO 1070 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF 1000 WRITE (*,*) 'Enter the data file name:' 1010 READ (*,1020) FILIN 1020 FORMAT (A20) INQUIRE (FILE=FILIN,EXIST=TEST) IF (TEST) GO TO 1040 CALL BEEP WRITE (*,1030) 1030 FORMAT (' Input file not found - Try again...') GO TO 1010 1040 IF (MODIFY.EQ.1) GO TO 1070 1050 WRITE (*,*) 'Enter the data file input format:' READ (*,*) FORM IF (MODIFY.EQ.1) GO TO 1070 1060 WRITE (*,*) 'Enter the output file name:' READ (*,1020) OUTFIL 1070 IST1 = 0 IST2 = 0 MODBETA = 0 MODBFM = 0 DO 1080 I = 1, 10 IF (I9(I).EQ.0) THEN NOFM = I-1 GO TO 1090 ENDIF IF (I9(I).EQ.8) IST1 = 1 IF (I9(I).EQ.9) IST2 = 1 1080 CONTINUE NOFM = 10 C C USER CONFIRMATION SECTION C 1090 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF WRITE (*,1100) 1100 FORMAT (/,6X,'You have entered the following information:',/) WRITE (*,1110) N,NV 1110 FORMAT (5X,' 1. Number of objects: ',I3,/,5X, * ' 2. Number of variables: ',I3) IF (MTYPE.EQ.1) THEN WRITE (*,1120) 1120 FORMAT (5X,' 3. Proximity measure: Program to compute ', * 'Euclidean distances') IF (ISTD.EQ.1) THEN WRITE (*,1130) 1130 FORMAT (5X,' 4. Distances standardized by range: Yes') ELSE WRITE (*,1140) 1140 FORMAT (5X,' 4. Distances standardized by range: No') ENDIF ENDIF IF (MTYPE.EQ.2) THEN WRITE (*,1150) 1150 FORMAT (5X,' 3. Proximity measure: Program to compute ', * 'Pearson correlations') ENDIF IF (MTYPE.EQ.0.OR.MTYPE.EQ.3) THEN IF (IS.EQ.1) THEN WRITE (*,1160) 1160 FORMAT (5X,' 3. Proximity measure: Dissimilarity measure') ELSE WRITE (*,1170) 1170 FORMAT (5X,' 3. Proximity measure: Similarity measure') ENDIF ENDIF WRITE (*,1180) (I9(I),I=1,NOFM) 1180 FORMAT (5X,' 5. Method codes requested : ',10(I3,2X)) IF (IST1.EQ.1) WRITE (*,1190) BETALW 1190 FORMAT (5X,' 6. Beta value for Lance & William''s method',' = ', * F5.2) IF (IST2.EQ.1) WRITE (*,1200) BETA 1200 FORMAT (5X,' 7. Beta value for Belbin, Faith & Milligan''s', * ' method = ',F5.2) IF (ICLUS1.NE.ICLUS2) THEN WRITE (*,1210) ICLUS1,ICLUS2 1210 FORMAT (5X,' 8. Detailed analysis for the',I3,' to',I3, * ' cluster solution') ELSE WRITE (*,1220) ICLUS1 1220 FORMAT (5X,' 8. Detailed analysis for the',I3, * ' cluster solution') ENDIF WRITE (*,1230) NCOL 1230 FORMAT (5X,' 9. Number of columns of output device:',I4) IF (IR.EQ.0) THEN WRITE (*,1240) 1240 FORMAT (5X,' 10. Test existence of cluster structure: No') ELSE WRITE (*,1250) 1250 FORMAT (5X,' 10. Test existence of cluster structure: Yes') ENDIF IF (INF.EQ.0) THEN WRITE (*,1260) 1260 FORMAT (5X,' 11. Internal influence for each point: No') ELSE WRITE (*,1270) 1270 FORMAT (5X,' 11. Internal influence for each point: Yes') ENDIF IF (IPLOT.EQ.0) THEN WRITE (*,1280) 1280 FORMAT (5X,' 12. Plot output selected: Skyline plot') ELSE IF (IPLOT.EQ.1) THEN WRITE (*,1290) 1290 FORMAT (5X,' 12. Plot output selected: Icicle plot') ELSE WRITE (*,1300) 1300 FORMAT (5X,' 12. Plot output selected: None') ENDIF ENDIF WRITE (*,1310) FILIN 1310 FORMAT (5X,' 13. Data file name: ',A20) WRITE (*,1320) FORM 1320 FORMAT (5X,' 14. Data file format: ',A40) WRITE (*,1330) OUTFIL 1330 FORMAT (5X,' 15. Output file name: ',A20) WRITE (*,1340) 1340 FORMAT (//,5X,' Select an item number to make a change', * ' or press enter to continue...') 1350 READ (*,50,ERR=1360) ICH IF (ICH.EQ.-99) GO TO 2970 GO TO 1370 1360 CALL BEEP WRITE (*,70) GO TO 1350 1370 IF (ICH.NE.0) THEN MODIFY = 1 IF (ICH.EQ.6) MODBETA = 1 IF (ICH.EQ.7) MODBFM = 1 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF GO TO (30,120,190,250,330,480,580,650,750,790,890,950,1000,1050 * ,1060) ICH ENDIF OPEN (5,FILE=FILIN,STATUS='UNKNOWN') OPEN (6,FILE=OUTFIL,STATUS='UNKNOWN') WRITE (6,1380) 1380 FORMAT (28X,'HCINFLU',//,5X, * ' Generalized Agglomerative Hierarchical Clustering',/,7X, * ' Algorithms with Internal Influence Measurement',//,15X, * ' Developed by Glenn W. Milligan',/,19X, * ' Revised by Richard Cheng',/,21X,'Ohio State University',/,18X * ,' Version 2.0 -- Oct. 1993',//) WRITE (6,1390) TITLE 1390 FORMAT (1X,20A4) NRMAX = 100 IF (MTYPE.EQ.1.OR.MTYPE.EQ.2) THEN WRITE (6,1400) FILIN,N,NV,MTYPE,(I9(I),I=1,NOFM) 1400 FORMAT (/,' User Input Specifications:',//, * ' Data file name = ',A20,//, * ' Number of objects = ',I3,//, * ' Number of variables = ',I3,//, * ' Proximity measure type = ',I3,//, * ' Method codes requested = ',10(I3,2X)) ELSE WRITE (6,1410) FILIN,N,NV,MTYPE,IS,(I9(I),I=1,NOFM) 1410 FORMAT (/,' User Input Specifications:',//, * ' Data file name = ',A20,//, * ' Number of objects = ',I3,//, * ' Number of variables = ',I3,//, * ' Proximity measure type = ',I3,//, * ' Nature of proximity measure = ',I3,//, * ' Method codes requested = ',10(I3,2X)) ENDIF IF (IST1.EQ.1) WRITE (6,1420) BETALW 1420 FORMAT (/,' Beta value for Lance & William''s method',' = ',F5.2) IF (IST2.EQ.1) WRITE (6,1430) BETA 1430 FORMAT (/,' Beta value for Belbin, Faith & Milligan''s', * ' method = ',F5.2) IF (ICLUS1.NE.ICLUS2) THEN WRITE (6,1440) ICLUS1,ICLUS2 1440 FORMAT (/,' Detailed analysis for the',I3,' to',I3, * ' cluster solution') ELSE WRITE (6,1450) ICLUS1 1450 FORMAT (/,' Detailed analysis for the',I3,' cluster solution') ENDIF WRITE (6,1460) FORM 1460 FORMAT (/,' The data file format is ',A40) IF (ISTD.EQ.1) THEN WRITE (6,1470) 1470 FORMAT (/,' Range standardization of variables selected') ENDIF IF (IPLOT.EQ.0) THEN WRITE (6,1480) 1480 FORMAT (/,1X,'Skyline plot output selected') ELSE IF (IPLOT.EQ.1) THEN WRITE (6,1490) 1490 FORMAT (/,1X,'Icicle plot output selected') ELSE WRITE (6,1500) 1500 FORMAT (/,1X,'Tree output omitted') ENDIF ENDIF IF (IR.EQ.1.AND.MTYPE.EQ.0) THEN WRITE (6,1510) 1510 FORMAT (' The program is configured to conduct the test of',/, * ' clustering tendency using euclidean distances or',/, * ' Pearson correlations only.',/, * ' If your proximity matrix was computed from',/, * ' a different measure, then the test may not',/, * ' be valid. You will have to modify the software',/, * ' to have the subroutine compute the measure',/, * ' that you used for the original data.',/) ENDIF c if(iiplot.eq.1) then c write (6,21) nlab, klev, nlev, kobjpo c21 format(" ICICLE PLOT CONTROL VALUES THAT WERE READ IN:",//, c & " MAXIMUM NUMBER OF LABEL CHARACTERS =",i2,//, c & " LEVEL OPTION =", i3,//, " NUMBER OF LEVELS =", i3,//, c & " OBJECT POSITION OPTION =", i3) c endif C LEV1 = N-ICLUS1 LEV2 = N-ICLUS2 S = IS N1 = N-1 N2 = N1-1 XN = N DO 1520 I = 1, N XR(I) = 0.0 DR(I) = 0.0 E(I,I) = 0.0 1520 CONTINUE IF (MTYPE.NE.0) THEN READ (5,FORM,ERR=2960,END=2950) ((DATA(I,J),I=1,NV),J=1,N) WRITE (6,1530) 1530 FORMAT (/,1X,'The first two data lines were read as follows:',/ * ) DO 1540 J = 1, 2 WRITE (6,FORM) (DATA(I,J),I=1,NV) 1540 CONTINUE WRITE (6,*) ' ' ENDIF IF (MTYPE.EQ.1) THEN IF (ISTD.EQ.1) THEN DO 1570 I = 1, NV XMAX = 0 XMIN = 9.9E99 DO 1550 J = 1, N IF (DATA(I,J).GT.XMAX) XMAX = DATA(I,J) IF (DATA(I,J).LT.XMIN) XMIN = DATA(I,J) 1550 CONTINUE RANGE = XMAX-XMIN DO 1560 J = 1, N DATA(I,J) = (DATA(I,J)-XMIN)/RANGE 1560 CONTINUE 1570 CONTINUE ENDIF DO 1600 I = 1, N1 II = I+1 DO 1590 J = II, N SUM = 0.0 DO 1580 K = 1, NV SUM = SUM+(DATA(K,J)-DATA(K,I))**2 1580 CONTINUE E(I,J) = SQRT(SUM) E(J,I) = E(I,J) 1590 CONTINUE 1600 CONTINUE ENDIF C C CALCULATE PAIRWISE CORRELATION C IF (MTYPE.EQ.2) THEN DO 1650 I = 1, N-1 SUMX = 0. SUMX2 = 0. DO 1610 K = 1, NV SUMX = SUMX+DATA(K,I) SUMX2 = SUMX2+DATA(K,I)**2 1610 CONTINUE DO 1640 J = I+1, N SUMY = 0. SUMY2 = 0. SUMXY = 0. DO 1620 K = 1, NV SUMY = SUMY+DATA(K,J) SUMY2 = SUMY2+DATA(K,J)**2 SUMXY = SUMXY+DATA(K,I)*DATA(K,J) 1620 CONTINUE SXY = SUMXY-(SUMX*SUMY/NV) SXX = SUMX2-((SUMX**2)/NV) SYY = SUMY2-((SUMY**2)/NV) IF (ABS(SXX).LT.0.000001) THEN CALL BEEP WRITE (*,1630) I 1630 FORMAT (' A correlation is undefined because', * ' data point ,',I3,' contains all constants.',/, * ' Program terminating.') WRITE (*,*) ' ' WRITE (*,*) 'Press enter to continue...' WRITE (*,*) IKEY = PAUSE() GO TO 2970 ENDIF IF (ABS(SYY).LT.0.000001) THEN CALL BEEP WRITE (*,1630) J WRITE (*,*) 'Press enter to continue...' WRITE (*,*) IKEY = PAUSE() GO TO 2970 ENDIF E(I,J) = SXY/(SQRT(SXX*SYY)) E(J,I) = E(I,J) 1640 CONTINUE 1650 CONTINUE ENDIF IF (MTYPE.EQ.3) THEN C C*********************************************************************** C Compute Proximity Measure C If a different distance or similarity measure is desired C replace the following block with an appropriate routine. C C*********************************************************************** C C The following 6 lines will need to be removed or commented: C CALL BEEP WRITE (*,*) ' ' WRITE (*,*) 'Source code has not been modified for' WRITE (*,*) 'user specified proximity measure -' WRITE (*,*) 'Program terminating.' GO TO 2970 ENDIF C C READ LOWER TRIANGLE MATRIX FROM DATA FILE C IF (MTYPE.EQ.0) THEN DO 1660 I = 2, N I1 = I-1 READ (5,FORM,ERR=2960,END=2950) (E(I,J),J=1,I1) DO 1660 J = 1, I1 E(J,I) = E(I,J) 1660 CONTINUE WRITE (6,1530) DO 1670 I = 2, 3 I1 = I-1 WRITE (6,FORM) (E(I,J),J=1,I1) 1670 CONTINUE WRITE (6,*) ' ' ENDIF IKOUNT = 1 K = I9(1) JR = 0 C C Random Number Seed-Value Can Be Changed As Long As It Is An Odd C Integer Value. The variable IX is the seed. C IX = 5131979 IF (NV.LT.1) NV = MAXVAR NVN = NV*N 1680 IF (IARGV.EQ.0) THEN ERR = CLEAR() ERR = HOME() ELSE CALL CLEAR_TEXT ENDIF WRITE (*,1690) MC(K) 1690 FORMAT (/,' Starting the ',A60) IF (IR.EQ.1) JR = 1 C C WRITE OUT METHOD CODE TO OUTPUT FILE C WRITE (6,1700) MC(K) 1700 FORMAT (1X,'Starting the ',A60) GO TO 1800 C C Do Random Permutations For Significance Testing On Internal Criterion C 1710 CONTINUE WRITE (*,1720) NR,NRMAX 1720 FORMAT ('+---> Simulating data set:',I4,' out of',I4) CALL RANDOM (IX,NVN,U) C C*********************************************************************** C Compute Distance Measure C If a different distance or similarity measure is desired C replace the following block with an appropriate routine. C All code should be completed within the 1750 DO-LOOP. C*********************************************************************** C IF (MTYPE1.EQ.1) THEN DO 1750 I = 1, N1 II = I+1 DO 1740 J = II, N SUM = 0.0 DO 1730 M = 1, NV IZORK = (I-1)*NV+M JZORK = (J-1)*NV+M SUM = SUM+((U(IZORK)-U(JZORK)))**2 1730 CONTINUE R(I,J) = SQRT(SUM) R(J,I) = R(I,J) D(I,J) = R(I,J) D(J,I) = D(I,J) 1740 CONTINUE 1750 CONTINUE ENDIF C PAIRWISE CORRELATION OF SIMULATED DATA C IF (MTYPE1.EQ.2) THEN DO 1790 I = 1, N-1 SUMX = 0. SUMX2 = 0. DO 1760 K = 1, NV SUMX = SUMX+U((I-1)*NV+K) SUMX2 = SUMX2+U((I-1)*NV+K)**2 1760 CONTINUE DO 1780 J = I+1, N SUMY = 0. SUMY2 = 0. SUMXY = 0. DO 1770 K = 1, NV SUMY = SUMY+U((J-1)*NV+K) SUMY2 = SUMY2+U((J-1)*NV+K)**2 SUMXY = SUMXY+U((I-1)*NV+K)*U((J-1)*NV+K) 1770 CONTINUE SXY = SUMXY-(SUMX*SUMY/NV) SXX = SUMX2-((SUMX**2)/NV) SYY = SUMY2-((SUMY**2)/NV) D(I,J) = SXY/(SQRT(SXX*SYY)) D(J,I) = D(I,J) R(I,J) = D(I,J) R(J,I) = R(I,J) 1780 CONTINUE 1790 CONTINUE ENDIF GO TO 1820 1800 CONTINUE DO 1810 I = 1, N DO 1810 J = 1, N D(I,J) = E(I,J)*S 1810 CONTINUE 1820 DO 1830 IIXP = 1, N W(IIXP) = .FALSE. LK(IIXP) = 0 P(IIXP) = 1. GM(IIXP) = 0.0 DO 1830 IIPJ = 1, N WD(IIPJ,IIXP) = .FALSE. 1830 CONTINUE C C Begin Clustering Program C DO 2010 II = 1, N1 DO 1850 NI = 2, N IF (W(NI)) GO TO 1850 NI1 = NI-1 DO 1840 NJ = 1, NI1 IF (.NOT.W(NJ)) GO TO 1860 1840 CONTINUE 1850 CONTINUE GO TO 2970 1860 X = D(NI,NJ) DO 1880 I = 2, N IF (W(I)) GO TO 1880 I1 = I-1 DO 1870 J = 1, I1 IF (W(J)) GO TO 1870 IF ((D(I,J)-X).GT.0.0) GO TO 1870 NI = I NJ = J X = D(I,J) 1870 CONTINUE 1880 CONTINUE XX(II) = X LKI(II) = NI LKJ(II) = NJ W(NJ) = .TRUE. DO 2000 L = 1, N IF (W(L)) GO TO 2000 IF (L.EQ.NI) GO TO 2000 C C Branch to Method C GO TO (1890,1900,1910,1920,1930,1940,1950,1960,1970,1980), K 1890 IF ((D(NI,L)-D(NJ,L)).LE.0.0) GO TO 1990 D(NI,L) = D(NJ,L) D(L,NI) = D(L,NJ) GO TO 1990 1900 IF ((D(NI,L)-D(NJ,L)).GE.0.0) GO TO 1990 D(NI,L) = D(NJ,L) D(L,NI) = D(L,NJ) GO TO 1990 1910 D(L,NI) = (P(NI)/(P(NI)+P(NJ)))*D(L,NI)+(P(NJ)/(P(NI)+P(NJ)) * )*D(L,NJ) D(NI,L) = D(L,NI) GO TO 1990 1920 D(L,NI) = .5*D(L,NI)+.5*D(L,NJ) D(NI,L) = D(L,NI) GO TO 1990 1930 D(L,NI) = (P(NI)/(P(NI)+P(NJ)))*D(L,NI)+(P(NJ)/(P(NI)+P(NJ)) * )*D(L,NJ)-((P(NI)*P(NJ))/((P(NI)+P(NJ))**2.))*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 1990 1940 D(L,NI) = .5*D(L,NI)+.5*D(L,NJ)-.25*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 1990 1950 D(L,NI) = ((P(NI)+P(L))/(P(NI)+P(NJ)+P(L)))*D(L,NI)+((P(NJ)+ * P(L))/(P(NI)+P(NJ)+P(L)))*D(L,NJ)-(P(L)/(P(NI)+P(NJ)+P(L) * ))*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 1990 1960 D(L,NI) = .5*(1.-BETALW)*D(L,NI)+.5*(1.-BETALW)*D(L,NJ)+ * BETALW*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 1990 1970 D(L,NI) = (1.-BETA)*(P(NI)/(P(NI)+P(NJ)))*D(L,NI)+(1.-BETA)* * (P(NJ)/(P(NI)+P(NJ)))*D(L,NJ)+BETA*(D(NI,NJ)) D(NI,L) = D(L,NI) GO TO 1990 1980 CONTINUE C C Specify Optional Clustering Method Here In Terms Of D(L,NI) C See Cormack (1971) JRSS Series A C Model your source code after that shown for the previous options C C D(NI,L)=D(L,NI) C C The following 6 lines will need to be removed or commented: C CALL BEEP WRITE (*,*) ' ' WRITE (*,*) 'Source code has not been modified for' WRITE (*,*) 'user specified clustering method -' WRITE (*,*) 'Program terminating.' GO TO 2970 1990 CONTINUE 2000 CONTINUE P(NI) = P(NI)+P(NJ) 2010 CONTINUE LK(1) = LKJ(N1) LK(2) = LKI(N1) DO 2050 II = 2, N1 JJ = N-II DO 2020 KK = 1, II IF (LKI(JJ).EQ.LK(KK)) GO TO 2030 2020 CONTINUE GO TO 2970 2030 DO 2040 IK = KK, II KI = II+KK-IK LK(KI+1) = LK(KI) 2040 CONTINUE LK(KK) = LKJ(JJ) 2050 CONTINUE DO 2100 I = 1, N1 DO 2060 JI = 1, N IF (LK(JI).EQ.LKI(I)) GO TO 2070 2060 CONTINUE GO TO 2970 2070 DO 2080 JJ = 1, N IF (LK(JJ).EQ.LKJ(I)) GO TO 2090 2080 CONTINUE GO TO 2970 2090 LKI(I) = JI LKJ(I) = JJ 2100 CONTINUE DO 2120 I = 2, N2 I1 = I-1 DO 2110 J = 1, I1 JJ = I-J IF (LKJ(I).NE.LKI(JJ)) GO TO 2110 LKJ(I) = LKJ(JJ) GO TO 2120 2110 CONTINUE 2120 CONTINUE XX(N) = XX(N1) IF (IPLOT.NE.0) GO TO 2300 C C Skyline Plot C IF (IR.EQ.1.AND.JR.EQ.0) GO TO 2300 C C Nloops Is The Number Of Loops Required To Complete The Plot. Ncol Is C The Width Of The Output Device User Entered Earlier C WRITE (*,2130) 2130 FORMAT (/,' Generating skyline plot, please wait...') NPTS = (NCOL-15)/2 NLOOP = N/NPTS RLOOP = FLOAT(N)/FLOAT(NPTS) IF ((RLOOP-NLOOP).GT.0.001) NLOOP = NLOOP+1 IF (NLOOP.GT.1) THEN WRITE (6,2140) 2140 FORMAT (/,' The entire Skyline plot does not fit on a single', * ' page,',/,' Please place the pages side by side from', * ' left to right.') ENDIF DO 2290 LOOP = 1, NLOOP M1 = (LOOP-1)*NPTS+1 M2 = MIN(LOOP*NPTS,N) IP1 = (LOOP-1)*2*NPTS+1 IP2 = MIN(LOOP*2*NPTS,2*N-1) WRITE (6,2150) 12 2150 FORMAT (A1) IF (NLOOP.GT.1) THEN WRITE (6,2160) LOOP,NLOOP 2160 FORMAT (' Tree output strip ',I3,' of ',I3,' strips',/) ENDIF DO 2170 I = M1, M2 IP(I) = LK(I)/100 IQ(I) = IP(I) 2170 CONTINUE IF (N.GE.100) THEN WRITE (6,2180) (IP(I),I=M1,M2) 2180 FORMAT (' Item Number:',2X,50(1X,I1)) ENDIF DO 2190 I = M1, M2 IP(I) = (LK(I)-100*IP(I))/10 2190 CONTINUE IF (N.GE.100) THEN WRITE (6,2200) (IP(I),I=M1,M2) 2200 FORMAT (15X,50(1X,I1)) ELSE WRITE (6,2210) (IP(I),I=M1,M2) 2210 FORMAT (' Item Number:',2X,50(1X,I1)) ENDIF DO 2220 I = M1, M2 IP(I) = LK(I)-100*IQ(I)-10*IP(I) 2220 CONTINUE WRITE (6,2230) (IP(I),I=M1,M2) 2230 FORMAT (15X,50(1X,I1)) M69 = N*2 DO 2240 I = 1, M69, 2 P(I+1) = BLANK P(I) = PERIOD 2240 CONTINUE WRITE (6,2250) 2250 FORMAT (' Level') DO 2280 I = 1, N1 JI = 2*LKI(I)-1 JJ = 2*LKJ(I)-1 DO 2260 KK = JJ, JI P(KK) = XXXX 2260 CONTINUE IF ((XX(I).EQ.XX(I+1)).AND.(I.NE.N1)) GO TO 2280 NSUBI = N-I WRITE (6,2270) NSUBI,(P(III),III=IP1,IP2) 2270 FORMAT (1X,I4,11X,100A1) 2280 CONTINUE 2290 CONTINUE WRITE (6,2150) 12 C C Initialize For Internal Criterion Statistics C 2300 CONTINUE C C IF (IR.EQ.1.AND.JR.EQ.0) GO TO 2370 C c- call ICICLE print out routine C IF (IPLOT.NE.1) GO TO 2370 C c+If klev=-1, then use log spacing C KLEV0 = MAX(KLEV,0) C c+supply level information to ICICLE C CALL ICILEV (KLEV0,LEVEL,NLEV) IF (.NOT.(KOBJPO.EQ.1)) GO TO 2310 C c+obj. information to ICICLE C CALL ICIPOS (KOBJPO,OBJPOS,N) c- Arrays lki, lkj (merge pointers) are entered as arguments c- in reverse alphabetical order so that objects appear in c- ascending order in the ICICLE plot. C 2310 CONTINUE CALL ICICLE (N,LKJ,LKI,XX,S,6,NLAB,LABEL,TITLE,LK) C c- call icilev if levels are to be returned C IF (.NOT.(KLEV.EQ.-1)) GO TO 2330 CALL ICILEV (KLEV,LEVEL,NWHAT) WRITE (6,*) ' ' WRITE (6,*) 'The levels used by Icicle are:' WRITE (6,2320) (LEVEL(I),I=1,NWHAT) 2320 FORMAT (F12.3) C c- call icipos if positions are to be returned C 2330 CONTINUE IF (.NOT.(KOBJPO.EQ.-1)) GO TO 2350 CALL ICIPOS (KOBJPO,OBJPOS,N) WRITE (6,*) ' ' WRITE (6,*) 'The object positions used by Icicle are:' WRITE (6,2340) (OBJPOS(I),I=1,N) 2340 FORMAT (I6) c- Go to next method. C 2350 CONTINUE C c+End of long LOOP over methods, indexed by k. Bracket depth 2. C c+reset all modes to default values C CALL ICIRES WRITE (6,2360) 12 2360 FORMAT (A1) 2370 CONTINUE C C DO 2380 IIXP = 2, N I1 = IIXP-1 DO 2380 IIPJ = 1, I1 D(IIPJ,IIXP) = 0.0 2380 CONTINUE NIC1 = (N*(N-1))/2 DO 2730 IIPX1 = 1, N2 JJX = LKJ(IIPX1) JX = LKI(IIPX1) JX1 = JX-1 DO 2400 IIXP2 = JJX, JX1 JJX1 = IIXP2+1 DO 2390 I1I = JJX1, JX JPL1 = MIN0(LK(I1I),LK(IIXP2)) JPL2 = MAX0(LK(I1I),LK(IIXP2)) IF (D(JPL1,JPL2).NE.0.0) GO TO 2390 D(JPL1,JPL2) = XX(IIPX1) WD(JPL1,JPL2) = .TRUE. 2390 CONTINUE 2400 CONTINUE C C Compute Cluster Centroids C IF (IIPX1.LT.LEV2) GO TO 2670 IF (IIPX1.GT.LEV1) GO TO 2670 IF (IR.EQ.1.AND.JR.EQ.0) GO TO 2670 INPX1 = N-IIPX1 WRITE (*,2410) INPX1 2410 FORMAT (/,' Conducting analyses of cluster level ',I2) DO 2420 I = 2, N IP(I) = 0 2420 CONTINUE NCR = 1 II = 1 J = 2 IP(1) = 1 2430 DO 2440 JJ = J, N IF (.NOT.WD(II,JJ)) GO TO 2440 IP(JJ) = NCR GO TO 2480 2440 CONTINUE 2450 NCR = NCR+1 IF (NCR.GT.INPX1) GO TO 2490 DO 2460 JJ = 1, N IF (IP(JJ).EQ.0) GO TO 2470 2460 CONTINUE GO TO 2970 2470 IP(JJ) = NCR 2480 II = JJ J = II+1 IF (II.GE.N) GO TO 2450 GO TO 2430 2490 CONTINUE WRITE (6,2500) INPX1,INPX1 2500 FORMAT (//,' Cluster membership of data points for the ',I3, * ' cluster solution',/' (Hierarchy level ',I3,')',/) WRITE (6,2510) 2510 FORMAT (' Point Assigned to cluster',/,1X,30('-')) DO 2530 I = 1, N WRITE (6,2520) I,IP(I) 2520 FORMAT (1X,I3,14X,I3) 2530 CONTINUE IF (MTYPE.EQ.0) GO TO 2640 WRITE (6,2540) INPX1 2540 FORMAT (///,' Centroids for the ',I3,' cluster solution') IF (ICLUS1.EQ.ICLUS2) THEN OPEN (7,FILE='CENT.DAT',STATUS='UNKNOWN') ENDIF DO 2630 I = 1, INPX1 GM(I) = 0.0 INCT = 0 DO 2550 J = 1, NV P(J) = 0.0 2550 CONTINUE DO 2570 J = 1, N IF (IP(J).NE.I) GO TO 2570 INCT = INCT+1 DO 2560 L = 1, NV P(L) = P(L)+DATA(L,J) 2560 CONTINUE 2570 CONTINUE DO 2580 L = 1, NV GM(L) = P(L)/INCT 2580 CONTINUE WRITE (6,2590) I 2590 FORMAT (//,' Centroids for cluster ',I3,/,30('-'),/) DO 2610 J = 1, NV WRITE (6,2600) J,GM(J) 2600 FORMAT (1X,'Variable ',I3,1X,F10.3) 2610 CONTINUE IF (ICLUS1.EQ.ICLUS2) THEN WRITE (7,2620) (GM(J),J=1,NV) 2620 FORMAT (30G15.5) ENDIF 2630 CONTINUE C C Compute Internal Influence Measure C 2640 IF (INF.EQ.1) THEN WRITE (*,2650) INPX1 2650 FORMAT (/,' Computing internal influence measure for', * ' cluster level',I3,//) CALL CHENG (K,INPX1,N,NV,IP) ENDIF IF (IR.EQ.1.AND.INPX1.EQ.ICLUS1) THEN WRITE (*,2660) 2660 FORMAT (/,' Conducting test of cluster structure:',//) ENDIF C C Compute Single Partition Criterion Statistic C 2670 CONTINUE SUM1 = 0.0 SUM2 = 0.0 SUM3 = 0.0 SUM4 = 0.0 ASSIGN 2680 TO IVY IF (IR.EQ.1.AND.JR.EQ.0) ASSIGN 2690 TO IVY DO 2710 I = 2, N II = I-1 DO 2710 JJ = 1, II GO TO IVY ,(2680,2690) 2680 Z = E(I,JJ) GO TO 2700 2690 Z = R(I,JJ) 2700 SUM3 = SUM3+Z SUM4 = SUM4+Z*Z IF (WD(JJ,I)) GO TO 2710 SUM1 = SUM1+1. SUM2 = SUM2+Z 2710 CONTINUE TOP = SUM2-SUM1*SUM3/NIC1 BOTTOM = SQRT((SUM1-SUM1*SUM1/NIC1)*(SUM4-SUM3*SUM3/NIC1)) IF (IR.NE.1.OR.JR.EQ.1) GO TO 2720 RCORR(IIPX1,NR) = TOP/BOTTOM XR(IIPX1) = XR(IIPX1)+RCORR(IIPX1,NR) GO TO 2730 2720 CONTINUE PCORR(IIPX1) = TOP/BOTTOM*S PX(IIPX1) = XX(IIPX1) 2730 CONTINUE IF (IR.EQ.1.AND.JR.EQ.0) GO TO 2880 PCORR(N1) = 0.0 PX(N1) = XX(N1) C C Stepsize Computation C GS(1) = 0.0 DO 2740 MWG = 2, N1 MWG1 = MWG-1 GS(MWG) = XX(MWG)-XX(MWG1) 2740 CONTINUE IF (JR.EQ.1) GO TO 2870 C C Output Statistics C 2750 CONTINUE WRITE (6,2150) 12 IF (IR.EQ.1) THEN WRITE (6,2760) 2760 FORMAT (' Output Statistics:',//, * ' Number of Proximity Stepsize Point P-Value', * ' Simulated',/, * ' Clusters Value Biserial', * ' Mean Std. Dev.',/,75('-')) ELSE WRITE (6,2770) 2770 FORMAT (' Output Statistics:',//, * ' Number of Proximity Stepsize Point Biserial',/, * ' Clusters Value Goodness-of-Fit',/,50('-') * ) ENDIF DO 2830 IRC = 1, N2 IF (IRC.EQ.1) GO TO 2780 IF (GS(IRC).EQ.0.0) GO TO 2830 IF (GS(IRCXX).NE.0.0) GO TO 2800 2780 IRCXX = IRC+1 2790 IRCXX = IRCXX+1 IF (IRCXX.EQ.N) GO TO 2840 IF (GS(IRCXX).EQ.0.0) GO TO 2790 2800 JRC = N-IRC PX(IRC) = PX(IRC)*IS GS(IRC) = GS(IRC)*IS IF (IR.NE.1) THEN WRITE (6,2810) JRC,PX(IRC),GS(IRC),PCORR(IRC) 2810 FORMAT (2X,I4,4X,G10.3,2X,G9.3,6X,F6.3) ELSE WRITE (6,2820) JRC,PX(IRC),GS(IRC),PCORR(IRC),SP(IRC),XR(IRC * ),DR(IRC) 2820 FORMAT (2X,I4,5X,G10.3,1X,G9.3,3X,F6.3,3X,F9.4,4X,F6.3,3X, * F8.4) ENDIF 2830 CONTINUE 2840 JRC = 1 PX(IRC) = PX(IRC)*IS GS(IRC) = GS(IRC)*IS IF (IR.NE.1) THEN WRITE (6,2810) JRC,PX(IRC),GS(IRC),PCORR(IRC) ELSE WRITE (6,2850) JRC,PX(IRC),GS(IRC),PCORR(IRC) 2850 FORMAT (2X,I4,5X,G10.3,1X,G9.3,3X,F6.3) ENDIF C C More Analyses To Run? C IKOUNT = IKOUNT+1 WRITE (6,2860) 2860 FORMAT (///,' End of Method') WRITE (6,2150) 12 IF (I9(IKOUNT).EQ.0) GO TO 2910 K = I9(IKOUNT) GO TO 1680 2870 JR = 0 NR = 1 GO TO 1710 C C Compute Significance Level Of Internal Criterion Statistic C 2880 CONTINUE NR = NR+1 IF (NR.LE.NRMAX) GO TO 1710 DO 2900 I = 1, N2 SP(I) = 0.000001**2 XR(I) = XR(I)/NRMAX DO 2890 J = 1, NRMAX X = RCORR(I,J) IF (X.GT.PCORR(I)) SP(I) = SP(I)+1 DR(I) = DR(I)+(XR(I)-X)*(XR(I)-X) 2890 CONTINUE DR(I) = SQRT(DR(I)/NRMAX) SP(I) = SP(I)/NRMAX 2900 CONTINUE GO TO 2750 2910 CONTINUE CALL BEEP WRITE (*,2920) 2920 FORMAT (//,' Job successfully completed. ') WRITE (*,*) ' ' WRITE (*,*) 'Continue with additional analyses?' WRITE (*,*) ' 1 = Yes' WRITE (*,*) ' 0 = No (Default)' 2930 READ (*,50,ERR=2940) MORE IF (MORE.EQ.-99.OR.MORE.EQ.0) GO TO 2970 IF (MORE.NE.1) GO TO 2940 REWIND 5 GO TO 1070 2940 CALL BEEP WRITE (*,920) GO TO 2930 C C Reset User's Palette C 2950 CALL BEEP WRITE (*,*) 'Program has reached the end of your data file' WRITE (*,*) 'but your format and data set size specifications' WRITE (*,*) 'indicate that more data is expected or required.' WRITE (*,*) 'Program terminating...' GO TO 2970 2960 CALL BEEP WRITE (*,*) 'invalid, Inconsistent, or incomplete data format' WRITE (*,*) 'Program terminating...' 2970 CONTINUE WRITE (*,*) WRITE (*,*) 'Press enter to return to system...' WRITE (*,*) IKEY = PAUSE() IF (IARGV.NE.0) GO TO 2980 ERR = SET_VIDEO_MODE(OMODE) 2980 CALL TEXT_MODE STOP END C*********************************************************************** C C Subroutine To Compute Internal Measure Of Influence For Each Point C C*********************************************************************** C SUBROUTINE CHENG(METHOD,ICREQ,N,NV,IP) C PARAMETER(MAXPTS=200,MAXVAR=30) C C MAXPTS = Maximum number of data points C MAXVAR = Maximum number of variables C IMPLICIT INTEGER*4 (I-N) INTEGER*4 FLAG(MAXPTS) COMMON E(MAXPTS,MAXPTS),BETA,BETALW DIMENSION P(MAXPTS),ICLUSN(MAXPTS), * NIJH2(MAXPTS,(MAXPTS-1)),ITCCH2((MAXPTS-1)),SMIJH2(MAXPTS), * ICLUSV(MAXPTS),IGRUP((MAXPTS-1)),IP(MAXPTS),HA2(MAXPTS), * D(MAXPTS,MAXPTS),XX(MAXPTS) LOGICAL W(MAXPTS) WRITE (6,10) 10 FORMAT (///,' Internal Influence Measure',//,' Point',4X, * ' Measure',/,' -------------------') DO 20 I = 1, N DO 20 J = 1, N D(I,J) = E(I,J) 20 CONTINUE C C Initialize Iclusv Array (For Cluster Membership) C DO 30 I = 1, N ICLUSV(I) = 0 30 CONTINUE C C Delete One Point At A Time, Deleting Point #0 Means Using All C DO 460 ID = 0, N WRITE (*,40) ID 40 FORMAT ('+---> Computing influence for data point:',I4) C C Nsub Is The Adjustment Value For The Number Of Data Points Used C In Calculating Hatmh2 C IF (ID.EQ.0) THEN NSUB = 0 ELSE NSUB = 1 ENDIF C C Initialization For Criterion Statistics C DO 60 I = 1, N P(I) = 1 SMIJH2(I) = 1 DO 50 J = 1, ICREQ NIJH2(I,J) = 0 50 CONTINUE 60 CONTINUE C C Set Up N(I,J) Array C IF (ID.EQ.0) THEN DO 70 I = 1, N NIJH2(I,IP(I)) = 1 70 CONTINUE ELSE DO 90 I = 1, N DO 80 J = 1, ICREQ IF (ICLUSV(I).EQ.IGRUP(J)) THEN NIJH2(I,J) = 1 GO TO 90 ENDIF 80 CONTINUE 90 CONTINUE ENDIF C C Compute N(.,J) Information For External Criterion Statistic C SMCWH2 = 0. IF (ID.EQ.0) THEN DO 100 I = 1, ICREQ ITCCH2(I) = 0 100 CONTINUE DO 120 I = 1, ICREQ DO 110 J = 1, N ITCCH2(I) = ITCCH2(I)+NIJH2(J,I) 110 CONTINUE SMCWH2 = SMCWH2+ITCCH2(I)**2 120 CONTINUE ELSE C C Check The Cluster Id Of The Point Being Removed And Reduced The C Corresponding N(.,J) Value By 1. C DO 130 K = 1, ICREQ IF (ICLUSV(ID).EQ.IGRUP(K)) THEN ITCCH2(K) = ITCCH2(K)-1 C C Memorize Cluster Number So As To Reset ITCCH2 Array C IMEM = K ENDIF 130 CONTINUE DO 140 K = 1, ICREQ SMCWH2 = SMCWH2+ITCCH2(K)**2 140 CONTINUE ENDIF C C Initialize For Clustering Process C DO 160 I = 1, N ICLUSN(I) = I IF (I.EQ.ID) THEN W(I) = .TRUE. ELSE W(I) = .FALSE. DO 150 J = 1, N D(I,J) = E(I,J) 150 CONTINUE ENDIF 160 CONTINUE C C Begin Clustering Process - C Find Smallest Remaining Distance At Each Iteration C N2 = N-2 DO 430 II = 1, N2-NSUB DO 180 NI = 2, N IF (NI.NE.ID) THEN IF (W(NI)) GO TO 180 NI1 = NI-1 DO 170 NJ = 1, NI1 IF (.NOT.W(NJ)) GO TO 190 170 CONTINUE ENDIF 180 CONTINUE STOP 190 X = D(NI,NJ) DO 210 I = 2, N IF (I.NE.ID) THEN IF (W(I)) GO TO 210 I1 = I-1 DO 200 J = 1, I1 IF (W(J)) GO TO 200 IF ((D(I,J)-X).GT.0.0) GO TO 200 NI = I NJ = J X = D(I,J) 200 CONTINUE ENDIF 210 CONTINUE XX(II) = X W(NJ) = .TRUE. C C Update D Array According To Method C DO 320 L = 1, N IF (L.NE.ID) THEN IF (W(L)) GO TO 320 IF (L.EQ.NI) GO TO 320 C C Branch To Method C GO TO (220,230,240,250,260,270,280,290,300,470), * METHOD 220 IF ((D(NI,L)-D(NJ,L)).LE.0.0) GO TO 310 D(NI,L) = D(NJ,L) D(L,NI) = D(L,NJ) GO TO 310 230 IF ((D(NI,L)-D(NJ,L)).GE.0.0) GO TO 310 D(NI,L) = D(NJ,L) D(L,NI) = D(L,NJ) GO TO 310 240 D(L,NI) = (P(NI)/(P(NI)+P(NJ)))*D(L,NI)+(P(NJ)/(P(NI)+ * P(NJ)))*D(L,NJ) D(NI,L) = D(L,NI) GO TO 310 250 D(L,NI) = .5*D(L,NI)+.5*D(L,NJ) D(NI,L) = D(L,NI) GO TO 310 260 D(L,NI) = (P(NI)/(P(NI)+P(NJ)))*D(L,NI)+(P(NJ)/(P(NI)+ * P(NJ)))*D(L,NJ)-((P(NI)*P(NJ))/((P(NI)+P(NJ))**2.)) * *D(NI,NJ) D(NI,L) = D(L,NI) GO TO 310 270 D(L,NI) = .5*D(L,NI)+.5*D(L,NJ)-.25*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 310 280 D(L,NI) = ((P(NI)+P(L))/(P(NI)+P(NJ)+P(L)))*D(L,NI)+(( * P(NJ)+P(L))/(P(NI)+P(NJ)+P(L)))*D(L,NJ)-(P(L)/(P(NI * )+P(NJ)+P(L)))*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 310 290 D(L,NI) = .5*(1.-BETALW)*D(L,NI)+.5*(1.-BETALW)*D(L,NJ * )+BETALW*D(NI,NJ) D(NI,L) = D(L,NI) GO TO 310 300 D(L,NI) = (1.-BETA)*(P(NI)/(P(NI)+P(NJ)))*D(L,NI)+(1.- * BETA)*(P(NJ)/(P(NI)+P(NJ)))*D(L,NJ)+BETA*(D(NI,NJ)) D(NI,L) = D(L,NI) 310 CONTINUE ENDIF 320 CONTINUE P(NI) = P(NI)+P(NJ) C C Update N(I,J), N(I,.) Arrays For External Criteria C SMIJH2(NI) = 0 DO 330 II1 = 1, ICREQ NIJH2(NI,II1) = NIJH2(NJ,II1)+NIJH2(NI,II1) SMIJH2(NI) = SMIJH2(NI)+NIJH2(NI,II1)*NIJH2(NI,II1) 330 CONTINUE C C Set Up Flags To Indicate Status Of Each Point C SMSQH2 = 0. SUMRW = 0. SMRCH2 = 0. INI = ICLUSN(NI) INJ = ICLUSN(NJ) DO 380 II1 = 1, N IF (II1.NE.ID) THEN IF (ICLUSN(II1)-INI) 340, 350, 340 340 IF (ICLUSN(II1)-INJ) 370, 360, 370 350 FLAG(II1) = 1 GO TO 380 360 FLAG(II1) = -1 GO TO 380 370 FLAG(II1) = 0 ENDIF 380 CONTINUE C C Update Values Computed From N(I,J) Array C DO 400 II1 = 1, N IF (II1.NE.ID) THEN IF (ICLUSN(II1).EQ.INJ) ICLUSN(II1) = INI IF (W(II1)) GO TO 400 SMSQH2 = SMIJH2(II1)+SMSQH2 SUMRW = P(II1)*P(II1)+SUMRW DO 390 II2 = 1, ICREQ SMRCH2 = (P(II1)*P(II1))*(ITCCH2(II2)*ITCCH2(II2))+ * SMRCH2 390 CONTINUE ENDIF 400 CONTINUE C C Calculate Internal Influence Measure C HATMH2 = (1./(N-NSUB-1.))*(SMRCH2/(N-NSUB)-SMCWH2-SUMRW+(N- * NSUB)) HA2(II) = (SMSQH2-(N-NSUB)-HATMH2)/(SMCWH2/2.+SUMRW/2.-(N- * NSUB)-HATMH2) INPX1 = N-NSUB-II IF (INPX1.NE.ICREQ) GO TO 430 IF (ID.NE.0) THEN IF (ABS(HA2(II)-1.).LT.0.0001) THEN WRITE (6,410) ID,HA2(II) 410 FORMAT (1X,I3,4X,F10.5) ELSE WRITE (6,420) ID,HA2(II) 420 FORMAT (1X,I3,4X,F10.5,5X,'<-----') ENDIF ENDIF GO TO 440 430 CONTINUE 440 IF (ID.EQ.0) THEN J = 1 DO 450 I = 1, N ICLUSV(I) = ICLUSN(I) IF (.NOT.W(I)) THEN ITCCH2(J) = P(I) IGRUP(J) = ICLUSN(I) J = J+1 ENDIF 450 CONTINUE ENDIF C C Reset Itcch2 Array C IF (ID.NE.0) THEN ITCCH2(IMEM) = ITCCH2(IMEM)+1 ENDIF 460 CONTINUE C C End Of Method, Return To Main Program C 470 RETURN END C*********************************************************************** C C RND -- Decent RND, from Knuth, volume 2, page 88, 11/21/81, Dr. N. B. C C RND is a simple and relatively fast random number routine. T C algorithm used is one from Knuth, "The Art of Computing", volume 2, pa C 88. As Knuth shows, the random numbers generated are of good qualit C C The algorithm is as follows: take a large integer XN. Replace C by XN times 5**15 plus 1 modulo 2**25, i.e., C C XN <-- (XN * 5**15 + 1) mod 2**35 . C C The floating point random number returned by the function is, then, R C = XN / 2**35. C C The numbers returned on calls where XN is not changed appear to b C independent and uniformly distributed on the interval [0, 1]. C C XN is the random number seed. It should be specified at the C call of the routine and not changed if the same stream of numbers is C desired. C C RAND is an array of up to 200 random numbers as generated by the C program. C C NTOT is the total number of numbers to be generated. C C Most of the arithmetic required is arithmetic on large who C numbers. This arithmetic can be performed in DOUBLE PRECISION floati C point as long as the whole numbers are sufficiently small. Using t C DMOD function three times keeps the numbers sufficiently small. C C Written in Fortran 66 by Dr. N. B. Waite on November 21, 198 C Modified by Glenn W. Milligan on February 15, 1982. C*********************************************************************** C SUBROUTINE RANDOM(XN,NTOT,RND) C PARAMETER(MAXPTS=200,MAXVAR=30) C C MAXPTS = Maximum number of data points C MAXVAR = Maximum number of variables C DOUBLE PRECISION XN,C DIMENSION RND(MAXPTS*MAXVAR) DATA C / 34359738368.0D0 / DO 10 I = 1, NTOT XN = DMOD(DMOD(DMOD(XN*3125,C)*3125,C)*3125+1,C) RND(I) = XN/C 10 CONTINUE RETURN END c#ICICLE, FORTRAN VERSION. c#The authors of this software are Joseph B Kruskal and James M Landwehr C c#Copyright (c) 1993 by AT&T. c#Permission to use, copy, modify, and distribute this software for any c#purpose without fee is hereby granted, provided that this entire notic c#is included in all copies of any software which is or includes a copy c#or modification of this software and in all copies of the supporting c#documentation for such software. c#THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIE c#WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY c#REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY c#OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. c c-An explanation of ICICLE plots may be found in The American c-Statistician, May 1983, Volume 37, Number 2, pages 162-168. c c#This file contains ICICLE proper, a program for printing ICICLE plots, c#in Fortran. It consists of five subroutines (there is no main routine) c#It is the Fortran expansion of the ICICLE.r file, c#except that the comments have been modified as appropriate. c c-This file was created by Joseph B Kruskal, James M Landwehr, and Jean c-E McRae at AT&T Bell Laboratories, Murray Hill, New Jersey, U.S.A. c-It was completed in November, 1984. It has been thoroughly tested c-by the authors. c c-Copyright November, 1984 by AT&T Bell Laboratories. c c-AT&T Bell Laboratories makes no warranties, express or implied, with c-respect to this program. In particular, Bell Laboratories makes no c-warranty of merchantability, fitness for a particular use, freedom c-from infringement of any patent, copyright or trademark, nor any c-warranty as to accuracy. Bell Laboratories assumes no obligation to c-furnish any assistance of any kind whatsoever, or to furnish any c-additional information or documentaion. c c# Explanation of the comment convention used: c# A comment that starts with "c-" is comment line in Ratfor version. c# A comment that starts with "c+" comes from line with statement on it, c# and belongs to the immediately following Fortran line, if any. c# A comment that starts with "c " corresponds to a Ratfor statement. c# A comment that starts with "c#" was added or altered after expansion. c c- This program contains five subroutines, "icicle", "icilev", c-"icipos", "icidis", "icires", separated by lines of #-#-#-... . c- "icicle" actually prints out the ICICLE plot. c-The other 4 subroutines are optional routines for "mode-setting" or c-immediately returning extra information. Mode-setting modifies the c-display mode used in subsequent calls to "icicle", as indicated here: c- "icilev" can be used to change what clustering LEVels are displayed, c- or to return to the calling routine those which were displayed. c- "icipos" can be used to change the POSitions used for the objects, c- or to return to the calling routine those which were used. c- "icidis" can be used to change the DISplay appearance. c- "icires" can be used to RESet all modes to default values. c c--#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# c c# The Ratfor DEFINE statements and their explanations have been deleted c# because there is nothing in Fortran that corresonds to DEFINE. c# The page width and all array sizes in this Fortran version are those c# which result from the DEFINE statements in the Ratfor version. c c--#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# c c-subroutine icicle, Fortran version. c c-Produces the ICICLE plot c SUBROUTINE ICICLE( NOBJ, JOI1, JOI2, JOILEV, SDISWW, NLABCM, * NLABCH, EXTLAB, TITLE,LK) c-EXPLANATION OF THE ARGUMENTS, IN ORDER c- nobj = input, number of objects c- joi1(*) = input; "join 1" is first input array of join pointers c- joi2(*) = input; "join 2" is second input array of join pointers c- joilev(*) = input; "join level" has levels at which clusters join c- sdisww = input; similarity-dissimilarity switch c- +1.0=dissim, -1.0=simil c- nlabcm = input; "no. of label char.s: max possible" is c- number of rows in external label array c- nlabch = input; "no. of label char.s" is c- number of characters in longest external label c- extlab(nlabcm,*) = input; external label = user labels for objects c- title = input; user title for this run c c- EXPLANATION OF SOME OTHER VARIABLES c- width = width per position in display. Default = min = 2 c- distot = width of total display (may be several pages wide) c- nchar1 = 1 + maximum number of characters in a label c- labsw = label switch. c- 0 = use object numbers, else use external labels c- objll(MAXN) = left-most object of cluster at present join level c- which contains indicated object c- objr(MAXN) = right-hand neighbor of indicated object c- objjoi(MAXN) = right-most object of left cluster making join c- NCHARMAX1 = 1+ max number of characters in internal labels c- wall(NCHARMAX1,MAXN) = array showing which potential walls are in use c- label(MAXN) = internal array of labels for the objects c- displa(DISPLAYWIDTHTOTAL)= display = one-line buffer for output c c- diswid = width of display on one page. c- nstrip = number of strips = how many pages wide the display is c c-Explanation of logarithmic spacing of display levels c- DISRATIO is compilation-time variable, defined below as 200 c- Two constants, disadd and disrat , are calculated. c- The additive constant disadd is chosen so that c- (top display level + disadd)/(bottom display level + disadd)=DISRATIO c- The multiplicative constant disrat is chosen so that c- (display level + disadd) is uniformly spaced on log scale. c- Thus each new display level is calculated by c- (display level + disadd)/(next display level + disadd) = disrat c PARAMETER(MAXPTS=200,MAXVAR=30) c c MAXPTS = Maximum number of data points c MAXVAR = Maximum number of variables c c-Common shared with the subroutines icilev, icipos, icidis c COMMON /ICLEV/ KDISLE,NDISLE,DISLEV COMMON /ICPOS/ KOBJPO,OBJPOS COMMON /ICDIS/ WIDTH, ECHAR, FCHAR, JCHAR, BCHAR c-Nonarray variables c INTEGER NOBJ, KDISLE, NDISLE, NDISL1, KOBJPO, NLABCM, NLABCH, * NCHAR INTEGER LABSW, WIDTH, DISTOT, DIGIT, LDIS, RDIS, NCHAR1 REAL SDISSW, DISLE, SDISWW, SDISLE CHARACTER ECHAR*1, JCHAR*1, FCHAR*1, BCHAR*1, SCHAR*1 c-Array variables which are arguments c INTEGER JOI1(MAXPTS), JOI2(MAXPTS), OBJPOS(MAXPTS*2), LK(MAXPTS) REAL JOILEV(MAXPTS) CHARACTER EXTLAB(6, MAXPTS)*1 DIMENSION TITLE(20) c-Array variables which are NOT arguments c INTEGER OBJLL(MAXPTS*2), OBJR(MAXPTS*2), OBJJOI(MAXPTS*2) REAL DISLEV(MAXPTS*2) CHARACTER LABEL(4,MAXPTS*2)*1, WALL(MAXPTS*2)*1 CHARACTER DISPLA(MAXPTS*6)*1 c-Default values for the modes c+switch which controls use of levels c DATA KDISLE /0/ c c+switch which controls positioning of objects c DATA KOBJPO /0/ c c- echar, fchar, jchar, bchar are the display characters c DATA WIDTH, ECHAR, FCHAR, JCHAR, BCHAR /2, "&", "&", "=", " "/ c-EXECUTABLE CODE STARTS HERE c c-PART 1. INITIALIZATION c+distot = display width total (all strips) c WRITE (*,10) 10 FORMAT (/,' Generating icicle plot, please wait...') WRITE (6,20) 20 FORMAT (/,' If the entire Icicle plot does not fit on a single', * ' page,',/,' please place the pages side by side from', * ' left to right.') SDISSW = 1.0 SSTORE = 1.0 IF (SDISWW.LT.0) SSTORE = -1.0 ICAUGH = 0 DO 30 IFIXX = 1, NOBJ IP1 = LK(IFIXX)/10 EXTLAB(1,IFIXX) = CHAR(IP1+48) IP2 = LK(IFIXX)-10*IP1 EXTLAB(2,IFIXX) = CHAR(IP2+48) 30 CONTINUE IF (NOBJ.LT.100) GO TO 50 DO 40 IFIXX = 1, NOBJ IP0 = LK(IFIXX)/100 EXTLAB(1,IFIXX) = CHAR(IP0+48) IP1 = (LK(IFIXX)-100*IP0)/10 EXTLAB(2,IFIXX) = CHAR(IP1+48) IP2 = LK(IFIXX)-100*IP0-10*IP1 EXTLAB(3,IFIXX) = CHAR(IP2+48) 40 CONTINUE 50 DISTOT = WIDTH*(NOBJ-1)+1 c+fnobj = floating number of objects c FNOBJ = FLOAT(NOBJ) c c+nobjdi = number of object digits c NOBJDI = 1+IFIX(ALOG10(FNOBJ+0.5)) IF (.NOT.(NLABCH.NE.0)) GO TO 60 NCHAR = NLABCH LABSW = 1 c c+nchar = number of characters c GO TO 70 c c else c 60 CONTINUE NCHAR = NOBJDI LABSW = 0 70 CONTINUE c c+labsw = label switch: 1=ext, 0=int c NCHAR1 = NCHAR+1 c+npos = number of object positions in display c NPOS = NOBJ c c- npos is always equal to nobj. Distinction is for conceptual clarity. c+njoile = number of join levels in display c NJOILE = NOBJ-1 c-Initialize handling of display levels. c switch c I23002 = (KDISLE) IF (.NOT.(I23002.EQ.(2).OR.I23002.EQ.(3))) GO TO 80 c c- Use join levels as display levels. c- Therefore, ignore dislev. Instead, use one level for each join c NDISL1 = NJOILE GO TO 130 80 CONTINUE IF (.NOT.(I23002.EQ.(1))) GO TO 90 c c- Use display levels provided in dislev c- Use number of display levels from ndisle c NDISL1 = NDISLE GO TO 130 90 CONTINUE IF (.NOT.(I23002.EQ.(0))) GO TO 120 c c- Logarithmic spacing of display levels. c IF (.NOT.(NDISLE.NE.0)) GO TO 100 NDISL1 = NDISLE GO TO 110 c c else c 100 CONTINUE NDISL1 = NJOILE 110 CONTINUE c c- Calculate disrat and disadd, based on rationale explained above c DISTOP = JOILEV(NJOILE) DISBOT = JOILEV(1) DISADD = (DISTOP-200.0*DISBOT)/(200.0-1.0) DISRAT = (1.0/200.0)**(1.0/FLOAT(NDISL1-1)) 120 CONTINUE 130 CONTINUE c c+For use in icilev c NDISLE = NDISL1 c c c-PART 2. CALCULATE objll(), objr() RECURSIVELY. ALSO CALC objjoi() c-Recursion follows the original merging process c-At each stage of recursion, for clusters defined at that stage: c- objll yields left-most object of largest current cluster containing c- argument object. c- objr yields right-hand neighbor of argument object in the largest c- current cluster containing argument object. However, c- if argument is right-most object, objr yields object itself c- objjoi yields right-most object of left cluster of argument level. c c-Initialize objll(), objr(). Each initial cluster has 1 object. c DO 140 IOBJ = 1, NOBJ OBJLL(IOBJ) = IOBJ OBJR(IOBJ) = IOBJ 140 CONTINUE c-Recursive calculation itself c DO 410 IJOILE = 1, NJOILE c c+Start of long do loop on ijoile c c+An object in first cluster being joined c I1OBJ = JOI1(IJOILE) c c+Find left-most object in first cluster c I1OBJL = OBJLL(I1OBJ) c c-Find right-most object in first cluster c+right-hand neighbor of object WITHIN cluster c I1OBJR = OBJR(I1OBJ) c c while c 150 IF (.NOT.(I1OBJR.NE.OBJR(I1OBJR))) GO TO 170 c c+Iterate till right end of cl reached c I1OBJR = OBJR(I1OBJR) IF (.NOT.(.NOT.(1.LE.I1OBJR.AND.I1OBJR.LE.NOBJ))) GO TO 160 c c- i1objr should lie between 1 and nobj c NERROR = 2 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 2.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 160 CONTINUE GO TO 150 c c endwhile c 170 CONTINUE c c-(Still in long do loop on ijoile) c+An object in second cluster being joined c I2OBJ = JOI2(IJOILE) c c+Find left-most object in second cluster c I2OBJL = OBJLL(I2OBJ) c c-Find right-most object in second cluster c+right-hand neighbor of object WITHIN cluster c I2OBJR = OBJR(I2OBJ) c c while c 180 IF (.NOT.(I2OBJR.NE.OBJR(I2OBJR))) GO TO 200 c c+Iterate until right end of cl reached c I2OBJR = OBJR(I2OBJR) IF (.NOT.(.NOT.(1.LE.I2OBJR.AND.I2OBJR.LE.NOBJ))) GO TO 190 c c- i2objr should lie between 1 and nobj c NERROR = 3 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 3.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 190 CONTINUE GO TO 180 c c endwhile c 200 CONTINUE c-(Still in long do loop on ijoile) c-Decide in which order to join the two clusters c- joisid=1 means first cluster on left c- joisid=2 means second cluster on left c-Present version of program: c- if using object position arrangement provided by calling routine c- (kobjpo=1), chooses case as required c- if creating object position arrangement (kobjpo=0) c- always uses case 1 c-Possible improvement of program: when creating arrangement (kobjpo=0), c-make choice of case depend on cluster sizes or other things, so as c-to obtain desirable ICICLE diagram. c c switch c I23020 = (KOBJPO) IF (.NOT.(I23020.EQ.(1))) GO TO 270 c c+Use arrangement given in objpos c JOISID = 0 DO 240 IPOS = 1, NPOS IF (.NOT.(OBJPOS(IPOS).EQ.I1OBJL)) GO TO 210 c c+first cluster on left c JOISID = 1 210 CONTINUE IF (.NOT.(OBJPOS(IPOS).EQ.I2OBJL)) GO TO 220 c c+second cluster on left c JOISID = 2 220 CONTINUE IF (.NOT.(JOISID.NE.0)) GO TO 230 GO TO 250 230 CONTINUE 240 CONTINUE 250 CONTINUE IF (.NOT.(JOISID.EQ.0)) GO TO 260 c c- joisid should be defined in preceding do loop c NERROR = 4 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 4.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 260 CONTINUE GO TO 290 270 CONTINUE IF (.NOT.(I23020.EQ.(0))) GO TO 280 c c+Create arrangement of objects if kobjpo=0 c+first cluster on left c JOISID = 1 280 CONTINUE 290 CONTINUE c c-(Still in long do loop on ijoile) c-Join the two clusters in order selected c switch c I23033 = (JOISID) c c+Start of switch on joisid c IF (.NOT.(I23033.EQ.(1))) GO TO 340 c c+first cluster on left c+Right-most object of left cluster c OBJJOI(IJOILE) = I1OBJR c c+Recursive updating of objr c OBJR(I1OBJR) = I2OBJL c c- Recursive updating of objll c IOBJ = I2OBJL c c repeat c 300 CONTINUE OBJLL(IOBJ) = I1OBJL IF (.NOT.(IOBJ.EQ.OBJR(IOBJ))) GO TO 310 GO TO 330 c c else c 310 CONTINUE IOBJ = OBJR(IOBJ) IF (.NOT.(.NOT.(1.LE.IOBJ.AND.IOBJ.LE.NOBJ))) GO TO 320 c c- iobj should lie between 1 and nobj c NERROR = 5 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 5.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 320 CONTINUE GO TO 300 330 CONTINUE GO TO 400 340 CONTINUE IF (.NOT.(I23033.EQ.(2))) GO TO 390 c c+second cluster on left c+Right-most object of left cluster c OBJJOI(IJOILE) = I2OBJR c c+Recursive updating of objr c OBJR(I2OBJR) = I1OBJL c c- Recursive updating of objll c IOBJ = I1OBJL c c repeat c 350 CONTINUE OBJLL(IOBJ) = I2OBJL IF (.NOT.(IOBJ.EQ.OBJR(IOBJ))) GO TO 360 GO TO 380 c c else c 360 CONTINUE IOBJ = OBJR(IOBJ) IF (.NOT.(.NOT.(1.LE.IOBJ.AND.IOBJ.LE.NOBJ))) GO TO 370 c c- iobj should lie between 1 and nobj c NERROR = 6 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 6.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 370 CONTINUE GO TO 350 380 CONTINUE 390 CONTINUE 400 CONTINUE c c+End of switch on joisid c 410 CONTINUE c c+End of long do loop on ijoile c-End of definition of objll, objr, objjoi c c+left-most object in entire set c LLLOBJ = OBJLL(1) c c c-PART 3. DEAL WITH OBJECT POSITION ARRAY, objpos() c switch c I23050 = (KOBJPO) c c+Start of switch on kojbpo c IF (.NOT.(I23050.EQ.(1))) GO TO 460 c c- If object position array was provided by calling program, check c- consistency with clustering provided by the calling program. c+To start with, assume consistent c KCONS = 0 c c for c IPOS = 1 420 IF (.NOT.(IPOS.LT.NPOS)) GO TO 440 IF (.NOT.(OBJR(OBJPOS(IPOS)).NE.OBJPOS(IPOS+1))) GO TO 430 KCONS = 1 GO TO 440 430 CONTINUE IPOS = IPOS+1 GO TO 420 c c endfor c 440 CONTINUE IF (.NOT.(KCONS.EQ.1)) GO TO 450 WRITE (6,*) '----------------------------------------' WRITE (6,*) 'ICICLE: GIVEN OBJECT POSITIONS NOT CONSISTENT WITH GI *VEN TREE,' WRITE (6,*) * 'SO THE ICICLE PLOT WILL USE MODIFIED OBJECT POSITIONS.' WRITE (6,*) '----------------------------------------' 450 CONTINUE GO TO 500 460 CONTINUE IF (.NOT.(I23050.EQ.(0))) GO TO 490 c c- Save object position c IOBJ = LLLOBJ c c for c IPOS = 1 470 IF (.NOT.(IPOS.LE.NPOS)) GO TO 480 OBJPOS(IPOS) = IOBJ IOBJ = OBJR(IOBJ) IPOS = IPOS+1 GO TO 470 c c endfor c 480 CONTINUE 490 CONTINUE 500 CONTINUE c c+End of switch on kobjpo c c-PART 4. PRINT ICICLE PLOT c-In general, display may be too wide to get on one page. c-It is broken up into a number of vertical strips, each one page wide. c-Entire display consists of several strips, to be placed side by side. c IF (.NOT.(KDISLE.NE.2)) GO TO 510 c c+diswid = display width which shows on one page c DISWID = 70 GO TO 520 c c else c 510 CONTINUE DISWID = 59 520 CONTINUE c c+nstrip = number of strips c NSTRIP = (DISTOT+DISWID-1)/DISWID c-Loop through number of strips necessary to show whole display. c-During loop c- istrip is number of the strip c- ldis is object position of left edge of strip c- rdis is object position of right edge of strip c-When generating display line for one strip, create display line for c-entire width, but only use it from ldis to rdis c ISTRIP = 0 c c for c LDIS = 1 530 IF (.NOT.(LDIS.LE.DISTOT)) GO TO 1140 c c+Start of long for loop through strips. Index is ldis. c LDTEST = LDIS+DISWID-1 RDIS = MIN(DISTOT,LDTEST) ISTRIP = ISTRIP+1 c-print header information at top of display c WRITE (6,540) 12 540 FORMAT (A1) IF (.NOT.(DISWID.LT.DISTOT)) GO TO 550 WRITE (6,560) ISTRIP,NSTRIP 550 CONTINUE 560 FORMAT ('Tree output strip',I3,' of',I3,' strips',/) c-initialize the label array, label(,) c IOBJ = LLLOBJ DO 620 IPOS = 1, NPOS c c for c ICHAR = NCHAR 570 IF (.NOT.(ICHAR.GE.1)) GO TO 600 I = IOBJ IF (.NOT.(LABSW.NE.0)) GO TO 580 LABEL(ICHAR,IPOS) = EXTLAB(ICHAR,IOBJ) GO TO 590 c c else c 580 CONTINUE DIGIT = MOD(I,10) I = (I-DIGIT)/10 c c write (*char,1) digit c1 format(i1) c LABEL(ICHAR,IPOS) = SCHAR 590 CONTINUE ICHAR = ICHAR-1 GO TO 570 c c endfor c 600 CONTINUE LABEL(NCHAR1,IPOS) = ECHAR IOBJ = OBJR(IOBJ) IF (.NOT.(.NOT.(1.LE.IOBJ.AND.IOBJ.LE.NOBJ))) GO TO 610 c c- iobj should be between 1 and nobj c NERROR = 7 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 7.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 610 CONTINUE 620 CONTINUE c+print object labels at top of columns c DO 760 ICHAR = 1, NCHAR DO 660 IPOS = 1, NPOS IDIS1 = WIDTH*(IPOS-1)+1 DISPLA(IDIS1) = LABEL(ICHAR,IPOS) IF (.NOT.(IPOS.EQ.NPOS)) GO TO 630 GO TO 670 630 CONTINUE c c for c IDIS = IDIS1+1 640 IF (.NOT.(IDIS.LE.WIDTH*IPOS)) GO TO 650 DISPLA(IDIS) = ' ' IDIS = IDIS+1 GO TO 640 c c endfor c 650 CONTINUE 660 CONTINUE 670 CONTINUE IF (.NOT.(KDISLE.NE.2)) GO TO 680 WRITE (6,710) (DISPLA(IDIS),IDIS=LDIS,RDIS) GO TO 700 c c else c 680 CONTINUE IF (ICAUGH.NE.0) GO TO 690 WRITE (6,720) (DISPLA(IDIS),IDIS=LDIS,RDIS) ICAUGH = ICAUGH+1 GO TO 700 690 WRITE (6,730) (DISPLA(IDIS),IDIS=LDIS,RDIS) ICAUGH = ICAUGH+1 IF (NOBJ.LT.100) ICAUGH = 0 IF (ICAUGH.EQ.3) ICAUGH = 0 700 CONTINUE 710 FORMAT (7X,114A1) 720 FORMAT ('Item Number:',6X,103A1) 730 FORMAT (18X,103A1) c- change blanks to fill characters, in case some labels shorter c DO 750 IPOS = 1, NPOS IF (.NOT.(LABEL(ICHAR,IPOS).EQ.' ')) GO TO 740 LABEL(ICHAR,IPOS) = FCHAR 740 CONTINUE 750 CONTINUE 760 CONTINUE c c+blank line c WRITE (6,770) 770 FORMAT ('Level Proximity',/,' Value') c-PRINT BODY OF DISPLAY c c-Initialize join level counter, ijoile, and wall. c IJOILE = NJOILE DO 780 IPOS = 1, NPOS c c+fill wall with join character c WALL(IPOS) = JCHAR 780 CONTINUE c c+blank out last entry of wall c WALL(NPOS) = BCHAR c-Initialize phase counter, ipha, and set number of phases, npha. c IPHA = 0 NPHA = NCHAR1 c-Initialize display level counter, idisle, and disl23. c IDISLE = NDISL1 c c+Special indicator for cases 2 and 3 c DISL23 = 0.0 c-Initialize display level, disle. c switch c I23096 = (KDISLE) IF (.NOT.(I23096.EQ.(0))) GO TO 790 c c+Start with root node join level c DISLE = JOILEV(NJOILE) GO TO 820 790 CONTINUE IF (.NOT.(I23096.EQ.(1))) GO TO 800 c c+Use display levels from calling program c DISLE = DISLEV(NDISL1) GO TO 820 800 CONTINUE IF (.NOT.(I23096.EQ.(2).OR.I23096.EQ.(3))) GO TO 810 c c+Use join levels c DISLE = JOILEV(NDISL1) 810 CONTINUE 820 CONTINUE c c+Save disle in dislev c DISLEV(NDISL1) = DISLE c repeat c 830 CONTINUE c c+Start of repeat loop to print one strip. c c- Update ijoile, wall. c while c 840 IF (.NOT.((JOILEV(IJOILE)-DISLE)*SDISSW+DISL23.GT.0..AND. * IJOILE.GE.1)) GO TO 940 c c- Normally: Continue until joilev moves past disle. disl23=0.0 c- Cases 2,3, on first entry: Skip loop. disl23=0.0 c- Cases 2,3, later entry: Do loop just once. disl23=1.0 c+Start of while loop to update ijoile and wall. c IF (.NOT.(KDISLE.EQ.2.OR.KDISLE.EQ.3)) GO TO 850 c c+Reset 2,3 indicator. c DISL23 = 1.0 c c- Find right-most object of left cluster making join c 850 CONTINUE IOBJ = OBJJOI(IJOILE) c c- Find its position by counting from left-most object c IPOS = 1 c c for c JOBJ = LLLOBJ 860 IF (.NOT.(JOBJ.NE.IOBJ)) GO TO 880 IPOS = IPOS+1 IF (.NOT.(.NOT.(1.LE.JOBJ.AND.JOBJ.LE.NOBJ))) GO TO 870 c c- jobj should be between 1 and nobj c NERROR = 8 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: REST OF OUTPUT ABORTED, NERROR = 8.' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' RETURN 870 CONTINUE JOBJ = OBJR(JOBJ) GO TO 860 c c endfor c 880 CONTINUE c- Blank out the join character making this join c WALL(IPOS) = BCHAR c- Blank out new singleton groups. c IF (.NOT.(WALL(IPOS+1).EQ.BCHAR)) GO TO 900 DO 890 ICHAR = 1, NCHAR1 LABEL(ICHAR,IPOS+1) = BCHAR 890 CONTINUE 900 CONTINUE IF (.NOT.((IPOS.GT.1.AND.WALL(IPOS-1).EQ.BCHAR).OR.IPOS.EQ.1)) GO * TO 920 DO 910 ICHAR = 1, NCHAR1 LABEL(ICHAR,IPOS) = BCHAR 910 CONTINUE c- Update ijoile. c 920 CONTINUE IJOILE = IJOILE-1 IF (.NOT.(DISL23.EQ.1.0)) GO TO 930 c c+Move just one join level in cases 2,3 c GO TO 940 930 CONTINUE GO TO 840 c c endwhile c 940 CONTINUE c c+End of while loop to update ijoile and wall. c c- Update ipha. c IPHA = IPHA+1 IF (.NOT.(IPHA.GT.NPHA)) GO TO 950 IPHA = 1 c- Fill displa from label and wall. c 950 CONTINUE DO 990 IPOS = 1, NPOS IDIS1 = WIDTH*(IPOS-1)+1 DISPLA(IDIS1) = LABEL(IPHA,IPOS) IF (.NOT.(IPOS.EQ.NPOS)) GO TO 960 GO TO 1000 960 CONTINUE c c for c IDIS = IDIS1+1 970 IF (.NOT.(IDIS.LE.WIDTH*IPOS)) GO TO 980 DISPLA(IDIS) = WALL(IPOS) IDIS = IDIS+1 GO TO 970 c c endfor c 980 CONTINUE 990 CONTINUE 1000 CONTINUE c- THE PRIMARY OUTPUT STATEMENT c- Print one line of display. c IF (.NOT.(KDISLE.NE.2)) GO TO 1010 IDFIX = NOBJ-IDISLE IF (SSTORE.LT.0) THEN SDISLE = DISLE*SSTORE ELSE SDISLE = DISLE ENDIF WRITE (6,1050) IDFIX,SDISLE,(DISPLA(IDIS),IDIS=LDIS,RDIS) GO TO 1040 c c else c 1010 CONTINUE IF (.NOT.(IJOILE.EQ.NJOILE)) GO TO 1020 IDFIX = NOBJ-IDISLE IF (SSTORE.LT.0) THEN SDISLE = DISLE*SSTORE ELSE SDISLE = DISLE ENDIF WRITE (6,1060) IDFIX,SDISLE,(DISPLA(IDIS),IDIS=LDIS,RDIS) GO TO 1030 c c else c 1020 CONTINUE IDFIX = NOBJ-IDISLE IF (SSTORE.LT.0) THEN SDISLE = DISLE*SSTORE ELSE SDISLE = DISLE ENDIF WRITE (6,1070) IDFIX,SDISLE,(DISPLA(IDIS),IDIS=LDIS,RDIS) 1030 CONTINUE 1040 CONTINUE 1050 FORMAT (1X,I3,2X,G11.5,1X,114A1) 1060 FORMAT (1X,I3,2X,G11.5,1X,103A1) 1070 FORMAT (1X,I3,2X,G11.5,1X,103A1) c- update idisle and apply stopping rule c IDISLE = IDISLE-1 IF (.NOT.(IDISLE.LT.1)) GO TO 1080 GO TO 1130 c- Update disle and save disle in dislev c 1080 CONTINUE c c switch c I23137 = (KDISLE) IF (.NOT.(I23137.EQ.(0))) GO TO 1090 c c+Calculate display level c DISLE = (DISLE+DISADD)*DISRAT-DISADD GO TO 1120 1090 CONTINUE IF (.NOT.(I23137.EQ.(1))) GO TO 1100 c c+Use display level from calling program c DISLE = DISLEV(IDISLE) GO TO 1120 1100 CONTINUE IF (.NOT.(I23137.EQ.(2).OR.I23137.EQ.(3))) GO TO 1110 c c+Use next join level as display level c DISLE = JOILEV(IJOILE-1) 1110 CONTINUE 1120 CONTINUE c c+Save disle c DISLEV(IDISLE) = DISLE c+End of repeat loop to print one strip. c GO TO 830 1130 CONTINUE LDIS = LDIS+DISWID GO TO 530 c c endfor c 1140 CONTINUE c c+End of long for loop through strips. Index is ldis. c RETURN END c--#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# c c-subroutine icilev, Fortran version. c c-sets the display levels according to user preference c SUBROUTINE ICILEV(SW, ARRAY, N) c-Argument description follows: c- sw = input. It controls handling of levels, as described below. c- It can be -1, 0 (default value) 1, 2, 3, or 4. c- It changes the display mode, except when sw=-1. c- array = input or output. It holds display levels. c- It is used only with sw values -1, 1. c- n = input. It indicates the number of display levels c- It is used only with sw values -1, 0, 1. c c- What icilev does for different sw values. c--1 = Does NOT change display mode. Instead, has immediate effect: c- display level values are returned in "array" ("array" is output) c- If n!=0, n lowest levels returned; if n=0, all levels returned c- 0 = Log spacing used for display levels. c- If n!=0, n levels used; if n=0, no. of levels=number of objects. c- 1 = n levels from "array" are used in display("array" is input). c- 2 = One level is used for each join. Also, TWO display levels are c- printed on each line, to show INTERVAL over which this c- clustering is relevant. c- 3 = Same, but as space saver only FIRST display level is printed. c PARAMETER(MAXPTS=200,MAXVAR=30) c c MAXPTS = Maximum number of data points c MAXVAR = Maximum number of variables c COMMON /ICLEV/ KDISLE, NDISLE, DISLEV INTEGER KDISLE, N, N1, NDISLE, SW REAL DISLEV(MAXPTS*2), ARRAY(N) c switch c I23141 = (SW) IF (.NOT.(I23141.EQ.(-1).OR.I23141.EQ.(1))) GO TO 30 IF (.NOT.(N.EQ.0)) GO TO 10 N1 = NDISLE N = NDISLE GO TO 20 c c else c 10 CONTINUE N1 = N 20 CONTINUE 30 CONTINUE c switch c I23145 = (SW) IF (.NOT.(I23145.EQ.(-1))) GO TO 50 DO 40 I = 1, N1 ARRAY(I) = DISLEV(I) 40 CONTINUE RETURN 50 CONTINUE IF (.NOT.(I23145.EQ.(1))) GO TO 70 DO 60 I = 1, N DISLEV(I) = ARRAY(I) 60 CONTINUE GO TO 90 70 CONTINUE IF (.NOT.(I23145.EQ.(0).OR.I23145.EQ.(2).OR.I23145.EQ.(3))) GO TO * 80 GO TO 90 80 CONTINUE c c default c+kdisle has invalid value c WRITE (6,*) '----------------------------------------' WRITE (6,*) 'ICICLE: CALL TO ICILEV IGNORED, BECAUSE 1ST ARGUMENT *IS INVALID' WRITE (6,*) '----------------------------------------' RETURN 90 CONTINUE c c+levels switch c KDISLE = SW c c+number of levels c NDISLE = N RETURN END c--#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# c c-subroutine icipos, Fortran version. c c- Allows the user to set the position order of the objects in the c- display (i.e., change the display mode) in subsequent calls to c- icicle, or to find what position order was used during the last call c- to icicle c SUBROUTINE ICIPOS(SW, ARRAY, N) c- Argument description follows: c- sw is input; controls handling of object positions (see below). c- array is input or output; gives positions of objects. c- It is used only with sw values -1 and 1. c- n is input; number of positions in "array". c- It is used only with sw values -1 and 1. c c- sw has three possible values (0 is default value): c- -1 : Does NOT change display mode. Instead, has immediate effect: c- ICICLE returns n object positions (from previous use) c- in "array" ("array" is output). c- 0 : ICICLE will supply its own object positions in future c- 1 : ICICLE reads in n object positions from "array" (input). It c- will use these object positions in future calls to icicle. c- However, if they are not consistent with hierarchical c- clustering, modifications are made as necessary (but the input c- "array" is never changed). c PARAMETER(MAXPTS=200,MAXVAR=30) c c MAXPTS = Maximum number of data points c MAXVAR = Maximum number of variables c COMMON /ICPOS/ KOBJPO, OBJPOS INTEGER KOBJPO,OBJPOS(MAXPTS*2),SW,ARRAY(N),N c switch c I23153 = (SW) IF (.NOT.(I23153.EQ.(-1))) GO TO 20 DO 10 I = 1, N ARRAY(I) = OBJPOS(I) 10 CONTINUE RETURN 20 CONTINUE IF (.NOT.(I23153.EQ.(1))) GO TO 40 DO 30 I = 1, N OBJPOS(I) = ARRAY(I) 30 CONTINUE GO TO 60 40 CONTINUE IF (.NOT.(I23153.EQ.(0))) GO TO 50 GO TO 60 50 CONTINUE c c default c+invalid value c WRITE (6,*) '----------------------------------------' WRITE (6,*) 'ICICLE: CALL TO ICIPOS IGNORED, BECAUSE 1ST ARGUMENT *IS INVALID' WRITE (6,*) '----------------------------------------' RETURN 60 CONTINUE KOBJPO = SW RETURN END c--#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# c c-subroutine icidis, Fortran version. c c-Enters ICICLE display parameters c SUBROUTINE ICIDIS(WIDTHI, ECHARI, FCHARI, JCHARI, BCHARI) c-If this subroutine is not called, ICICLE will use default values of c- display parameters. c-width = width per position = number of columns from one object column c-to next in display, measured "center to center". minimum=default=2. c-echar = end character. Default = "&". c- Separates vertical repetitions of labels. c-fchar = fill character. Default = "&". c- Fills out shorter labels to length of longer ones. c-jchar = join character. Default = "=". c- Joins objects which belong to the same cluster. c-bchar = blank character. Default = " ". c- Used for white space which separates objects in c- distinct clusters. c c+display parameters c COMMON /ICDIS/ WIDTH, ECHAR, FCHAR, JCHAR, BCHAR INTEGER WIDTH, WIDTHI CHARACTER ECHARI*1, JCHARI*1, FCHARI*1, BCHARI*1 CHARACTER ECHAR*1, JCHAR*1, FCHAR*1, BCHAR*1 IF (.NOT.(.NOT.(2.LE.WIDTHI.AND.WIDTHI.LE.5))) GO TO 10 WRITE (6,*) '----------------------------------------' WRITE (6,*) 'ICICLE: WIDTH PER POSITION HAS BEEN CHANGED TO BETWEE *N 2 AND 5' WRITE (6,*) '----------------------------------------' 10 CONTINUE WIDTH = MIN(MAX(2,WIDTHI),5) IF (.NOT.(WIDTHI.GT.3)) GO TO 20 WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' WRITE (6,*) 'ICICLE: IF WIDTH PER POSITION > 3, ARRAY ''DISPLA'' M *AY OVERFLOW' WRITE (6,*) '%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+%+' 20 CONTINUE ECHAR = ECHARI JCHAR = JCHARI FCHAR = FCHARI BCHAR = BCHARI RETURN END c--#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-# c c-subroutine icires, Fortran version. c c-Resets all modes to their default values. c SUBROUTINE ICIRES INTEGER DUMMYI REAL DUMMYR(1) CALL ICILEV (0,DUMMYR,0) CALL ICIPOS (0,DUMMYR,DUMMYI) c c+default values for all arguments c CALL ICIDIS (2,'&','&','=',' ') RETURN END