4 static void f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
5 double *expect, double *percnt, double *emin, double
6 *prt, double *pre, double *fact, int *ico, int
7 *iro, int *kyy, int *idif, int *irn, int *key,
8 int *ldkey, int *ipoin, double *stp, int *ldstp,
9 int *ifrq, double *dlp, double *dsp, double *tm,
10 int *key2, int *iwk, double *rwk);
11 static void f3xact(int *nrow, int *irow, int *ncol, int *icol,
12 double *dlp, int *mm, double *fact, int *ico, int
13 *iro, int *it, int *lb, int *nr, int *nt, int
14 *nu, int *itc, int *ist, double *stv, double *alen,
16 static void f4xact(int *nrow, int *irow, int *ncol, int *icol,
17 double *dsp, double *fact, int *icstk, int *ncstk,
18 int *lstk, int *mstk, int *nstk, int *nrstk, int
19 *irstk, double *ystk, const double *tol);
20 static void f5xact(double *pastp, const double *tol, int *kval, int *key,
21 int *ldkey, int *ipoin, double *stp, int *ldstp,
22 int *ifrq, int *npoin, int *nr, int *nl, int
23 *ifreq, int *itop, int *ipsh);
24 static void f6xact(int *nrow, int *irow, int *iflag, int *kyy,
25 int *key, int *ldkey, int *last, int *ipn);
26 static void f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
28 static void f8xact(int *irow, int *is, int *i1, int *izero, int *myNew);
29 static double f9xact(int *n, int *mm, int *ir, double *fact);
30 static void f10act(int *nrow, int *irow, int *ncol, int *icol,
31 double *val, int *xmin, double *fact, int *nd,
33 static void f11act(int *irow, int *i1, int *i2, int *myNew);
34 static void prterr(int icode, char *mes);
35 static int iwork(int iwkmax, int *iwkpt, int number, int itype);
36 // void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
37 // double *expect, double *percnt, double *emin, double *prt,
38 // double *pre, /* myNew in C : */ int *workspace);
39 static void isort(int *n, int *ix);
40 static double gammds(double *y, double *p, int *ifault);
41 static double alogam(double *x, int *ifault);
44 /* The only public function : */
45 void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
46 double *expect, double *percnt, double *emin, double *prt,
47 double *pre, /* myNew in C : */ int *workspace) {
50 ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
51 THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
52 VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
53 -----------------------------------------------------------------------
55 Purpose: Computes Fisher's exact test probabilities and a hybrid
56 approximation to Fisher exact test probabilities for a
57 contingency table using the network algorithm.
58 Usage: CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
61 NROW - The number of rows in the table. (Input)
62 NCOL - The number of columns in the table. (Input)
63 TABLE - NROW by NCOL matrix containing the contingency
65 LDTABL - Leading dimension of TABLE exactly as specified
66 in the dimension statement in the calling
68 EXPECT - Expected value used in the hybrid algorithm for
69 deciding when to use asymptotic theory
70 probabilities. (Input)
71 If EXPECT <= 0.0 then asymptotic theory probabilities
72 are not used and Fisher exact test probabilities are
73 computed. Otherwise, if PERCNT or more of the cells in
74 the remaining table have estimated expected values of
75 EXPECT or more, with no remaining cell having expected
76 value less than EMIN, then asymptotic chi-squared
77 probabilities are used. See the algorithm section of the
78 manual document for details.
79 Use EXPECT = 5.0 to obtain the 'Cochran' condition.
80 PERCNT - Percentage of remaining cells that must have
81 estimated expected values greater than EXPECT
82 before asymptotic probabilities can be used. (Input)
83 See argument EXPECT for details.
84 Use PERCNT = 80.0 to obtain the 'Cochran' condition.
85 EMIN - Minimum cell estimated expected value allowed for
86 asymptotic chi-squared probabilities to be used. (Input)
87 See argument EXPECT for details.
88 Use EMIN = 1.0 to obtain the 'Cochran' condition.
89 PRT - Probability of the observed table for fixed
90 marginal totals. (Output)
91 PRE - Table p-value. (Output)
92 PRE is the probability of a more extreme table,
93 where `extreme' is in a probabilistic sense.
94 If EXPECT < 0 then the Fisher exact probability
95 is returned. Otherwise, an approximation to the
96 Fisher exact probability is computed based upon
97 asymptotic chi-squared probabilities for ``large''
98 table expected values. The user defines ``large''
99 through the arguments EXPECT, PERCNT, and EMIN.
102 1. For many problems one megabyte or more of workspace can be
103 required. If the environment supports it, the user should begin
104 by increasing the workspace used to 200,000 units.
105 2. In FEXACT, LDSTP = 30*LDKEY. The proportion of table space used
106 by STP may be changed by changing the line MULT = 30 below to
108 3. FEXACT may be converted to single precision by setting IREAL = 3,
109 and converting all DOUBLE PRECISION specifications (except the
110 specifications for RWRK, IWRK, and DWRK) to REAL. This will
111 require changing the names and specifications of the intrinsic
112 functions ALOG, AMAX1, AMIN1, EXP, and REAL. In addition, the
113 machine specific constants will need to be changed, and the name
114 DWRK will need to be changed to RWRK in the call to F2XACT.
115 4. Machine specific constants are specified and documented in F2XACT.
116 A missing value code is specified in both FEXACT and F2XACT.
117 5. Although not a restriction, is is not generally practical to call
118 this routine with large tables which are not sparse and in
119 which the 'hybrid' algorithm has little effect. For example,
120 although it is feasible to compute exact probabilities for the
125 computing exact probabilities for a similar table which has been
126 enlarged by the addition of an extra row (or column) may not be
128 -----------------------------------------------------------------------
131 /* CONSTANT Parameters : */
133 /* To increase the length of the table of paste path lengths relative
134 to the length of the hash table, increase MULT.
137 /* AMISS is a missing value indicator which is returned when the
138 probability is not defined.
140 const double amiss = -12345.;
142 Set IREAL = 4 for DOUBLE PRECISION
143 Set IREAL = 3 for SINGLE PRECISION
148 /* System generated locals */
150 /* Local variables */
151 int nco, nro, ntot, numb, iiwk, irwk;
152 int i, j, k, kk, ldkey, ldstp, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
153 int i3a, i3b, i3c, i9a, iwkmax, iwkpt;
155 /* Workspace Allocation (freed at end) */
157 iwkmax = 2 * (int) (*workspace / 2);
158 // equiv = (double *) R_alloc(iwkmax / 2, sizeof(double));
159 equiv = (double *) calloc(iwkmax / 2, sizeof(double));
161 /* The check could never happen with Calloc!
162 equiv = Calloc(iwkmax / 2, double);
164 prterr(0, "Can not allocate specified workspace");
168 #define iwrk ((int *)equiv)
169 #define rwrk ((float *)equiv)
171 /* Parameter adjustments */
172 table -= *ldtabl + 1;
178 prterr(1, "NROW must be less than or equal to LDTABL.");
181 for (i = 1; i <= *nrow; ++i) {
182 for (j = 1; j <= *ncol; ++j) {
183 if (table[i + j * *ldtabl] < 0.)
184 prterr(2, "All elements of TABLE must be positive.");
185 ntot = (int) (ntot + table[i + j * *ldtabl]);
189 prterr(3, "All elements of TABLE are zero.\n"
190 "PRT and PRE are set to missing values.");
196 nco = max(*nrow, *ncol);
197 nro = *nrow + *ncol - nco;/* = min(*nrow, *ncol) */
198 k = *nrow + *ncol + 1;
202 i1 = iwork(iwkmax, &iwkpt, ikh, i_real);
203 i2 = iwork(iwkmax, &iwkpt, nco, i_int);
204 i3 = iwork(iwkmax, &iwkpt, nco, i_int);
205 i3a = iwork(iwkmax, &iwkpt, nco, i_int);
206 i3b = iwork(iwkmax, &iwkpt, nro, i_int);
207 i3c = iwork(iwkmax, &iwkpt, nro, i_int);
208 ikh = max(k * 5 + (kk << 1), nco * 7 + 800);
209 iiwk= iwork(iwkmax, &iwkpt, ikh, i_int);
210 ikh = max(nco + 401, k);
211 irwk= iwork(iwkmax, &iwkpt, ikh, i_real);
214 What follows below splits the remaining amount iwkmax - iwkpt of
215 (int) workspace into hash tables as follows.
217 INT 2 * ldkey i4 i5 i11
218 REAL 2 * ldkey i8 i9 i10
221 Hence, we need ldkey times
222 3 * 2 + 3 * 2 * s + 2 * mult * s + 6 * mult
223 chunks of integer memory, where s = sizeof(REAL) / sizeof(INT).
224 If doubles are used and are twice as long as ints, this gives
226 so that the value of ldkey can be obtained by dividing available
227 (int) workspace by this number.
229 In fact, because iwork() can actually s * n + s - 1 int chunks
230 when allocating a REAL, we use ldkey = available / numb - 1.
233 Can we always assume that sizeof(double) / sizeof(int) is 2?
236 if (i_real == 4) { /* Double precision reals */
237 numb = 18 + 10 * mult;
238 } else { /* Single precision reals */
239 numb = (mult << 3) + 12;
241 ldkey = (iwkmax - iwkpt) / numb - 1;
242 ldstp = mult * ldkey;
243 ikh = ldkey << 1; i4 = iwork(iwkmax, &iwkpt, ikh, i_int);
244 ikh = ldkey << 1; i5 = iwork(iwkmax, &iwkpt, ikh, i_int);
245 ikh = ldstp << 1; i6 = iwork(iwkmax, &iwkpt, ikh, i_real);
246 ikh = ldstp * 6; i7 = iwork(iwkmax, &iwkpt, ikh, i_int);
247 ikh = ldkey << 1; i8 = iwork(iwkmax, &iwkpt, ikh, i_real);
248 ikh = ldkey << 1; i9 = iwork(iwkmax, &iwkpt, ikh, i_real);
249 ikh = ldkey << 1; i9a = iwork(iwkmax, &iwkpt, ikh, i_real);
250 ikh = ldkey << 1; i10 = iwork(iwkmax, &iwkpt, ikh, i_int);
252 /* To convert to double precision, change RWRK to DWRK in the next CALL.
294 -----------------------------------------------------------------------
296 Purpose: Computes Fisher's exact test for a contingency table,
297 routine with workspace variables specified.
298 Usage: F2XACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
299 EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF,
300 IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ,
301 DLP, DSP, TM, KEY2, IWK, RWK)
302 -----------------------------------------------------------------------
305 f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
306 double *expect, double *percnt, double *emin, double *prt,
307 double *pre, double *fact, int *ico, int *iro, int *kyy,
308 int *idif, int *irn, int *key, int *ldkey, int *ipoin,
309 double *stp, int *ldstp, int *ifrq, double *dlp, double *dsp,
310 double *tm, int *key2, int *iwk, double *rwk)
312 /* IMAX is the largest representable int on the machine. */
313 const int imax = SINT_MAX;
314 // const int imax = 2147483647; //xx: I DONÂ¥T like this, and
315 // thanks to the hint from Jason Turner I don't do it anymore. (R.D-U).
317 /* AMISS is a missing value indicator which is returned when the
318 probability is not defined. */
319 const double amiss = -12345.;
321 /* TOL is chosen as the square root of the smallest relative spacing. */
323 const static double tol = 3.45254e-7;
325 static double tol = 3.45254e-7;
327 /* EMX is a large positive value used in comparing expected values. */
328 const static double emx = 1e30;
330 /* Local variables {{any really need to be static ???}} */
331 static int kval, kmax, jkey, last, ipsh, itmp, itop, jstp, ntot,
332 jstp2, jstp3, jstp4, i, ii, j, k, n, iflag, ncell, ifreq, chisq,
333 ikkey, ikstp, ikstp2, k1, kb, kd, ks,
334 i31, i32, i33, i34, i35, i36, i37, i38, i39,
335 i41, i42, i43, i44, i45, i46, i47, i48, i310, i311,
336 nco, nrb, ipn, ipo, itp, nro, nro2;
337 static double dspt, dd, df,ddf, drn,dro, emn, obs, obs2, obs3,
346 /* Parameter adjustments */
347 table -= *ldtabl + 1;
365 /* Check table dimensions */
367 prterr(1, "NROW must be less than or equal to LDTABL.");
369 prterr(4, "NCOL must be at least 2");
371 /* Initialize KEY array */
372 for (i = 1; i <= *ldkey << 1; ++i) {
376 /* Initialize parameters */
389 /* nco := max(nrow, ncol) : */
394 /* Initialize pointers for workspace */
408 k = *nrow + *ncol + 1;
418 /* Compute row marginals and total */
420 for (i = 1; i <= *nrow; ++i) {
422 for (j = 1; j <= *ncol; ++j) {
423 if (table[i + j * *ldtabl] < -1e-4)
424 prterr(2, "All elements of TABLE must be positive.");
425 iro[i] += (int) table[i + j * *ldtabl];
431 prterr(3, "All elements of TABLE are zero.\n"
432 "PRT and PRE are set to missing values.");
437 /* Column marginals */
438 for (i = 1; i <= *ncol; ++i) {
440 for (j = 1; j <= *nrow; ++j)
441 ico[i] += (int) table[j + i * *ldtabl];
445 isort(nrow, &iro[1]);
446 isort(ncol, &ico[1]);
448 /* Determine row and column marginals.
449 Define max(nrow,ncol) =: nco >= nro := min(nrow,ncol)
452 Swap marginals if necessary to ico[1:nco] & iro[1:nro]
457 for (i = 1; i <= nco; ++i) {
467 /* Get multiplers for stack */
469 for (i = 2; i <= nro; ++i) {
470 /* Hash table multipliers */
471 if (iro[i - 1] + 1 <= imax / kyy[i - 1]) {
472 kyy[i] = kyy[i - 1] * (iro[i - 1] + 1);
477 /* Maximum product */
478 if (iro[nro - 1] + 1 <= imax / kyy[nro - 1]) {
479 kmax = (iro[nro] + 1) * kyy[nro - 1];
482 prterr(5, "The hash table key cannot be computed because "
484 "is larger than the largest representable int.\n"
485 "The algorithm cannot proceed.\n"
486 "Reduce the workspace size, or use `exact = FALSE'.");
490 /* Compute log factorials */
493 if(ntot >= 2) fact[2] = log(2.);
494 /* MM: old code assuming log() to be SLOW */
495 for (i = 3; i <= ntot; i += 2) {
496 fact[i] = fact[i - 1] + log((double) i);
499 fact[j] = fact[i] + fact[2] + fact[j / 2] - fact[j / 2 - 1];
501 /* Compute obs := observed path length */
504 for (j = 1; j <= nco; ++j) {
506 for (i = 1; i <= nro; ++i) {
508 dd += fact[(int) table[j + i * *ldtabl]];
509 ntot += (int) table[j + i * *ldtabl];
511 dd += fact[(int) table[i + j * *ldtabl]];
512 ntot += (int) table[i + j * *ldtabl];
515 obs += fact[ico[j]] - dd;
517 /* Denominator of observed table: DRO */
518 dro = f9xact(&nro, &ntot, &iro[1], fact);
519 *prt = exp(obs - dro);
520 /* Initialize pointers */
525 jstp2 = *ldstp * 3 + 1;
526 jstp3 = (*ldstp << 2) + 1;
527 jstp4 = *ldstp * 5 + 1;
530 ikstp2 = *ldstp << 1;
535 ifrq[ikstp2 + 1] = -1;
543 /* IDIF is the difference in going to the daughter */
544 for (i = 1; i <= nro; ++i)
547 /* Generate the first daughter */
550 ntot = min(n, iro[kd]);
556 while (n > 0 && kd != 1);
565 for (i = kb + 1; i <= nco; ++i)
570 /* Arc to daughter length=ICO(KB) */
571 for (i = 1; i <= nro; ++i)
572 irn[i] = iro[i] - idif[i];
577 if (irn[1] > irn[2]) {
582 } else if (nro == 3) {
586 if (irn[2] > irn[3]) {
599 } else if (ii > irn[2]) {
602 } else if (irn[2] > irn[3]) {
608 for (j = 2; j <= nro; ++j) {
612 while (ii < irn[i]) {
621 /* Adjust start for zero */
622 for (i = 1; i <= nro; ++i) {
633 /* Some table values */
634 ddf = f9xact(&nro, &n, &idif[1], fact);
635 drn = f9xact(&nro2, &ntot, &irn[nrb], fact) - dro + ddf;
638 kval = irn[1] + irn[2] * kyy[2];
639 for (i = 3; i <= nro; ++i) {
640 kval += irn[i] * kyy[i];
642 /* Get hash table entry */
643 i = kval % (*ldkey << 1) + 1;
644 /* Search for unused location */
645 for (itp = i; itp <= *ldkey << 1; ++itp) {
657 for (itp = 1; itp <= i - 1; ++itp) {
669 prterr(6, "LDKEY is too small.\n"
670 "It is not possible to give the value of LDKEY required,\n"
671 "but you could try doubling LDKEY (and possibly LDSTP).");
673 prterr(6, "LDKEY is too small for this problem.\n"
674 "Try increasing the size of the workspace.");
680 ipn = ipoin[ipo + ikkey];
681 pastp = stp[ipn + ikstp];
682 ifreq = ifrq[ipn + ikstp];
683 /* Compute shortest and longest path */
685 obs2 = obs - fact[ico[kb + 1]] - fact[ico[kb + 2]] - ddf;
686 for (i = 3; i <= k1; ++i) {
687 obs2 -= fact[ico[kb + i]];
690 dspt = obs - obs2 - ddf;
691 /* Compute longest path */
693 f3xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dlp[itp],
694 &ntot, fact, &iwk[i31], &iwk[i32], &iwk[i33],
695 &iwk[i34], &iwk[i35], &iwk[i36], &iwk[i37],
696 &iwk[i38], &iwk[i39], &rwk[i310], &rwk[i311], &tol);
697 dlp[itp] = min(0., dlp[itp]);
698 /* Compute shortest path */
700 f4xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dsp[itp], fact,
701 &iwk[i47], &iwk[i41], &iwk[i42], &iwk[i43],
702 &iwk[i44], &iwk[i45], &iwk[i46], &rwk[i48], &tol);
703 dsp[itp] = min(0., dsp[itp] - dspt);
704 /* Use chi-squared approximation? */
705 if ((irn[nrb] * ico[kb + 1]) > ntot * emn) {
707 for (i = 0; i < nro2; ++i)
708 for (j = 1; j <= k1; ++j)
709 if (irn[nrb + i] * ico[kb + j] >= ntot * *expect)
712 if (ncell * 100 >= k1 * nro2 * *percnt) {
714 for (i = 0; i < nro2; ++i)
715 tmp += (fact[irn[nrb + i]] -
716 fact[irn[nrb + i] - 1]);
718 for (j = 1; j <= k1; ++j)
719 tmp += (nro2 - 1) * (fact[ico[kb + j]] -
720 fact[ico[kb + j] - 1]);
721 df = (double) ((nro2 - 1) * (k1 - 1));
722 tmp += df * 1.83787706640934548356065947281;
723 tmp -= (nro2 * k1 - 1) * (fact[ntot] - fact[ntot - 1]);
724 tm[itp] = (obs - dro) * -2. - tmp;
726 /* tm(itp) set to a flag value */
733 obs3 = obs2 - dlp[itp];
735 if (tm[itp] == -9876.) {
742 obs2 = obs - drn - dro;
747 /* Process node with new PASTP */
750 *pre += (double) ifreq * exp(pastp + drn);
751 } else if (pastp < obs2) {
753 df = (double) ((nro2 - 1) * (k1 - 1));
755 pv = pgamma(fmax2(0., tmp + (pastp + drn) * 2.) / 2.,
756 df / 2., /*scale = */ 1.,
757 /*lower_tail = */FALSE, /*log_p = */ FALSE);
759 d1 = max(0., tmp + (pastp + drn) * 2.) / 2.;
761 pv = 1. - gammds(&d1, &d2, &ifault);
763 *pre += (double) ifreq * exp(pastp + drn) * pv;
765 /* Put daughter on queue */
767 f5xact(&d1, &tol, &kval, &key[jkey], ldkey, &ipoin[jkey],
768 &stp[jstp], ldstp, &ifrq[jstp], &ifrq[jstp2],
769 &ifrq[jstp3], &ifrq[jstp4], &ifreq, &itop, &ipsh);
773 /* Get next PASTP on chain */
774 ipn = ifrq[ipn + ikstp2];
776 pastp = stp[ipn + ikstp];
777 ifreq = ifrq[ipn + ikstp];
780 /* Generate a new daughter node */
781 f7xact(&kmax, &iro[1], &idif[1], &kd, &ks, &iflag);
787 /* Go get a new mother from stage K */
790 f6xact(&nro, &iro[1], &iflag, &kyy[1], &key[ikkey + 1], ldkey,
792 /* Update pointers */
795 /* else iflag == 3 : no additional nodes to process */
801 jkey = *ldkey - jkey + 2;
802 jstp = *ldstp - jstp + 2;
803 jstp2 = (*ldstp << 1) + jstp;
804 for (i = 1; i <= *ldkey << 1; ++i)
811 -----------------------------------------------------------------------
813 Purpose: Computes the shortest path length for a given table.
814 Usage: F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, IRO,
815 IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, TOL)
817 NROW - The number of rows in the table. (Input)
818 IROW - Vector of length NROW containing the row sums
819 for the table. (Input)
820 NCOL - The number of columns in the table. (Input)
821 ICOL - Vector of length K containing the column sums
822 for the table. (Input)
823 DLP - The longest path for the table. (Output)
824 MM - The total count in the table. (Output)
825 FACT - Vector containing the logarithms of factorials. (Input)
826 ICO - Work vector of length MAX(NROW,NCOL).
827 IRO - Work vector of length MAX(NROW,NCOL).
828 IT - Work vector of length MAX(NROW,NCOL).
829 LB - Work vector of length MAX(NROW,NCOL).
830 NR - Work vector of length MAX(NROW,NCOL).
831 NT - Work vector of length MAX(NROW,NCOL).
832 NU - Work vector of length MAX(NROW,NCOL).
833 ITC - Work vector of length 400.
834 IST - Work vector of length 400.
835 STV - Work vector of length 400.
836 ALEN - Work vector of length MAX(NROW,NCOL).
837 TOL - Tolerance. (Input)
838 -----------------------------------------------------------------------
842 f3xact(int *nrow, int *irow, int *ncol, int *icol, double *dlp,
843 int *mm, double *fact, int *ico, int *iro, int *it,
844 int *lb, int *nr, int *nt, int *nu, int *itc, int *ist,
845 double *stv, double *alen, const double *tol)
847 /* Initialized data */
848 static int ldst = 200;
852 /* Local variables */
856 static int n11, n12, ii, nn, ks, ic1, ic2, nc1, nn1;
859 static int nct, ipn, irl, key, lev, itp, nro;
861 static int nrt, kyy, nc1s;
863 /* Parameter adjustments */
878 for (i = 0; i <= *ncol; ++i) {
881 for (i = 1; i <= 400; ++i) {
887 *dlp -= fact[icol[1]];
888 for (i = 2; i <= *ncol; ++i) {
889 *dlp -= fact[icol[i]];
897 *dlp = *dlp - fact[irow[1]] - fact[irow[2]];
898 for (i = 3; i <= *nrow; ++i) {
899 *dlp -= fact[irow[i]];
905 if (*nrow * *ncol == 4) {
906 n11 = (irow[1] + 1) * (icol[1] + 1) / (*mm + 2);
908 *dlp = *dlp - fact[n11] - fact[n12] - fact[icol[1] - n11]
909 - fact[icol[2] - n12];
912 /* Test for optimal table */
915 if (irow[*nrow] <= irow[1] + *ncol) {
916 f10act(nrow, &irow[1], ncol, &icol[1], &val, &xmin, fact,
917 &lb[1], &nu[1], &nr[1]);
920 if (icol[*ncol] <= icol[1] + *nrow) {
921 f10act(ncol, &icol[1], nrow, &irow[1], &val, &xmin, fact,
922 &lb[1], &nu[1], &nr[1]);
930 /* Setup for dynamic programming */
933 if (*nrow >= *ncol) {
936 for (i = 1; i <= *nrow; ++i) {
941 for (i = 2; i <= *ncol; ++i) {
943 nt[i] = nt[i - 1] - ico[i];
950 for (i = 2; i <= *nrow; ++i) {
952 nt[i] = nt[i - 1] - ico[i];
954 for (i = 1; i <= *ncol; ++i)
957 /* Initialize pointers */
965 LnewNode: /* Setup to generate new node */
971 lb[1] = (int) ((double) ((nrt + 1) * (nct + 1)) /
972 (double) (nn + nr1 * nc1s + 1) - *tol) - 1;
973 nu[1] = (int) ((double) ((nrt + nc1s) * (nct + nr1)) /
974 (double) (nn + nr1 + nc1s)) - lb[1] + 1;
977 LoopNode: /* Generate a node */
989 alen[lev] = alen[lev - 1] + fact[lb[lev]];
996 lb[lev] = (int) ((double) ((nrt + 1) * (nct + 1)) /
997 (double) (nn1 + nr1 * nc1 + 1) - *tol);
998 nu[lev] = (int) ((double) ((nrt + nc1) * (nct + nr1)) /
999 (double) (nn1 + nr1 + nc1) - lb[lev] + 1);
1000 nr[lev] = nrt - lb[lev];
1003 alen[nco] = alen[lev] + fact[nr[lev]];
1006 v = val + alen[nco];
1008 /* Only 1 row left */
1009 v = v + fact[ico[1] - lb[1]] + fact[ico[2] - lb[2]];
1010 for (i = 3; i <= nco; ++i) {
1011 v += fact[ico[i] - lb[i]];
1016 } else if (nro == 3 && nco == 2) {
1017 /* 3 rows and 2 columns */
1018 nn1 = nn - iro[irl] + 2;
1019 ic1 = ico[1] - lb[1];
1020 ic2 = ico[2] - lb[2];
1021 n11 = (iro[irl + 1] + 1) * (ic1 + 1) / nn1;
1022 n12 = iro[irl + 1] - n11;
1023 v = v + fact[n11] + fact[n12] + fact[ic1 - n11]
1029 /* Column marginals are new node */
1030 for (i = 1; i <= nco; ++i) {
1031 it[i] = ico[i] - lb[i];
1033 /* Sort column marginals */
1035 if (it[1] > it[2]) {
1040 } else if (nco == 3) {
1044 if (it[2] > it[3]) {
1057 } else if (ii > it[2]) {
1060 } else if (it[2] > it[3]) {
1066 isort(&nco, &it[1]);
1068 /* Compute hash value */
1069 key = it[1] * kyy + it[2];
1070 for (i = 3; i <= nco; ++i) {
1071 key = it[i] + key * kyy;
1074 //PROBLEM "Bug in FEXACT: gave negative key" RECOVER(NULL_ENTRY);
1075 printf("Bug in FEXACT: gave negative key \n"); //xx:another one of my ugly kluges (R.D-U)
1078 ipn = key % ldst + 1;
1080 /* Find empty position */
1081 for (itp = ipn, ii = ks + ipn; itp <= ldst; ++itp, ++ii) {
1084 } else if (ist[ii] == key) {
1089 for (itp = 1, ii = ks + 1; itp <= ipn - 1; ++itp, ++ii) {
1092 } else if (ist[ii] == key) {
1097 prterr(30, "Stack length exceeded in f3xact.\n"
1098 "This problem should not occur.");
1100 L180: /* Push onto stack */
1108 L190: /* Marginals already on stack */
1109 stv[ii] = min(v, stv[ii]);
1114 L200: /* Pop item from stack */
1117 itp = itc[nitc + k] + k;
1122 /* Compute marginals */
1123 for (i = nco; i >= 2; --i) {
1128 /* Set up nt array */
1129 nt[1] = nn - ico[1];
1130 for (i = 2; i <= nco; ++i)
1131 nt[i] = nt[i - 1] - ico[i];
1133 /* Test for optimality (L90) */
1135 if (iro[nro] <= iro[irl] + nco) {
1136 f10act(&nro, &iro[irl], &nco, &ico[1], &val, &xmin, fact,
1137 &lb[1], &nu[1], &nr[1]);
1139 if (!xmin && ico[nco] <= ico[1] + nro)
1140 f10act(&nco, &ico[1], &nro, &iro[irl], &val, &xmin, fact,
1141 &lb[1], &nu[1], &nr[1]);
1149 } else if (nro > 2 && nst > 0) {
1150 /* Go to next level */
1165 -----------------------------------------------------------------------
1167 Purpose: Computes the longest path length for a given table.
1168 Usage: CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK,
1169 NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK,
1172 NROW - The number of rows in the table. (Input)
1173 IROW - Vector of length NROW containing the row sums for the
1175 NCOL - The number of columns in the table. (Input)
1176 ICOL - Vector of length K containing the column sums for the
1178 DSP - The shortest path for the table. (Output)
1179 FACT - Vector containing the logarithms of factorials. (Input)
1180 ICSTK - NCOL by NROW+NCOL+1 work array.
1181 NCSTK - Work vector of length NROW+NCOL+1.
1182 LSTK - Work vector of length NROW+NCOL+1.
1183 MSTK - Work vector of length NROW+NCOL+1.
1184 NSTK - Work vector of length NROW+NCOL+1.
1185 NRSTK - Work vector of length NROW+NCOL+1.
1186 IRSTK - NROW by MAX(NROW,NCOL) work array.
1187 YSTK - Work vector of length NROW+NCOL+1.
1188 TOL - Tolerance. (Input)
1189 -----------------------------------------------------------------------
1193 f4xact(int *nrow, int *irow, int *ncol, int *icol, double *dsp,
1194 double *fact, int *icstk, int *ncstk, int *lstk, int *mstk,
1195 int *nstk, int *nrstk, int *irstk, double *ystk, const double *tol)
1197 /* System generated locals */
1200 /* Local variables */
1201 int i, j, k, l, m, n, mn, ic1, ir1, ict, irt, istk, nco, nro;
1204 /* Parameter adjustments */
1217 /* Take care of the easy cases first */
1219 for (i = 1; i <= *ncol; ++i) {
1220 *dsp -= fact[icol[i]];
1225 for (i = 1; i <= *nrow; ++i) {
1226 *dsp -= fact[irow[i]];
1230 if (*nrow * *ncol == 4) {
1231 if (irow[2] <= icol[2]) {
1232 *dsp = *dsp - fact[irow[2]] - fact[icol[1]]
1233 - fact[icol[2] - irow[2]];
1235 *dsp = *dsp - fact[icol[2]] - fact[irow[1]]
1236 - fact[irow[2] - icol[2]];
1240 /* initialization before loop */
1241 for (i = 1; i <= *nrow; ++i) {
1242 irstk[i + *nrow] = irow[*nrow - i + 1];
1244 for (j = 1; j <= *ncol; ++j) {
1245 icstk[j + *ncol] = icol[*ncol - j + 1];
1260 ir1 = irstk[istk * *nrow + 1];
1261 ic1 = icstk[istk * *ncol + 1];
1268 } else if (ir1 < ic1) {
1289 irt = irstk[i + istk * *nrow];
1290 ict = icstk[j + istk * *ncol];
1299 f11act(&irstk[istk * *nrow + 1], &i, &nro,
1300 &irstk[(istk + 1) * *nrow + 1]);
1301 f11act(&icstk[istk * *ncol + 1], &j, &nco,
1302 &icstk[(istk + 1) * *ncol + 1]);
1303 } else if (irt > ict) {
1305 f11act(&icstk[istk * *ncol + 1], &j, &nco,
1306 &icstk[(istk + 1) * *ncol + 1]);
1308 f8xact(&irstk[istk * *nrow + 1], &ikh, &i,
1309 &nro, &irstk[(istk + 1) * *nrow + 1]);
1312 f11act(&irstk[istk * *nrow + 1], &i, &nro,
1313 &irstk[(istk + 1) * *nrow + 1]);
1315 f8xact(&icstk[istk * *ncol + 1], &ikh, &j,
1316 &nco, &icstk[(istk + 1) * *ncol + 1]);
1320 for (k = 1; k <= nco; ++k) {
1321 y += fact[icstk[k + (istk + 1) * *ncol]];
1326 for (k = 1; k <= nro; ++k) {
1327 y += fact[irstk[k + (istk + 1) * *nrow]];
1340 } while(1);/* end do */
1345 if (*dsp - amx <= *tol) {
1355 if (*dsp - amx <= *tol) {
1364 if (l > mstk[istk]) goto L100;
1371 if (irstk[l + istk * *nrow] <
1372 irstk[l - 1 + istk * *nrow]) goto L60;
1375 if (icstk[l + istk * *ncol] <
1376 icstk[l - 1 + istk * *ncol]) goto L60;
1382 -----------------------------------------------------------------------
1384 Purpose: Put node on stack in network algorithm.
1385 Usage: CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP,
1386 LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP,
1389 PASTP - The past path length. (Input)
1390 TOL - Tolerance for equivalence of past path lengths. (Input)
1391 KVAL - Key value. (Input)
1392 KEY - Vector of length LDKEY containing the key values. (in/out)
1393 LDKEY - Length of vector KEY. (Input)
1394 IPOIN - Vector of length LDKEY pointing to the
1395 linked list of past path lengths. (in/out)
1396 STP - Vector of length LSDTP containing the
1397 linked lists of past path lengths. (in/out)
1398 LDSTP - Length of vector STP. (Input)
1399 IFRQ - Vector of length LDSTP containing the past path
1400 frequencies. (in/out)
1401 NPOIN - Vector of length LDSTP containing the pointers to
1402 the next past path length. (in/out)
1403 NR - Vector of length LDSTP containing the right object
1404 pointers in the tree of past path lengths. (in/out)
1405 NL - Vector of length LDSTP containing the left object
1406 pointers in the tree of past path lengths. (in/out)
1407 IFREQ - Frequency of the current path length. (Input)
1408 ITOP - Pointer to the top of STP. (Input)
1409 IPSH - Option parameter. (Input)
1410 If IPSH is true, the past path length is found in the
1411 table KEY. Otherwise the location of the past path
1412 length is assumed known and to have been found in
1413 a previous call. ==>>>>> USING "static" variables
1414 -----------------------------------------------------------------------
1418 f5xact(double *pastp, const double *tol, int *kval, int *key, int *ldkey,
1419 int *ipoin, double *stp, int *ldstp, int *ifrq, int *npoin,
1420 int *nr, int *nl, int *ifreq, int *itop, int *ipsh)
1422 /* Local variables */
1423 static int itmp, ird, ipn, itp;
1424 double test1, test2;
1426 /* Parameter adjustments */
1437 /* Convert KVAL to int in range 1, ..., LDKEY. */
1438 ird = *kval % *ldkey + 1;
1439 /* Search for an unused location */
1440 for (itp = ird; itp <= *ldkey; ++itp) {
1441 if (key[itp] == *kval) {
1448 for (itp = 1; itp <= ird - 1; ++itp) {
1449 if (key[itp] == *kval) {
1456 /* Return if KEY array is full */
1458 prterr(6, "LDKEY is too small for this problem.\n"
1459 "It is not possible to estimate the value of LDKEY "
1461 "but twice the current value may be sufficient.");
1463 prterr(6, "LDKEY is too small for this problem.\n"
1464 "Try increasing the size of the workspace.");
1471 /* Return if STP array full */
1472 if (*itop > *ldstp) {
1474 prterr(7, "LDSTP is too small for this problem.\n"
1475 "It is not possible to estimate the value of LDSTP "
1477 "but twice the current value may be sufficient.");
1479 prterr(7, "LDSTP is too small for this problem.\n"
1480 "Try increasing the size of the workspace.");
1482 /* Update STP, etc. */
1486 stp[*itop] = *pastp;
1487 ifrq[*itop] = *ifreq;
1491 /* Find location, if any, of pastp */
1494 test1 = *pastp - *tol;
1495 test2 = *pastp + *tol;
1498 if (stp[ipn] < test1) {
1503 } else if (stp[ipn] > test2) {
1509 ifrq[ipn] += *ifreq;
1512 /* Return if STP array full */
1514 if (*itop > *ldstp) {
1516 prterr(7, "LDSTP is too small for this problem.\n"
1517 "It is not possible to estimate the value of LDSTP "
1519 "but twice the current value may be sufficient.");
1521 prterr(7, "LDSTP is too small for this problem.\n"
1522 "Try increasing the size of the workspace.");
1525 /* Find location to add value */
1530 if (stp[ipn] < test1) {
1538 } else if (stp[ipn] > test2) {
1547 /* Update STP, etc. */
1548 npoin[*itop] = npoin[itmp];
1549 npoin[itmp] = *itop;
1550 stp[*itop] = *pastp;
1551 ifrq[*itop] = *ifreq;
1557 -----------------------------------------------------------------------
1559 Purpose: Pop a node off the stack.
1560 Usage: CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, IPN)
1562 NROW - The number of rows in the table. (Input)
1563 IROW - Vector of length nrow containing the row sums on
1565 IFLAG - Set to 3 if there are no additional nodes to process.
1567 KYY - Constant mutlipliers used in forming the hash
1569 KEY - Vector of length LDKEY containing the hash table
1571 LDKEY - Length of vector KEY. (Input)
1572 LAST - Index of the last key popped off the stack. (In/out)
1573 IPN - Pointer to the linked list of past path lengths. (Output)
1574 -----------------------------------------------------------------------
1577 f6xact(int *nrow, int *irow, int *iflag, int *kyy, int *key, int
1578 *ldkey, int *last, int *ipn)
1582 /* Parameter adjustments */
1590 if (*last <= *ldkey) {
1591 if (key[*last] < 0) {
1594 /* Get KVAL from the stack */
1597 for (j = *nrow; j >= 2; --j) {
1598 irow[j] = kval / kyy[j];
1599 kval -= irow[j] * kyy[j];
1611 -----------------------------------------------------------------------
1613 Purpose: Generate the new nodes for given marinal totals.
1614 Usage: CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG)
1616 NROW - The number of rows in the table. (Input)
1617 IMAX - The row marginal totals. (Input)
1618 IDIF - The column counts for the new column. (in/out)
1619 K - Indicator for the row to decrement. (in/out)
1620 KS - Indicator for the row to increment. (in/out)
1621 IFLAG - Status indicator. (Output)
1622 If IFLAG is zero, a new table was generated. For
1623 IFLAG = 1, no additional tables could be generated.
1624 -----------------------------------------------------------------------
1628 f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
1634 /* Parameter adjustments */
1640 /* Find node which can be incremented, ks */
1644 } while (idif[*ks] == imax[*ks]);
1646 /* Find node to decrement (>ks) */
1647 if (idif[*k] > 0 && *k > *ks) {
1651 } while(imax[*k] == 0);
1655 /* Find node to increment (>=ks) */
1656 while (idif[m] >= imax[m]) {
1662 if (idif[m] == imax[m]) {
1669 /* Check for finish */
1670 for (k1 = *k + 1; k1 <= *nrow; ++k1) {
1679 /* Reallocate counts */
1681 for (i = 1; i <= *k; ++i) {
1689 m = min(mm, imax[*k]);
1692 } while (mm > 0 && *k != 1);
1694 /* Check that all counts reallocated */
1711 } while (idif[*ks] >= imax[*ks]);
1716 -----------------------------------------------------------------------
1718 Purpose: Routine for reducing a vector when there is a zero
1720 Usage: CALL F8XACT (IROW, IS, I1, IZERO, NEW)
1722 IROW - Vector containing the row counts. (Input)
1723 IS - Indicator. (Input)
1724 I1 - Indicator. (Input)
1725 IZERO - Position of the zero. (Input)
1726 NEW - Vector of new row counts. (Output)
1727 -----------------------------------------------------------------------
1731 f8xact(int *irow, int *is, int *i1, int *izero, int *myNew)
1735 /* Parameter adjustments */
1740 for (i = 1; i < *i1; ++i)
1743 for (i = *i1; i <= *izero - 1; ++i) {
1744 if (*is >= irow[i + 1])
1746 myNew[i] = irow[i + 1];
1760 -----------------------------------------------------------------------
1762 Purpose: Computes the log of a multinomial coefficient.
1763 Usage: F9XACT(N, MM, IR, FACT)
1765 N - Length of IR. (Input)
1766 MM - Number for factorial in numerator. (Input)
1767 IR - Vector of length N containing the numbers for
1768 the denominator of the factorial. (Input)
1769 FACT - Table of log factorials. (Input)
1770 F9XACT - The log of the multinomal coefficient. (Output)
1771 -----------------------------------------------------------------------
1775 f9xact(int *n, int *mm, int *ir, double *fact)
1781 for (k = 0; k < *n; k++)
1787 -----------------------------------------------------------------------
1789 Purpose: Computes the shortest path length for special tables.
1790 Usage: F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, NE, M)
1792 NROW - The number of rows in the table. (Input)
1793 IROW - Vector of length NROW containing the row totals. (Input)
1794 NCOL - The number of columns in the table. (Input)
1795 ICO - Vector of length NCOL containing the column totals.(Input)
1796 VAL - The shortest path. (Output)
1797 XMIN - Set to true if shortest path obtained. (Output)
1798 FACT - Vector containing the logarithms of factorials. (Input)
1799 ND - Workspace vector of length NROW. (Input)
1800 NE - Workspace vector of length NCOL. (Input)
1801 M - Workspace vector of length NCOL. (Input)
1803 Chapter: STAT/LIBRARY Categorical and Discrete Data Analysis
1804 -----------------------------------------------------------------------
1808 f10act(int *nrow, int *irow, int *ncol, int *icol, double *val,
1809 int *xmin, double *fact, int *nd, int *ne, int *m)
1811 /* Local variables */
1812 int i, is, ix, nrw1;
1814 /* Parameter adjustments */
1822 for (i = 1; i <= *nrow - 1; ++i)
1825 is = icol[1] / *nrow;
1826 ix = icol[1] - *nrow * is;
1832 for (i = 2; i <= *ncol; ++i) {
1833 ix = icol[i] / *nrow;
1836 ix = icol[i] - *nrow * ix;
1842 for (i = *nrow - 2; i >= 1; --i)
1847 for (i = *nrow; i >= 2; --i) {
1848 ix = ix + is + nd[nrw1 - i] - irow[i];
1853 for (i = 1; i <= *ncol; ++i) {
1856 *val = *val + is * fact[ix + 1] + (*nrow - is) * fact[ix];
1864 -----------------------------------------------------------------------
1866 Purpose: Routine for revising row totals.
1867 Usage: CALL F11ACT (IROW, I1, I2, NEW)
1869 IROW - Vector containing the row totals. (Input)
1870 I1 - Indicator. (Input)
1871 I2 - Indicator. (Input)
1872 NEW - Vector containing the row totals. (Output)
1873 -----------------------------------------------------------------------
1876 f11act(int *irow, int *i1, int *i2, int *myNew)
1880 /* Parameter adjustments */
1884 for (i = 1; i <= (*i1 - 1); ++i) myNew[i] = irow[i];
1885 for (i = *i1; i <= *i2; ++i) myNew[i] = irow[i + 1];
1891 -----------------------------------------------------------------------
1893 Purpose: Print an error message and stop.
1894 Usage: prterr(icode, mes)
1896 icode - Integer code for the error message. (Input)
1897 mes - Character string containing the error message. (Input)
1898 -----------------------------------------------------------------------
1901 prterr(int icode, char *mes)
1903 // PROBLEM "FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY);
1904 // printf("FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY));
1905 printf("FEXACT error %d.\n", icode); //xx:another one of my ugly kluges
1910 -----------------------------------------------------------------------
1912 Purpose: Routine for allocating workspace.
1913 Usage: iwork (iwkmax, iwkpt, number, itype)
1915 iwkmax - Maximum (int) amount of workspace. (Input)
1916 iwkpt - Amount of (int) workspace currently allocated. (in/out)
1917 number - Number of elements of workspace desired. (Input)
1918 itype - Workspace type. (Input)
1923 iwork(): Index in rwrk, dwrk, or iwrk of the beginning of
1924 the first free element in the workspace array. (Output)
1925 -----------------------------------------------------------------------
1928 iwork(int iwkmax, int *iwkpt, int number, int itype)
1933 if (itype == 2 || itype == 3)
1938 *iwkpt += (number << 1);
1941 if (*iwkpt > iwkmax)
1942 prterr(40, "Out of workspace.");
1949 void isort(int *n, int *ix)
1952 -----------------------------------------------------------------------
1954 Purpose: Shell sort for an int vector.
1955 Usage: CALL ISORT (N, IX)
1957 N - Lenth of vector IX. (Input)
1958 IX - Vector to be sorted. (in/out)
1959 -----------------------------------------------------------------------
1961 static int ikey, i, j, m, il[10], kl, it, iu[10], ku;
1963 /* Parameter adjustments */
1979 /* Find element in first half */
1983 if (ix[ikey] > ix[i]) {
1987 /* Find element in second half */
1990 if (ix[j] > ix[ikey]) {
2003 /* Save upper and lower subscripts of the array yet to be sorted */
2005 if (j - kl < ku - j) {
2019 prterr(20, "This should never occur.");
2021 /* Use another segment */
2032 double gammds(double *y, double *p, int *ifault)
2035 -----------------------------------------------------------------------
2037 Purpose: Cumulative distribution for the gamma distribution.
2038 Usage: PGAMMA (Q, ALPHA,IFAULT)
2040 Q - Value at which the distribution is desired. (Input)
2041 ALPHA - Parameter in the gamma distribution. (Input)
2042 IFAULT - Error indicator. (Output)
2045 1 An argument is misspecified.
2046 2 A numerical error has occurred.
2047 PGAMMA - The cdf for the gamma distribution with parameter alpha
2048 evaluated at Q. (Output)
2049 -----------------------------------------------------------------------
2051 Algorithm AS 147 APPL. Statist. (1980) VOL. 29, P. 113
2053 Computes the incomplete gamma integral for positive parameters Y, P
2054 using and infinite series.
2057 static double a, c, f, g;
2060 /* Checks for the admissibility of arguments and value of F */
2063 if (*y <= 0. || *p <= 0.) {
2069 ALOGAM is natural log of gamma function no need to test ifail as
2070 an error is impossible
2074 f = exp(*p * log(*y) - alogam(&a, &ifail) - *y);
2096 -----------------------------------------------------------------------
2098 Purpose: Value of the log-gamma function.
2099 Usage: ALOGAM (X, IFAULT)
2101 X - Value at which the log-gamma function is to be evaluated.
2103 IFAULT - Error indicator. (Output)
2107 ALGAMA - The value of the log-gamma function at XX. (Output)
2108 -----------------------------------------------------------------------
2110 Algorithm ACM 291, Comm. ACM. (1966) Vol. 9, P. 684
2112 Evaluates natural logarithm of gamma(x) for X greater than zero.
2115 double alogam(double *x, int *ifault)
2117 /* Initialized data */
2118 //printf("alogam x = %f\t%d\n",*x,*ifault);
2119 static double a1 = .918938533204673;
2120 static double a2 = 5.95238095238e-4;
2121 static double a3 = 7.93650793651e-4;
2122 static double a4 = .002777777777778;
2123 static double a5 = .083333333333333;
2125 /* Local variables */
2126 static double f, y, z;
2151 //printf("returning %f\n",(f + (y - .5) * log(y) - y + a1 + (((-a2 * z + a3) * z - a4) * z + a5) / y));
2152 return(f + (y - .5) * log(y) - y + a1 +
2153 (((-a2 * z + a3) * z - a4) * z + a5) / y);
2157 #endif /* not USING_R */