SUBROUTINE enum(r, c, n, m, ifault) ! Code converted using TO_F90 by Alan Miller ! Date: 2002-02-18 Time: 10:00:55 ! ALGORITHM AS 205 APPL. STATIST. (1984) VOL 33, NO. 3. ! ENUMERATES ALL R*C CONTINGENCY TABLES WITH GIVEN ROW TOTALS N(I) ! AND COLUMN TOTALS M(J) AND CALCULATES THE HYPERGEOMETRIC ! PROBABILITY OF EACH TABLE. ! FOR TABLES HAVING TWO OR MORE ROW SUMS REPEATED, EQUIVALENT ! TABLES DIFFERING ONLY BY A ROW PERMUTATION ARE NOT SEPARATELY ! ENUMERATED. A REPRESENTATIVE OF EACH EQUIVALENCE CLASS IS ENUMERATED ! AND THE MULTIPLICITY OF EACH CLASS CALCULATED. ! FOR EACH TABLE ENUMERATED, SUBROUTINE EVAL IS CALLED TO CARRY OUT ! CALCULATIONS ON THE TABLE. IMPLICIT NONE INTEGER, INTENT(IN OUT) :: r INTEGER, INTENT(IN OUT) :: c INTEGER, INTENT(IN OUT) :: n(10) INTEGER, INTENT(IN OUT) :: m(10) INTEGER, INTENT(OUT) :: ifault INTERFACE SUBROUTINE eval(table, r, c, n, m, prob, mult) IMPLICIT NONE INTEGER, INTENT(IN) :: table(10,10) INTEGER, INTENT(IN) :: r INTEGER, INTENT(IN) :: c INTEGER, INTENT(IN) :: n(10) INTEGER, INTENT(IN) :: m(10) REAL, INTENT(IN) :: prob INTEGER, INTENT(IN) :: mult END SUBROUTINE eval END INTERFACE INTEGER :: table(10,10), bound(10,10), z(10) INTEGER :: reps(10), mult(10), repsr, repsc, maxrc, multr, multc INTEGER :: ntotal, ntot, rm, cm, left, rowbnd, rowsum INTEGER :: i, j, ii, ij, iip, iim, ip, jj, jjm, jkeep, jm, jnext, jp, k, KEEP REAL :: prob(10,10), flogm(201), factlm(201) REAL :: prob0 LOGICAL :: rept(10), reptc(10), ieqim ! LOCAL VARIABLES - ! TABLE(I,J) - (I,J)-TH ENTRY OF CURRENT TABLE ! NTOTAL - TOTAL OF TABLE ENTRIES ! BOUND(I,J) - CURRENT UPPER BOUND ON TABLE(I,J) TO SATISFY ! ROW AND COLUMN TOTALS ! REPT(I) - LOGICAL = TRUE IF ROW TOTALS N(I), N(I-1) ARE EQUAL ! = FALSE OTHERWISE ! REPS(I) - NUMBER OF PREVIOUS ROWS EQUAL TO ROW I ! MULT(I) - MAXIMUM NUMBER OF EQUIVALENT TABLES ! GIVEN FIRST I ROWS ! Z(J) - LOWER BOUND ON SUM OF ENTRIES USED BY ALGORITHM C ! PROB(I,J) - PARTIAL SUM OF TERMS IN LOG(P) INTEGER, PARAMETER :: maxt = 10, nmax = 201, zero = 0.0, one = 1 ! MAXT - MAXIMUM DIMENSION OF TABLE ! NMAX - MAXIMUM NUMBER OF OBSERVATIONS + 1 ! CHECK INPUT VALUES ifault = 1 IF (r > MAXT .OR. c > MAXT .OR. r <= 0 .OR. c <= 0) RETURN ifault = 3 ntotal = 0 DO i = 1, r IF (n(i) <= 0) RETURN ntotal = ntotal + n(i) END DO ntot = 0 DO j = 1, c IF (m(j) <= 0) RETURN ntot = ntot + m(j) END DO ifault = 2 IF (ntot /= ntotal) RETURN ifault = 4 IF (ntotal >= nmax) RETURN ifault = 0 ! INITIALISE FLOGM(K)=LOG(K-1), FACTLM(K)=LOG(K-1 FACTORIAL) flogm(1) = zero factlm(1) = zero DO k = 1, ntotal flogm(k+1) = LOG(REAL(k)) factlm(k+1) = factlm(k) + flogm(k+1) END DO ! CONSTANTS rm = r - 1 cm = c - 1 ! SORT ROWS AND COLUMNS INTO ASCENDING ORDER DO i = 1, rm ip = i + 1 DO ii = ip, r IF (n(i) > n(ii)) THEN KEEP = n(i) n(i) = n(ii) n(ii) = KEEP END IF END DO END DO DO j = 1, cm jp = j + 1 DO jj = jp, c IF (m(j) > m(jj)) THEN KEEP = m(j) m(j) = m(jj) m(jj) = KEEP END IF END DO END DO ! CALCULATE MULTIPLICITIES OF ROWS AND COLUMNS ! REPTC(J) = .TRUE. IF COLUMNS J AND J-1 HAVE THE SAME TOTAL ! REPT(I) = .TRUE. IF ROWS I AND I-1 HAVE THE SAME TOTAL multc = one repsc = one reptc(1) = .false. DO j = 2, c reptc(j) = (m(j) == m(j-1)) IF (.NOT.reptc(j)) THEN repsc = one ELSE repsc = repsc + one multc = multc * repsc END IF END DO multr = one repsr = one rept(1) = .false. DO i = 2, r rept(i) = (n(i) == n(i-1)) IF (.NOT.rept(i)) THEN repsr = one ELSE repsr = repsr + one multr = multr * repsr END IF END DO ! IF COLUMN MULTIPLICITY EXCEEDS ROW MULTIPLICITY TRANSPOSE TABLE IF (multr < multc) THEN maxrc = MAX(r,c) DO ij = 1, maxrc KEEP = n(ij) n(ij) = m(ij) m(ij) = KEEP END DO KEEP = r r = c c = KEEP rm = r - 1 cm = c - 1 DO i = 1, r rept(i) = reptc(i) END DO multr = multc END IF ! SET UP INITIAL TABLE ! MAXIMUM MULTIPLICITY mult(1) = multr reps(1) = one ! CONSTANT TERM IN PROBABILITY prob0 = -factlm(ntotal+1) DO i = 1, r ii = n(i) prob0 = prob0 + factlm(ii+1) END DO DO j = 1, c jj = m(j) prob0 = prob0 + factlm(jj+1) END DO ! CALCULATE BOUNDS ON ROW 1 DO j = 1, c bound(1,j) = m(j) END DO ! FOR EACH I FIND GREATEST I-TH ROW SATISFYING BOUNDS DO i = 1, r IF (i /= 1) prob0 = prob(i-1,c) left = n(i) ! ELEMENTS OF ROW I ieqim = rept(i) DO j = 1, cm ij = MIN(left,bound(i,j)) table(i,j) = ij IF (j == 1) prob(i,j) = prob0 - factlm(ij+1) IF (j /= 1) prob(i,j) = prob(i,j-1) - factlm(ij+1) left = left - table(i,j) IF (i < r) bound(i+1,j) = bound(i,j) - table(i,j) IF (left == 0) GO TO 160 IF (ieqim) ieqim = table(i,j) == table(i-1,j) END DO table(i,c) = left prob(i,c) = prob(i,cm) - factlm(left+1) IF (i < r) bound(i+1,c) = bound(i,c) - left GO TO 180 160 jp = j + 1 DO jj = jp, c table(i,jj) = 0 prob(i,jj) = prob(i,jj-1) bound(i+1,jj) = bound(i,jj) END DO 180 IF (i /= 1) THEN mult(i) = mult(i-1) reps(i) = one IF (ieqim) THEN reps(i) = reps(i-1) + one mult(i) = mult(i) / reps(i) END IF END IF END DO ! CALL EVAL FOR TABLE 1 CALL eval(table, r, c, n, m, prob(r,c), mult(r)) ! COMMENCE ENUMERATION OF REMAINING TABLES ! START OF MAIN LOOP 200 i = r 210 i = i - 1 ! IF I = 0 NO MORE TABLES ARE POSSIBLE IF (i == 0) RETURN j = cm left = table(i,c) rowbnd = bound(i,c) ! TRY TO DECREASE ELEMENT (I,J) 220 IF (table(i,j) <= 0 .OR. left >= rowbnd) THEN ! ELEMENT (I,J) CANNOT BE DECREASED - TRY (I,J-1) IF (j == 1) GO TO 210 left = left + table(i,j) rowbnd = rowbnd + bound(i,j) j = j - 1 GO TO 220 END IF ! DECREASE ELEMENT (I,J) ij = table(i,j) prob(i,j) = prob(i,j) + flogm(ij+1) table(i,j) = table(i,j) - 1 bound(i+1,j) = bound(i+1,j) + 1 ! IF ROW I WAS THE SAME AS ROW I-1 IT IS NO LONGER IF (reps(i) /= one) THEN reps(i) = one mult(i) = mult(i-1) END IF ! COMPLETE ROW I WITH THE LARGEST POSSIBLE VALUES ii = i iip = ii + 1 iim = ii - 1 jnext = j + 1 left = left + 1 GO TO 320 ! FILL UP REMAINING ROWS 230 ii = ii + 1 ! THE LAST ROW IS TREATED SEPARATELY IF (ii == r) GO TO 370 iip = ii + 1 iim = ii - 1 IF (.NOT.rept(ii)) THEN ! ROW TOTAL N(II) IS NOT A REPEAT - MAKE ROW II AS LARGE AS POSSIBLE left = n(ii) jnext = 1 ELSE ! REPEATED ROW TOTALS ! (I) IF ROW II-1 SATISFIES THE BOUNDS ON ROW II REPEAT IT DO j = 1, c IF (table(iim,j) > bound(ii,j)) GO TO 250 ij = table(iim,j) table(ii,j) = ij bound(iip,j) = bound(ii,j) - table(ii,j) IF (j == 1) prob(ii,j) = prob(iim,c) - factlm(ij+1) IF (j /= 1) prob(ii,j) = prob(ii,j-1) - factlm(ij+1) END DO ! ROW II IS A REPEAT OF ROW II-1 reps(ii) = reps(iim) + one mult(ii) = mult(iim) / reps(ii) GO TO 230 ! ELEMENT J OF ROW II-1 WAS TOO BIG ! CONSTRUCT THE SEQUENCE Z(J) OF LOWER BOUNDS 250 IF (j <= 1) THEN ! IF J=1 THE BOUNDS ARE SATISFIED AUTOMATICALLY ij = bound(ii,1) table(ii,1) = ij prob(ii,1) = prob(iim,c) - factlm(ij+1) jnext = 2 left = n(ii) - table(ii,1) bound(iip,1) = 0 ELSE z(j) = n(ii) jm = j - 1 IF (j /= c) THEN jp = j + 1 DO jj = jp, c z(j) = z(j) - bound(ii,jj) END DO END IF DO jjm = 1, jm jj = j - jjm z(jj) = z(jj+1) - bound(ii,jj+1) END DO ! (II) IF THE CUMULATIVE TOTALS OF ROW II-1 ALL EXCEED THE BOUNDS Z(J) ! MAKE ELEMENT (II,J) EQUAL TO ITS BOUND rowsum = 0 jkeep = 0 DO jj = 1, jm rowsum = rowsum + table(iim,jj) IF (rowsum < z(jj)) GO TO 300 IF (rowsum > z(jj) .AND. table(iim,jj) > 0) jkeep = jj END DO table(ii,j) = bound(ii,j) bound(iip,j) = 0 ij = table(ii,j) prob(ii,j) = prob(ii,jm) - factlm(ij+1) reps(ii) = one mult(ii) = mult(iim) ! COMPLETE ROW II WITH THE LARGEST POSSIBLE ELEMENTS jnext = jp left = n(ii) DO jj = 1, j left = left - table(ii,jj) END DO GO TO 320 ! (III) THE CUMULATIVE SUMS VIOLATE THE BOUNDS ! IF NO ELEMENT OF ROW II-1 CAN BE CHANGED TO SATISFY THE BOUNDS ! NO SUITABLE ROW II IS POSSIBLE ! IN THAT CASE GO BACK AND TRY DECREASING ROW II-1 300 IF (jkeep == 0) THEN i = ii GO TO 210 END IF ! ELEMENT (II,JKEEP) CAN BE DECREASED bound(iip,jkeep) = bound(iip,jkeep) + 1 ij = table(ii,jkeep) prob(ii,jkeep) = prob(ii,jkeep) + flogm(ij+1) table(ii,jkeep) = table(ii,jkeep) - 1 ! COMPLETE THE ROW jnext = jkeep + 1 left = n(ii) DO jj = 1, jkeep left = left - table(ii,jj) END DO END IF END IF ! ROW II IS COMPLETE UP TO ELEMENT JNEXT-1 ! MAKE THE REMAINING ELEMENTS AS LARGE AS POSSIBLE ! (THIS SECTION OF CODE IS USED FOR EVERY ROW, REPEATED OR NOT) 320 IF (jnext /= c) THEN DO j = jnext, cm table(ii,j) = MIN(left,bound(ii,j)) left = left - table(ii,j) bound(iip,j) = bound(ii,j) - table(ii,j) ij = table(ii,j) IF (j == 1) prob(ii,j) = prob(iim,c) - factlm(ij+1) IF (j /= 1) prob(ii,j) = prob(ii,j-1) - factlm(ij+1) IF (left == 0) GO TO 340 END DO END IF table(ii,c) = left prob(ii,c) = prob(ii,cm) - factlm(left+1) bound(iip,c) = bound(ii,c) - left GO TO 360 340 jp = j + 1 DO jj = jp, c table(ii,jj) = 0 prob(ii,jj) = prob(ii,jj-1) bound(iip,jj) = bound(ii,jj) END DO 360 reps(ii) = one IF (ii > 1) mult(ii) = mult(iim) GO TO 230 ! THE FINAL ROW 370 IF (.NOT.rept(r)) THEN ! NOT A REPEAT - SET ROW R EQUAL TO ITS BOUNDS ij = bound(r,1) table(r,1) = ij prob(r,1) = prob(rm,c) - factlm(ij+1) DO j = 2, c ij = bound(r,j) table(r,j) = ij prob(r,j) = prob(r,j-1) - factlm(ij+1) END DO mult(r) = mult(rm) ELSE ! ROW TOTAL R IS A REPEAT - ENSURE THAT IT IS LESS THAN ROW R-1 DO j = 1, c IF (bound(r,j) > table(rm,j)) GO TO 400 ij = bound(r,j) table(r,j) = ij IF (j == 1) prob(r,j) = prob(rm,c) - factlm(ij+1) IF (j /= 1) prob(r,j) = prob(r,j-1) - factlm(ij+1) IF (table(r,j) /= table(rm,j)) GO TO 410 END DO ! ROW R IS A REPEAT OF ROW R-1 reps(r) = reps(rm) + one mult(r) = mult(rm) / reps(r) GO TO 430 ! IF ROW R WOULD BE BIGGER THAN ROW R-1 GO BACK AND TRY ! DECREASING ROW R-2 400 i = rm GO TO 210 ! ROW R IS ALREADY LESS THEN ROW R-1 SO NO MORE CHECKS ARE NEEDED 410 jp = j + 1 DO jj = jp, c ij = bound(r,jj) table(r,jj) = ij prob(r,jj) = prob(r,jj-1) - factlm(ij+1) END DO mult(r) = mult(rm) END IF ! THE TABLE IS COMPLETE - CALL SUBROUTINE EVAL 430 CALL eval(table, r, c, n, m, prob(r,c), mult(r)) ! END OF MAIN LOOP GO TO 200 END SUBROUTINE enum SUBROUTINE eval(table, r, c, n, m, prob, mult) IMPLICIT NONE INTEGER, INTENT(IN) :: table(10,10) INTEGER, INTENT(IN) :: r INTEGER, INTENT(IN) :: c INTEGER, INTENT(IN) :: n(10) INTEGER, INTENT(IN) :: m(10) REAL, INTENT(IN) :: prob INTEGER, INTENT(IN) :: mult RETURN END SUBROUTINE eval