C THE INPUT DECK IS OF THE WILLIAMS TYPE
C GLOSSARY
C AD IS AN ARRAY OF ALL FLOATING POINT DATA
C ID IS AN ARRAY OF ALL FIXED POINT DATA.
C TITLE IS A HEADING FOR RECOGNITION OF THE JOB.
C IBUF IS AN INTEGER BUFFER ON INPUT.
C BUF IS THE CORRESPONDING FLOATING POINT BUFFER.
C ISIZE RECORDS THE NUMBER OF STATES POSSIBLE FOR EACH MULTISTATE
C ATTRIBUTE.
C ICON IS THE CONTINGENCY TABLE.
C NV IS AN ARRAY RECORDING THE NUMBER OF VALUES IN EACH CLASS OF A
C CLASSIFICATION.
C X IS THE ARRAY OF VALUES IN A CLASSIFICATION OR IN A CORRELATION.
C Y IS THE SECOND ARRAY IN A CORRELATION.
C NENT IS THE NUMBER OF ENTITIES.
C NB IS THE NUMBER OF BINARY VALUES.
C NN IS THE NUNBER OF NUMERIC VALUES.
C ND IS THE NUMBER OF DISORDERED MULTISTATE ATTRIBUTES.
C NO IS THE NUMBER OF ORDERED MULTISTATE ATTRIBUTES.
C NM IS THE NUMBER OF MULTISTATE VARIABLES.
C I,J,AND K ARE USED AS DO LOOP SUBSCRIPTS.
C IRUN,JRUN,AND KRUN ARE COMPUTED VALUES FOR ACCESSING AD AND ID.
C DATA DECK HAS LAST 8 SPACES ON EACH RESERVED FOR IDENTIFICATION
C
C **************************************************************
C DIMENSION ID(NENT*(NB+ND)),AD(NENT*(NN+NO)),TITLE(10),IBUF(NB)
C DIMENSION BUF(NN),ISIZE(ND+NO)
C DIMENSION RANGE(NM),IN(NENT),S(NENT),LZ(NENT),DMIN(NENT)
C DIMENSION MINCOL(NENT),LAB(NENT)
C DIMENSION EDSQ(NENT/2*(NENT-1)),ED(NENT)
C **************************************************************
C **************************************************************
DIMENSION ID(4720),AD(2000),TITLE(13),IBUF(50),BUF(10),ISIZE(40)
DIMENSION RANGE(40),IN(84),S(84),LZ(84),DMIN(84),MINCOL(84),
1LAB(84),EDSQ(2700),ED(84)
EQUIVALENCE(ID,EDSQ),(AD,S),(IBUF,LAB)
C
CJA MOD
C
DIMENSION DUMTTL(8),IOPFIL(34)
COMMON KI,KO,IDSK1,IDSK2
C ICON IPWR IDISS IPDISS ISORT
DATA IOPFIL /4*0 , 4*0 , 4,3*0 , 4*0 , 4*0 ,
1 2,3*0, 2*0 , 1,3*0 , 4*99999 /
C IDEN ITT IPCA ICD
DATA DUMTTL /'ENT1/ATTR DATA ENT1/ENT1 ANAL. (TAXAN)'/
C
C **************************************************************
C READ AND PRINT HEADING.
C
CALL INIT
CALL DATPUT('TAXAN','DATFL',"150000127)
1100 READ(KI,1101,END=3000) TITLE
1101 FORMAT(13A5)
201 FORMAT(1H1,13A5)
C
C READ CONTROL CARDS
C
READ(KI,103) (IBUF(IND),IND=1,16)
103 FORMAT(16A5)
DECODE(80,104,IBUF(1)) NENT,NB,NN,ND,NO,J1,J2,PLOT,ORDIN
104 FORMAT(10I)
IF(PLOT .GT. 0) IOPFIL(21)=0
IF(ORDIN .GT. 0) IOPFIL(27)=0
C
CJA MODS
C
IDSK1=20
IDSK2=21
NROWSZ=NB+NN+ND+NO
NTIMZ=1
CALL OFILE(IDSK1,'AVGE')
WRITE(IDSK1)NENT,NTIMZ,NROWSZ,TITLE
WRITE(IDSK1)IOPFIL
CALL RELEAS(IDSK1)
CALL OFILE(IDSK1,'DENI1')
CALL OFILE(IDSK2,'DEN1')
WRITE(IDSK1)NENT
WRITE(IDSK2)NENT
WRITE(IDSK2)DUMTTL,TITLE
WRITE(IDSK1)DUMTTL,TITLE
C
C
C
WRITE(KO,201)TITLE
WRITE(KO,105) NENT,NB,NN,ND,NO
105 FORMAT('0CONTROL PARAMETERS ',10I5)
C
C MULTISTATE SIZES
C
NM = ND + NO
IF(NM.EQ.0)GO TO 1401
READ(KI,1102)(ISIZE(I),I=1,NM)
1102 FORMAT(40(A1,1X))
203 FORMAT('0CLASSES IN MULTISTATES ',40(A1,1X))
WRITE(KO,203) (ISIZE(I),I = 1 , NM)
DO 301 I=1,NM
IF(ISIZE(I).EQ.' ') GO TO 302
IF(ISIZE(I) - '*')303,302,303
302 ISIZE(I)=-9
GO TO 301
303 DECODE(1,107,ISIZE(I)),ISIZE(I)
301 CONTINUE
C
C THE TRUE DATA CAN NOW BE READ IN.
C THESE ARE STORED IN TWO ONE DIMENSIONAL ARRAYS.
C
1401 NM = NN + NO
NQ = NB + ND
C
C LOOP FOR NENT ENTITIES
C
DO 40 J = 1,NENT
C
C BINARIES
C
WRITE(KO,1422)J
1422 FORMAT('0ENTITY ',I2,' DATA')
IF(NB.LE.0)GO TO 11
READ(KI,1103) (IBUF(I), I=1,NB)
WRITE(KO,11030)(IBUF(I),I=1,NB)
11030 FORMAT(' A ',60A1)
1103 FORMAT(60A1)
205 FORMAT(10X,70(A1,1X))
JNN = (J-1)*NQ
DO 10 I= 1,NB
C CALCULATION OF PLACE IN ARRAY FOR EACH PIECE OF DATA.
IRUN = I + JNN
C DETERMINE IF THE NUMBER IS KNOWN (I.E. A NON NEGATIVE NUMBER.)
IF(IBUF(I).EQ.' ')GO TO 7
2 IF(IBUF(I)- '*') 8,7,8
C IF VALUE IS UNKNOWN RECORD IT AS -9
7 ID(IRUN) = -9
GO TO 10
C KNOWN BINARY DATA ARE STORED AS ONE OR TWO.
8 DECODE(1,107,IBUF(I))ID(IRUN)
ID(IRUN)=ID(IRUN) + 1
10 CONTINUE
C
C NUMERICS
C
11 IF(NN.LE.0)GO TO 21
READ(KI,1105)(BUF(I),I=1,NN)
WRITE(KO,11050) (BUF(I),I=1,NN)
11050 FORMAT(' B ',30F)
1105 FORMAT(30F)
206 FORMAT(1X,15F8.3,10X,15F8.3)
JNN = (J-1)*NM
DO 20 I = 1,NN
IRUN = I + JNN
12 IF(BUF(I)-.00001)15,15,16
C A NEGATIVE OR ZERO VALUE IS CONSIDERED AS MISSING DATA.
15 AD(IRUN) = -999.
GO TO 20
16 AD(IRUN) = BUF(I)
C MISSING DATA IS STORED AS -999.
IA = IRUN
20 CONTINUE
C
C DISORDERED MULTISTATE
C
21 IF(ND.LE.0) GO TO 31
READ(KI,106) (IBUF(I),I=1,ND)
WRITE(KO,1060) (IBUF(I), I=1,ND)
1060 FORMAT(' C ',(40(A1,1X)))
106 FORMAT(40(A1,1X))
JNN = (J-1)*NQ + NB
DO 30 I = 1 , ND
IRUN = I + JNN
IF(IBUF(I).EQ.' ')GO TO 25
22 IF(IBUF(I) - '*')25,25,27
25 ID (IRUN) = -9
GO TO 30
27 DECODE(1,107,IBUF(I))ID(IRUN)
107 FORMAT(I1)
30 CONTINUE
IS = IRUN
C
C ORDERED MULTISTATE
C
31 IF(NO.LE.0) GO TO 40
READ(KI,106) (IBUF(I), I=1,NO)
WRITE(KO,10600)(IBUF(I),I=1,NO)
10600 FORMAT(' D ',(40(A1,1X)))
JNN = NN + (J-1) * NM
DO 40 I = 1 , NO
IRUN = I + JNN
IF(IBUF(I).EQ.' ')GO TO 36
IF(IBUF(I) - '*')35,36,35
36 AD(IRUN) =-999.
GO TO 40
35 DECODE(1,107,IBUF(I))IAA
AD(IRUN) = IAA
40 CONTINUE
C
C THIS CONCLUDES INPUT
C
AREL = 0.0
C COMPUTE RANGES
IF(NM.LE.0) GO TO 52
DO 50 I = 1, NM
BIG = 0.0
SMALL = 1000.0
SUM = 0.0
SS = 0.0
C = 0.0
IRUN = I
DO 45 J = 1,NENT
C FIND VALUE
ADI = AD(IRUN)
IF (ADI + 900.) 45,45,141
141 SUM = SUM + ADI
SS = SS + ADI *ADI
C = C + 1.
IF (ADI - BIG) 42,42,41
41 BIG = ADI
42 IF(ADI - SMALL) 43,45,45
43 SMALL = ADI
45 IRUN = IRUN + NM
IF(J2.EQ.2) RANGE(I) = SQRT(2.*(SS-(SUM*SUM)/C)/(C-1.))
50 IF(J2.EQ.1)RANGE(I) = BIG - SMALL
52 IF(NO.LE.0) GO TO 2000
C
C CALCULATE ADDRESS IN EDSQ OF ZEROTH ELEMENT OF EACH ROW.
C STORE THESE IN ARRAY LZ
2000 LZ(1)=0
WRITE(KO,1450)
1450 FORMAT(1H1,'DISSIMILARITY MATRIX',/)
C CALCULATE THE SIMILARITY MATRIX (LOWER TRIANGLE)
DO 2080 I = 2, NENT
LZ(I) = LZ(I-1) + I-2
IM = I-1
DM = 1000.0
MCOL = 0
DO 2060 J=1,IM
DIF = 0.0
AREL = 0.0
IF(NQ.LE.0) GO TO 2030
INN = IM*NQ
JNN = (J-1) * NQ
DO 2020 K=1,NQ
IRUN = K + INN
JRUN = K + JNN
II = ID(IRUN)
JJ = ID(JRUN)
IF(II + JJ.LE.0) GO TO 2020
IF(II.NE.JJ) DIF = DIF + 1.0
AREL = AREL + 1.0
2020 CONTINUE
2030 IF(NM.LE.0) GO TO 2050
INN = IM*NM
JNN = (J-1)*NM
DO 2040 K = 1,NM
IF(RANGE(K).LE.0.01) GO TO 2040
IRUN = K + INN
JRUN = K + JNN
AA = AD(IRUN)
BB = AD(JRUN)
IF(AA + BB.LE.0.0) GO TO 2040
DIF = DIF + ((AA - BB)/RANGE(K))**2
AREL = AREL + 1.0
2040 CONTINUE
C IF NO COMPARISON SET DISTANCE = 100000.0
2050 AA = 100000.0
IF(AREL.GT.0.0) AA = DIF/AREL
ED(J) = AA
C FIND SMALLEST ELEMENT IN ROW I, STORE IT IN
C DMIN(I), STORE ITS COLUMN ADDRESS IN MINCOL(I)
IF(AA.GE.DM) GO TO 2060
DM = AA
MCOL = J
2060 CONTINUE
DMIN(I) = DM
MINCOL(I) = MCOL
WRITE(13) (ED(J), J = 1,IM)
C
CJA AGAIN
C
WRITE(IDSK2)(ED(J),J=1,IM)
WRITE(KO,20600) (ED(J), J=1,IM)
20600 FORMAT((15F8.4))
2080 CONTINUE
C READ MATRIX BACK INTO CORE
REWIND 13
WRITE(KO,1451)
1451 FORMAT(1H1,'FUSIONS ON DENDROGRAM',/)
DO 2090 I = 2,NENT
IRUN = LZ(I) + 1
IM = IRUN + I-2
2090 READ(13) (EDSQ(J), J = IRUN,IM)
C INITIALIZE LAB(I) = CLUSTER LABEL CORRESPONDING TO ROW I
C INTIALIZE IN = CLUSTER SIZE, S = SSWC
DO 2100 I = 1,NENT
IN(I) = 1
S(I) = 0.0
2100 LAB(I) = I
NNENT = NENT
C
C FUSION CYCLE BEGINS. FIND SMALLEST DIJ
2110 NNENT = NNENT + 1
IF(NNENT.GE.NENT+NENT) GO TO 1100
DIJ = 1000.0
DO 2120 K=2,NENT
IF(DMIN(K).GE.DIJ) GO TO 2120
DIJ = DMIN(K)
I = K
J = MINCOL(K)
2120 CONTINUE
IF (DIJ.GE.1000.0) GO TO 1100
INI = IN(I)
INJ = IN(J)
WRITE(KO,1001) LAB(J), LAB(I), NNENT, DIJ, INJ, INI
1001 FORMAT(I4,' + ',I4,' = ', I4, F8.3, 2I8,2X,F6.3,I8)
WRITE(IDSK1)LAB(J),LAB(I),NNENT,DIJ
C FUSION
A = INI
B = INJ
C = A + B
SI=S(I)
SINEW = (C-1.)*DIJ
PIJ = C*SINEW
S(I) = SINEW
SINI = A*SI
SJNJ = A*S(J)
IN(I) = INI + INJ
IN(J) = 0
LAB(I) = NNENT
C DELETE OBSOLETE ROW MINIMA
DMIN(J) = 1000.0
C PATCH BY C J MCGOVERN JAN 1977
DMIN(I) = 1000.0
C END PATCH
DO 2130 K=J,NENT
IF(IN(K).LE.0) GO TO 2130
IF(MINCOL(K).NE.I.AND.MINCOL(K).NE.J) GO TO 2130
C PATCH BY C J MCGOVERN JAN 1977
IF (I .EQ. K) GO TO 2130
C END PATCH
LMIN = 0
DMK = 1000.0
KM = K-1
IF(KM.EQ.0) GO TO 2129
DO 2128 L = 1,KM
IF(L.EQ.I) GO TO 2128
IF(IN(L).LE.0) GO TO 2128
KL = LZ(K) + L
IF(EDSQ(KL).GE.DMK) GO TO 2128
DMK = EDSQ(KL)
LMIN = L
2128 CONTINUE
2129 DMIN(K) = DMK
MINCOL(K) = LMIN
2130 CONTINUE
C RECOMPUTE DISTANCES IN ROW I AND COLUMN I
DO 2320 K = 1,NENT
C EXTRACT DIK AMD DJK UNLESS K IS DISCARDED
D = IN(K)
IF(D.LE.0.0) GO TO 2320
IF(K-J) 2140,2320,2150
2140 JK = LZ(J) + K
GO TO 2160
2150 JK = LZ(K) + J
2160 DJK = EDSQ(JK)
IF(K-I) 2170,2320,2180
2170 IK = LZ(I) + K
NROW = I
NCOL = K
GO TO 2190
2180 IK = LZ(K) + I
NROW = K
NCOL = I
2190 DIK = EDSQ(IK)
GO TO (2210,2220,2230,2240,2250,2260,2270), J1
C 1 - NEAREST NEIGHBOUR
2210 IF(DIK - DJK) 2212,2212,2214
2212 DNEW = DIK
GO TO 2300
2214 DNEW = DJK
GO TO 2300
C 2 - FURTHEST NEIGHBOUR
2220 IF(DIK - DJK) 2214,2212,2212
C 3 - WEIGHTED AVERAGE
2230 DNEW = (A*DIK + B*DJK)/C
GO TO 2300
C 4 - SIMPLE AVERAGE
2240 DNEW = (DIK + DJK)*0.5
GO TO 2300
C 5 - CENTROID
2250 DNEW = (A*DIK + B*DJK - A*B*DIJ/C)/C
GO TO 2300
C 6 - INCREMENTAL SS
2260 DNEW = ((A+D)*DIK + (B+D)*DJK -D*DIJ)/(C+D)
GO TO 2300
C 7 - VARIANCE
2270 SKNK =D*S(K)
PIJK = C+D
PIJK = PIJK*PIJK-PIJK
PIK = A+D
PIK = (PIK*PIK - PIK)*DIK
PJK = B+D
PJK = (PJK*PJK - PJK)*DJK
DNEW = (PIK+PJK+PIJ-SINI-SJNJ-SKNK)/PIJK
C FINALLY, UPDATE DMIN IF REQUIRED
2300 EDSQ(IK) = DNEW
IF(DNEW.GE.DMIN(NROW)) GO TO 2320
DMIN(NROW) = DNEW
MINCOL(NROW) = NCOL
2320 CONTINUE
GO TO 2110
3000 CALL RELEAS(IDSK1)
CALL RELEAS(IDSK2)
CALL RUN7
END