1 /*
2
3 This file is a C++ translation of the Mehta-Patel network
4 algorithm for computing association p-values on NxM contigency
5 table
6
7 Original source from
8
9 http://portal.acm.org/citation.cfm?id=214326&jmp=indexterms&dl=GUIDE&dl=ACM
10
11 Cleaned up and dapted for C++. Included a bug fix for negative keys
12 in hash tables in f3xact_.
13
14 Sergei L. Kosakovsky Pond,
15 October 2003
16
17 Function prototype is declared in matrix.h
18
19 */
20
21 #include <math.h>
22
23 #include "global_things.h"
24 #include "hy_strings.h"
25
26 using namespace hy_global;
27
28
29 /* ----------------------------------------------------------------------- */
30 // Function Prototypes
31
32 void allocate_fexact_keys (long, long);
33 void free_fexact_keys (void);
34
35 int fexact_(long, long, double *, double , double , double , double *, double *);
36 int isort_ (long*, long *);
37 double alogam_(double*, long *);
38 double gammds_(double*, double *, long *);
39 int f2xact_(long *, long *, double *,
40 double *, double *, double *,
41 double *, double *, double *, long *, long *, long *, long *, long *,
42 long *, long *, long *, double *, long *,
43 long *, double *, double *, double *,
44 long *, long *, double *);
45
46
47 int f3xact_(long *, long *, long *, long *, double *, long *, double *,
48 long *, long *, long *, long *, long *, long *, long *, long *, long *, double *, double *, double *);
49 int f4xact_(long *, long *, long *, long *, double *, double *, long *,
50 long *, long *, long *, long *, long *, long *, double *, double *);
51 int f5xact_(double *, double *, long *, long *, long *, long *,
52 double *, long *, long *, long *, long *, long *, long *, long *, bool *);
53 int f6xact_(long *, long *, long *, long *, long *, long *, long *, long *);
54 int f7xact_(long *, long *, long *, long *, long *, long *);
55 int f8xact_(long *, long *, long *, long *, long *);
56 double f9xact_(long *, long *, long *, double *);
57 int f10act_(long, long *, long , long * , double * , bool * , double * , long * , long * , long *);
58 void f11act_(long *, long, long, long *);
59 long i_dnnt (double *);
60 /* ----------------------------------------------------------------------- */
61
62 #define Minimum(a,b) (((a)<(b))?(a):(b))
63 #define Maximum(a,b) (((a)>(b))?(a):(b))
64
65 /* ----------------------------------------------------------------------- */
i_dnnt(double * x)66 inline long i_dnnt(double *x)
67 {
68 return( (*x)>=0 ? (long) (*x + .5) : (long) (*x - .5) );
69 }
70
71 long *fexact_i4 = nil,
72 *fexact_i5 = nil,
73 *fexact_i7 = nil,
74 *fexact_i10 = nil;
75
76 double *fexact_i6 = nil,
77 *fexact_i8 = nil,
78 *fexact_i9 = nil,
79 *fexact_i9a = nil;
80
81 long fexact_ldkey = 0,
82 fexact_ldstp = 0;
83
84 /* ----------------------------------------------------------------------- */
85
allocate_fexact_keys(long ldkey,long mult)86 void allocate_fexact_keys (long ldkey, long mult)
87 {
88 fexact_ldkey = ldkey;
89 fexact_ldstp = ldkey * mult;
90 long i__1 = ldkey << 1;
91
92 fexact_i4 = (long*) MemAllocate(i__1*sizeof (long));
93 fexact_i5 = (long*) MemAllocate(i__1*sizeof (long));
94 i__1 = fexact_ldstp << 1;
95 fexact_i6 = (double*) MemAllocate(i__1*sizeof (double));
96 i__1 = fexact_ldstp * 6;
97 fexact_i7 = (long*) MemAllocate(i__1*sizeof (long));
98 i__1 = ldkey << 1;
99 fexact_i8 = (double*) MemAllocate(i__1*sizeof (double));
100 fexact_i9 = (double*) MemAllocate(i__1*sizeof (double));
101 fexact_i9a = (double*) MemAllocate(i__1*sizeof (double));
102 fexact_i10 = (long*) MemAllocate(i__1*sizeof (long));
103 }
104
105 /* ----------------------------------------------------------------------- */
106
free_fexact_keys(void)107 void free_fexact_keys (void)
108 {
109 free (fexact_i4);
110 free (fexact_i5);
111 free (fexact_i6);
112 free (fexact_i7);
113 free (fexact_i8);
114 free (fexact_i9);
115 free (fexact_i9a);
116 free (fexact_i10);
117
118 fexact_i4 = nil;
119 }
120
121
122
123 /* ----------------------------------------------------------------------- */
124 /* Name: isort_ */
125
126 /* Purpose: Shell sort for an long vector. */
127
128 /* Usage: isort_ (N, IX) */
129
130 /* Arguments: */
131 /* n - Lenth of vector IX. (Input) */
132 /* ix - Vector to be sorted. (Input/output) */
133 /* ----------------------------------------------------------------------- */
134
isort_(long * n,long * ix)135 int isort_(long *n, long *ix)
136 {
137 long i__,
138 j,
139 m,
140 il[10],
141 kl,
142 it,
143 iu[10],
144 ku,
145 ikey;
146 /* Parameter adjustments */
147 --ix;
148
149 /* Function Body */
150 m = 1;
151 i__ = 1;
152 j = *n;
153 L10:
154 if (i__ >= j) {
155 goto L40;
156 }
157
158 kl = i__;
159 ku = j;
160 ikey = i__;
161 ++j;
162 /* Find element in first half */
163 L20:
164 ++i__;
165 if (i__ < j) {
166 if (ix[ikey] > ix[i__]) {
167 goto L20;
168 }
169 }
170 /* Find element in second half */
171 L30:
172 --j;
173 if (ix[j] > ix[ikey]) {
174 goto L30;
175 }
176 /* Interchange */
177 if (i__ < j) {
178 it = ix[i__];
179 ix[i__] = ix[j];
180 ix[j] = it;
181 goto L20;
182 }
183 it = ix[ikey];
184 ix[ikey] = ix[j];
185 ix[j] = it;
186 /* Save upper and lower subscripts of */
187 /* the array yet to be sorted */
188 if (m < 11) {
189 if (j - kl < ku - j) {
190 il[m - 1] = j + 1;
191 iu[m - 1] = ku;
192 i__ = kl;
193 --j;
194 } else {
195 il[m - 1] = kl;
196 iu[m - 1] = j - 1;
197 i__ = j + 1;
198 j = ku;
199 }
200 ++m;
201 goto L10;
202 } else {
203 HandleApplicationError ("Internal error in shell sort");
204 }
205 /* Use another segment */
206 L40:
207 --m;
208 if (m == 0) {
209 goto L9000;
210 }
211 i__ = il[m - 1];
212 j = iu[m - 1];
213 goto L10;
214
215 L9000:
216 return 0;
217 } /* isort_ */
218
219 /* ----------------------------------------------------------------------- */
220 /* Name: GAMMDS */
221
222 /* Purpose: Cumulative distribution for the gamma distribution. */
223
224 /* Usage: PGAMMA (Q, ALPHA,IFAULT) */
225
226 /* Arguments: */
227 /* Q - Value at which the distribution is desired. (Input) */
228 /* ALPHA - Parameter in the gamma distribution. (Input) */
229 /* IFAULT - Error indicator. (Output) */
230 /* IFAULT DEFINITION */
231 /* 0 No error */
232 /* 1 An argument is misspecified. */
233 /* 2 A numerical error has occurred. */
234 /* PGAMMA - The cdf for the gamma distribution with parameter alpha */
235 /* evaluated at Q. (Output) */
236 /* ----------------------------------------------------------------------- */
237
gammds_(double * y,double * p,long * ifault)238 double gammds_(double *y, double *p, long *ifault)
239 {
240 /* Initialized data */
241
242 double e = 1e-6;
243 double zero = 0.;
244 double one = 1.;
245
246 /* System generated locals */
247 double ret_val, d__1, d__2;
248
249
250 /* Local variables */
251 double a, c__, f;
252 long ifail;
253
254 *ifault = 1;
255 ret_val = zero;
256 if (*y <= zero || *p <= zero) {
257 return ret_val;
258 }
259 *ifault = 2;
260
261 d__2 = *p + one;
262 d__1 = *p * log(*y) - alogam_(&d__2, &ifail) - *y;
263 f = exp(d__1);
264 if (f == zero) {
265 return ret_val;
266 }
267 *ifault = 0;
268
269 /* Series begins */
270
271 c__ = one;
272 ret_val = one;
273 a = *p;
274 L10:
275 a += one;
276 c__ = c__ * *y / a;
277 ret_val += c__;
278 if (c__ / ret_val > e) {
279 goto L10;
280 }
281 ret_val *= f;
282 return ret_val;
283 } /* gammds_ */
284
285 /* ----------------------------------------------------------------------- */
286 /* Name: ALOGAM */
287
288 /* Purpose: Value of the log-gamma function. */
289
290 /* Usage: ALOGAM (X, IFAULT) */
291
292 /* Arguments: */
293 /* X - Value at which the log-gamma function is to be evaluated. */
294 /* (Input) */
295 /* IFAULT - Error indicator. (Output) */
296 /* IFAULT DEFINITION */
297 /* 0 No error */
298 /* 1 X .LT. 0 */
299 /* ALGAMA - The value of the log-gamma function at XX. (Output) */
300 /* ----------------------------------------------------------------------- */
301
alogam_(double * x,long * ifault)302 double alogam_(double *x, long *ifault)
303 {
304 /* Initialized data */
305
306 static double a1 = .918938533204673;
307 static double a2 = 5.95238095238e-4;
308 static double a3 = 7.93650793651e-4;
309 static double a4 = .002777777777778;
310 static double a5 = .083333333333333;
311 static double half = .5;
312 static double zero = 0.;
313 static double one = 1.;
314 static double seven = 7.;
315
316 /* System generated locals */
317 double ret_val;
318
319 /* Local variables */
320 double f, y, z__;
321
322 ret_val = zero;
323 *ifault = 1;
324 if (*x < zero) {
325 return ret_val;
326 }
327 *ifault = 0;
328 y = *x;
329 f = zero;
330 if (y >= seven) {
331 goto L30;
332 }
333 f = y;
334 L10:
335 y += one;
336 if (y >= seven) {
337 goto L20;
338 }
339 f *= y;
340 goto L10;
341 L20:
342 f = -log(f);
343 L30:
344 z__ = one / (y * y);
345 ret_val = f + (y - half) * log(y) - y + a1 + (((-a2 * z__ + a3) * z__ -
346 a4) * z__ + a5) / y;
347 return ret_val;
348 } /* alogam_ */
349
350 /* ----------------------------------------------------------------------- */
351 /* Name: F11ACT */
352
353 /* Purpose: Routine for revising row totals. */
354
355 /* Arguments: */
356 /* IROW - Vector containing the row totals. (Input) */
357 /* I1 - Indicator. (Input) */
358 /* I2 - Indicator. (Input) */
359 /* N - Vector containing the row totals. (Input) */
360 /* ----------------------------------------------------------------------- */
f11act_(long * irow,long i1,long i2,long * n)361 void f11act_(long *irow, long i1, long i2, long *n)
362 {
363 for (long i = 0; i < i1-1; i++) {
364 n[i] = irow[i];
365 }
366 {
367 for (long i = i1-1; i < i2; i++) {
368 n[i] = irow[i + 1];
369 }
370 }
371 } /* f11act_ */
372
373 /* ----------------------------------------------------------------------- */
374 /* Name: F10ACT */
375
376 /* Purpose: Computes the shortest path length for special tables. */
377
378 /* Usage: CALL F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, */
379 /* NE, M) */
380
381 /* Arguments: */
382 /* NROW - The number of rows in the table. (Input) */
383 /* IROW - Vector of length NROW containing the row totals. (Input) */
384 /* NCOL - The number of columns in the table. (Input) */
385 /* ICO - Vector of length NCOL containing the column totals. */
386 /* (Input) */
387 /* VAL - The shortest path. (Output) */
388 /* XMIN - Set to true if shortest path obtained. (Output) */
389 /* FACT - Vector containing the logarithms of factorials. */
390 /* (Input) */
391 /* ND - Workspace vector of length NROW. */
392 /* NE - Workspace vector of length NCOL. */
393 /* M - Workspace vector of length NCOL. */
394
395 /* Chapter: STAT/LIBRARY Categorical and Discrete Data Analysis */
396 /* ----------------------------------------------------------------------- */
f10act_(long nrow,long * irow,long ncol,long * icol,double * val,bool * xmin,double * fact,long * nd,long * ne,long * m)397 int f10act_(long nrow, long *irow, long ncol, long *icol, double *val, bool *xmin, double *fact, long *nd, long *ne, long *m)
398 {
399 /* System generated locals */
400 long is, ix;
401
402 /* Parameter adjustments */
403
404 long i__,
405 i__1;
406
407 --m;
408 --ne;
409 --nd;
410 --icol;
411 --irow;
412
413 /* Function Body */
414
415 //for (long i1 = 0; i1 < nrow; i1++)
416 // nd[i1] = 0;
417
418 i__1 = nrow - 1;
419 for (i__ = 1; i__ <= i__1; ++i__) {
420 nd[i__] = 0;
421 }
422
423 is = icol[1] / nrow;
424 ne[1] = is;
425
426 //is = icol[0]/nrow;
427 //ne[0] = is;
428
429 ix = icol[1] - nrow * is;
430 m[1] = ix;
431
432 //ix = icol[0] - nrow*is;
433 //m[0] = ix;
434
435 //if (ix != 0)
436 //++nd[ix-1];
437
438 if (ix != 0) {
439 ++nd[ix];
440 }
441
442 i__1 = ncol;
443 for (i__ = 2; i__ <= i__1; ++i__) {
444 ix = icol[i__] / nrow;
445 ne[i__] = ix;
446 is += ix;
447 ix = icol[i__] - nrow * ix;
448 m[i__] = ix;
449 if (ix != 0) {
450 ++nd[ix];
451 }
452 }
453
454 /*for (long i2 = 1; i2 < ncol; i2++)
455 {
456 long t = icol[i2]/nrow;
457 ne[i2] = t;
458 is += t;
459 t = icol[i2] - nrow*t;
460 m[i2] = t;
461 if (t != 0)
462 nd[t-1]++;
463 }*/
464
465 for (i__ = nrow - 2; i__ >= 1; --i__) {
466 nd[i__] += nd[i__ + 1];
467 }
468
469 //for (long i3 = nrow-3; i3>=0; i3--)
470 // nd[i3] += nd [i3+1];
471
472 ix = 0;
473 long nrw1 = nrow + 1;
474 for (i__ = nrow; i__ >= 2; --i__) {
475 ix = ix + is + nd[nrw1 - i__] - irow[i__];
476 if (ix < 0) {
477 return 0;
478 }
479 }
480
481 /*ix = 0;
482 for (long i4 = nrow; i4>=1; i4--)
483 {
484 ix += is + nd[nrow-i4+1] - irow[i4];
485 if (ix<0)
486 return 0;
487 }*/
488
489 i__1 = ncol;
490 for (i__ = 1; i__ <= i__1; ++i__) {
491 ix = ne[i__];
492 is = m[i__];
493 *val = *val + is * fact[ix + 1] + (nrow - is) * fact[ix];
494 }
495 *xmin = true;
496
497 /* double temp = 0.0;
498
499 for (long i5 = 0; i5 < ncol; i5++)
500 {
501 ix = ne[i5];
502 is = m[i5];
503
504 temp += is *fact[ix+1] + (nrow-is)*fact[ix];
505 }
506 *val += temp;
507 *xmin = true;*/
508
509 return 0;
510 }
511
512 /* ----------------------------------------------------------------------- */
513 /* Name: F9XACT */
514
515 /* Purpose: Computes the log of a multinomial coefficient. */
516
517 /* Usage: F9XACT(N, MM, IR, FACT) */
518
519 /* Arguments: */
520 /* N - Length of IR. (Input) */
521 /* MM - Number for factorial in numerator. (Input) */
522 /* IR - Vector of length N containing the numebers for the */
523 /* denominator of the factorial. (Input) */
524 /* FACT - Table of log factorials. (Input) */
525 /* F9XACT - The log of the multinomal coefficient. (Output) */
526 /* ----------------------------------------------------------------------- */
f9xact_(long * n,long * mm,long * ir,double * fact)527 double f9xact_(long *n, long *mm, long *ir, double *fact)
528 {
529 /* System generated locals */
530 long i__1,
531 k;
532
533 double ret_val;
534
535 /* Parameter adjustments */
536 --ir;
537
538 /* Function Body */
539 ret_val = fact[*mm];
540 i__1 = *n;
541 for (k = 1; k <= i__1; ++k) {
542 ret_val -= fact[ir[k]];
543 }
544
545 return ret_val;
546 } /* f9xact_ */
547
548 /* ----------------------------------------------------------------------- */
549 /* Name: F8XACT */
550
551 /* Purpose: Routine for reducing a vector when there is a zero */
552 /* element. */
553
554 /* Usage: CALL F8XACT (IROW, IS, I1, IZERO, NEW) */
555
556 /* Arguments: */
557 /* IROW - Vector containing the row counts. (Input) */
558 /* IS - Indicator. (Input) */
559 /* I1 - Indicator. (Input) */
560 /* IZERO - Position of the zero. (Input) */
561 /* NEW - Vector of new row counts. (Output) */
562 /* ----------------------------------------------------------------------- */
f8xact_(long * irow,long * is,long * i1,long * izero,long * new__)563 int f8xact_(long *irow, long *is, long *i1, long *izero, long *new__)
564 {
565 /* System generated locals */
566 long i__1,
567 i__;
568
569 /* Parameter adjustments */
570 --new__;
571 --irow;
572
573 /* Function Body */
574 i__1 = *i1 - 1;
575 for (i__ = 1; i__ <= i__1; ++i__) {
576 new__[i__] = irow[i__];
577 /* L10: */
578 }
579
580 i__1 = *izero - 1;
581 for (i__ = *i1; i__ <= i__1; ++i__) {
582 if (*is >= irow[i__ + 1]) {
583 goto L30;
584 }
585 new__[i__] = irow[i__ + 1];
586 /* L20: */
587 }
588
589 i__ = *izero;
590 L30:
591 new__[i__] = *is;
592 L40:
593 ++i__;
594 if (i__ > *izero) {
595 return 0;
596 }
597 new__[i__] = irow[i__];
598 goto L40;
599 } /* f8xact_ */
600
601 /* ----------------------------------------------------------------------- */
602 /* Name: F7XACT */
603
604 /* Purpose: Generate the new nodes for given marinal totals. */
605
606 /* Usage: CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG) */
607
608 /* Arguments: */
609 /* NROW - The number of rows in the table. (Input) */
610 /* IMAX - The row marginal totals. (Input) */
611 /* IDIF - The column counts for the new column. (Input/output) */
612 /* K - Indicator for the row to decrement. (Input/output) */
613 /* KS - Indicator for the row to increment. (Input/output) */
614 /* IFLAG - Status indicator. (Output) */
615 /* If IFLAG is zero, a new table was generated. For */
616 /* IFLAG = 1, no additional tables could be generated. */
617 /* ----------------------------------------------------------------------- */
f7xact_(long * nrow,long * imax,long * idif,long * k,long * ks,long * iflag)618 int f7xact_(long *nrow, long *imax, long *idif, long *k, long *ks, long *iflag)
619 {
620 /* System generated locals */
621 long i__1, i__2, i__, m, k1, mm;
622
623 /* Parameter adjustments */
624 --idif;
625 --imax;
626
627 /* Function Body */
628 *iflag = 0;
629 /* Find node which can be */
630 /* incremented, ks */
631 if (*ks == 0) {
632 L10:
633 ++(*ks);
634 if (idif[*ks] == imax[*ks]) {
635 goto L10;
636 }
637 }
638 /* Find node to decrement (>ks) */
639 /* L20: */
640 if (idif[*k] > 0 && *k > *ks) {
641 --idif[*k];
642 L30:
643 --(*k);
644 if (imax[*k] == 0) {
645 goto L30;
646 }
647 m = *k;
648 /* Find node to increment (>=ks) */
649 L40:
650 if (idif[m] >= imax[m]) {
651 --m;
652 goto L40;
653 }
654 ++idif[m];
655 /* Change ks */
656 if (m == *ks) {
657 if (idif[m] == imax[m]) {
658 *ks = *k;
659 }
660 }
661 } else {
662 /* Check for finish */
663 L50:
664 i__1 = *nrow;
665 for (k1 = *k + 1; k1 <= i__1; ++k1) {
666 if (idif[k1] > 0) {
667 goto L70;
668 }
669 /* L60: */
670 }
671 *iflag = 1;
672 goto L9000;
673 /* Reallocate counts */
674 L70:
675 mm = 1;
676 i__1 = *k;
677 for (i__ = 1; i__ <= i__1; ++i__) {
678 mm += idif[i__];
679 idif[i__] = 0;
680 /* L80: */
681 }
682 *k = k1;
683 L90:
684 --(*k);
685 /* Computing MIN */
686 i__1 = mm, i__2 = imax[*k];
687 m = Minimum(i__1,i__2);
688 idif[*k] = m;
689 mm -= m;
690 if (mm > 0 && *k != 1) {
691 goto L90;
692 }
693 /* Check that all counts */
694 /* reallocated */
695 if (mm > 0) {
696 if (k1 != *nrow) {
697 *k = k1;
698 goto L50;
699 }
700 *iflag = 1;
701 goto L9000;
702 }
703 /* Get ks */
704 --idif[k1];
705 *ks = 0;
706 L100:
707 ++(*ks);
708 if (*ks > *k) {
709 goto L9000;
710 }
711 if (idif[*ks] >= imax[*ks]) {
712 goto L100;
713 }
714 }
715
716 L9000:
717 return 0;
718 } /* f7xact_ */
719
720
721 /* ----------------------------------------------------------------------- */
722 /* Name: F6XACT */
723
724 /* Purpose: Pop a node off the stack. */
725
726 /* Usage: CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, */
727 /* IPN) */
728
729 /* Arguments: */
730 /* NROW - The number of rows in the table. (Input) */
731 /* IROW - Vector of length nrow containing the row sums on output. */
732 /* (Output) */
733 /* IFLAG - Set to 3 if there are no additional nodes to process. */
734 /* (Output) */
735 /* KYY - Constant mutlipliers used in forming the hash table key. */
736 /* (Input) */
737 /* KEY - Vector of length LDKEY containing the hash table keys. */
738 /* (Input/output) */
739 /* LDKEY - Length of vector KEY. (Input) */
740 /* LAST - Index of the last key popped off the stack. */
741 /* (Input/output) */
742 /* IPN - Pointer to the linked list of past path lengths. */
743 /* (Output) */
744 /* ----------------------------------------------------------------------- */
f6xact_(long * nrow,long * irow,long * iflag,long * kyy,long * key,long * ldkey,long * last,long * ipn)745 int f6xact_(long *nrow, long *irow, long *iflag, long *kyy, long *key, long *ldkey, long *last, long *ipn)
746 {
747 long j,
748 kval;
749
750
751 /* Parameter adjustments */
752 --key;
753 --kyy;
754 --irow;
755
756 /* Function Body */
757 L10:
758 ++(*last);
759 if (*last <= *ldkey) {
760 if (key[*last] < 0) {
761 goto L10;
762 }
763 /* Get KVAL from the stack */
764 kval = key[*last];
765 key[*last] = -9999;
766 for (j = *nrow; j >= 2; --j) {
767 irow[j] = kval / kyy[j];
768 kval -= irow[j] * kyy[j];
769 /* L20: */
770 }
771 irow[1] = kval;
772 *ipn = *last;
773 } else {
774 *last = 0;
775 *iflag = 3;
776 }
777 return 0;
778 } /* f6xact_ */
779
780 /* ----------------------------------------------------------------------- */
781 /* Name: F5XACT */
782
783 /* Purpose: Put node on stack in network algorithm. */
784
785 /* Usage: CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP, */
786 /* LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP, */
787 /* IPSH) */
788
789 /* Arguments: */
790 /* PASTP - The past path length. (Input) */
791 /* TOL - Tolerance for equivalence of past path lengths. (Input) */
792 /* KVAL - Key value. (Input) */
793 /* KEY - Vector of length LDKEY containing the key values. */
794 /* (Input/output) */
795 /* LDKEY - Length of vector KEY. (Input) */
796 /* IPOIN - Vector of length LDKEY pointing to the linked list */
797 /* of past path lengths. (Input/output) */
798 /* STP - Vector of length LSDTP containing the linked lists */
799 /* of past path lengths. (Input/output) */
800 /* LDSTP - Length of vector STP. (Input) */
801 /* IFRQ - Vector of length LDSTP containing the past path */
802 /* frequencies. (Input/output) */
803 /* NPOIN - Vector of length LDSTP containing the pointers to */
804 /* the next past path length. (Input/output) */
805 /* NR - Vector of length LDSTP containing the right object */
806 /* pointers in the tree of past path lengths. */
807 /* (Input/output) */
808 /* NL - Vector of length LDSTP containing the left object */
809 /* pointers in the tree of past path lengths. */
810 /* (Input/output) */
811 /* IFREQ - Frequency of the current path length. (Input) */
812 /* ITOP - Pointer to the top of STP. (Input) */
813 /* IPSH - Option parameter. (Input) */
814 /* If IPSH is true, the past path length is found in the */
815 /* table KEY. Otherwise the location of the past path */
816 /* length is assumed known and to have been found in */
817 /* a previous call. */
818 /* ----------------------------------------------------------------------- */
f5xact_(double * pastp,double * tol,long * kval,long * key,long * ldkey,long * ipoin,double * stp,long * ldstp,long * ifrq,long * npoin,long * nr,long * nl,long * ifreq,long * itop,bool * ipsh)819 int f5xact_(double *pastp, double *tol, long *kval, long *key, long *ldkey, long *ipoin,
820 double *stp, long *ldstp, long *ifrq, long *npoin, long *nr, long *nl, long *ifreq, long *itop, bool *ipsh)
821 {
822 /* System generated locals */
823 static long itp;
824
825 long i__1,
826 ird,
827 ipn,
828 itmp;
829
830 double test1,
831 test2;
832
833 _String errMsg;
834
835 /* Parameter adjustments */
836 --nl;
837 --nr;
838 --npoin;
839 --ifrq;
840 --stp;
841 --ipoin;
842 --key;
843
844 /* Function Body */
845 if (*ipsh) {
846 /* Convert KVAL to long in range */
847 /* 1, ..., LDKEY. */
848 ird = *kval % *ldkey + 1;
849 /* Search for an unused location */
850 i__1 = *ldkey;
851 for (itp = ird; itp <= i__1; ++itp) {
852 if (key[itp] == *kval) {
853 goto L40;
854 }
855 if (key[itp] < 0) {
856 goto L30;
857 }
858 /* L10: */
859 }
860
861 i__1 = ird - 1;
862 for (itp = 1; itp <= i__1; ++itp) {
863 if (key[itp] == *kval) {
864 goto L40;
865 }
866 if (key[itp] < 0) {
867 goto L30;
868 }
869 /* L20: */
870 }
871 /* Return if KEY array is full */
872 HandleApplicationError("Fisher Exact:LDKEY is too small for this problem. It is not possible to estimate the value of LDKEY required, but twice the current value may be sufficient.");
873 return 0;
874 /* Update KEY */
875 L30:
876 key[itp] = *kval;
877 ++(*itop);
878 ipoin[itp] = *itop;
879 /* Return if STP array full */
880 if (*itop > *ldstp) {
881 HandleApplicationError("Fisher Exact: LDSTP is too small for this problem. It is not possible to estimate the value of LDSTP required, but twice the current value may be sufficient.");
882 return 0;
883 }
884 /* Update STP, etc. */
885 npoin[*itop] = -1;
886 nr[*itop] = -1;
887 nl[*itop] = -1;
888 stp[*itop] = *pastp;
889 ifrq[*itop] = *ifreq;
890 goto L9000;
891 }
892 /* Find location, if any, of pastp */
893 L40:
894 ipn = ipoin[itp];
895 test1 = *pastp - *tol;
896 test2 = *pastp + *tol;
897
898 L50:
899 if (stp[ipn] < test1) {
900 ipn = nl[ipn];
901 if (ipn > 0) {
902 goto L50;
903 }
904 } else if (stp[ipn] > test2) {
905 ipn = nr[ipn];
906 if (ipn > 0) {
907 goto L50;
908 }
909 } else {
910 ifrq[ipn] += *ifreq;
911 goto L9000;
912 }
913 /* Return if STP array full */
914 ++(*itop);
915 if (*itop > *ldstp) {
916 HandleApplicationError ("Fisher Exact: LDSTP is too small for this problem. It is not possible to estimate the value of LDSTP required, but twice the current value may be sufficient.");
917 goto L9000;
918 }
919 /* Find location to add value */
920 ipn = ipoin[itp];
921 itmp = ipn;
922 L60:
923 if (stp[ipn] < test1) {
924 itmp = ipn;
925 ipn = nl[ipn];
926 if (ipn > 0) {
927 goto L60;
928 } else {
929 nl[itmp] = *itop;
930 }
931 } else if (stp[ipn] > test2) {
932 itmp = ipn;
933 ipn = nr[ipn];
934 if (ipn > 0) {
935 goto L60;
936 } else {
937 nr[itmp] = *itop;
938 }
939 }
940 /* Update STP, etc. */
941 npoin[*itop] = npoin[itmp];
942 npoin[itmp] = *itop;
943 stp[*itop] = *pastp;
944 ifrq[*itop] = *ifreq;
945 nl[*itop] = -1;
946 nr[*itop] = -1;
947
948 L9000:
949 return 0;
950 } /* f5xact_ */
951
952 /* ----------------------------------------------------------------------- */
953 /* Name: F4XACT */
954
955 /* Purpose: Computes the longest path length for a given table. */
956
957 /* Usage: CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK, */
958 /* NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK, */
959 /* TOL) */
960
961 /* Arguments: */
962 /* NROW - The number of rows in the table. (Input) */
963 /* IROW - Vector of length NROW containing the row sums for the */
964 /* table. (Input) */
965 /* NCOL - The number of columns in the table. (Input) */
966 /* ICOL - Vector of length K containing the column sums for the */
967 /* table. (Input) */
968 /* DSP - The shortest path for the table. (Output) */
969 /* FACT - Vector containing the logarithms of factorials. (Input) */
970 /* ICSTK - NCOL by NROW+NCOL+1 work array. */
971 /* NCSTK - Work vector of length NROW+NCOL+1. */
972 /* LSTK - Work vector of length NROW+NCOL+1. */
973 /* MSTK - Work vector of length NROW+NCOL+1. */
974 /* NSTK - Work vector of length NROW+NCOL+1. */
975 /* NRSTK - Work vector of length NROW+NCOL+1. */
976 /* IRSTK - NROW by MAX(NROW,NCOL) work array. */
977 /* YSTK - Work vector of length NROW+NCOL+1. */
978 /* TOL - Tolerance. (Input) */
979 /* ----------------------------------------------------------------------- */
f4xact_(long * nrow,long * irow,long * ncol,long * icol,double * dsp,double * fact,long * icstk,long * ncstk,long * lstk,long * mstk,long * nstk,long * nrstk,long * irstk,double * ystk,double * tol)980 int f4xact_(long *nrow, long *irow, long *ncol,
981 long *icol, double *dsp, double *fact, long *icstk,
982 long *ncstk, long *lstk, long *mstk, long *nstk,
983 long *nrstk, long *irstk, double *ystk, double *tol)
984 {
985 /* System generated locals */
986 long icstk_dim1,
987 icstk_offset,
988 irstk_dim1,
989 irstk_offset,
990 i__1,
991 i__,
992 j,
993 k,
994 l,
995 m,
996 n,
997 mn,
998 ic1,
999 ir1,
1000 ict,
1001 nco,
1002 irt,
1003 nro,
1004 istk;
1005
1006 double y,
1007 amx;
1008
1009
1010 /* Parameter adjustments */
1011 irstk_dim1 = *nrow;
1012 irstk_offset = 1 + irstk_dim1;
1013 irstk -= irstk_offset;
1014 --irow;
1015 icstk_dim1 = *ncol;
1016 icstk_offset = 1 + icstk_dim1;
1017 icstk -= icstk_offset;
1018 --icol;
1019 --ncstk;
1020 --lstk;
1021 --mstk;
1022 --nstk;
1023 --nrstk;
1024 --ystk;
1025
1026 /* Function Body */
1027 if (*nrow == 1) {
1028 i__1 = *ncol;
1029 for (i__ = 1; i__ <= i__1; ++i__) {
1030 *dsp -= fact[icol[i__]];
1031 /* L10: */
1032 }
1033 goto L9000;
1034 }
1035
1036 if (*ncol == 1) {
1037 i__1 = *nrow;
1038 for (i__ = 1; i__ <= i__1; ++i__) {
1039 *dsp -= fact[irow[i__]];
1040 /* L20: */
1041 }
1042 goto L9000;
1043 }
1044
1045 if (*nrow * *ncol == 4) {
1046 if (irow[2] <= icol[2]) {
1047 *dsp = *dsp - fact[irow[2]] - fact[icol[1]] - fact[icol[2] - irow[
1048 2]];
1049 } else {
1050 *dsp = *dsp - fact[icol[2]] - fact[irow[1]] - fact[irow[2] - icol[
1051 2]];
1052 }
1053 goto L9000;
1054 }
1055 /* initialization before loop */
1056 i__1 = *nrow;
1057 for (i__ = 1; i__ <= i__1; ++i__) {
1058 irstk[i__ + irstk_dim1] = irow[*nrow - i__ + 1];
1059 /* L30: */
1060 }
1061
1062 i__1 = *ncol;
1063 for (j = 1; j <= i__1; ++j) {
1064 icstk[j + icstk_dim1] = icol[*ncol - j + 1];
1065 /* L40: */
1066 }
1067
1068 nro = *nrow;
1069 nco = *ncol;
1070 nrstk[1] = nro;
1071 ncstk[1] = nco;
1072 ystk[1] = 0.f;
1073 y = 0.f;
1074 istk = 1;
1075 l = 1;
1076 amx = 0.f;
1077
1078 L50:
1079 ir1 = irstk[istk * irstk_dim1 + 1];
1080 ic1 = icstk[istk * icstk_dim1 + 1];
1081 if (ir1 > ic1) {
1082 if (nro >= nco) {
1083 m = nco - 1;
1084 n = 2;
1085 } else {
1086 m = nro;
1087 n = 1;
1088 }
1089 } else if (ir1 < ic1) {
1090 if (nro <= nco) {
1091 m = nro - 1;
1092 n = 1;
1093 } else {
1094 m = nco;
1095 n = 2;
1096 }
1097 } else {
1098 if (nro <= nco) {
1099 m = nro - 1;
1100 n = 1;
1101 } else {
1102 m = nco - 1;
1103 n = 2;
1104 }
1105 }
1106
1107 L60:
1108 if (n == 1) {
1109 i__ = l;
1110 j = 1;
1111 } else {
1112 i__ = 1;
1113 j = l;
1114 }
1115
1116 irt = irstk[i__ + istk * irstk_dim1];
1117 ict = icstk[j + istk * icstk_dim1];
1118 mn = irt;
1119 if (mn > ict) {
1120 mn = ict;
1121 }
1122 y += fact[mn];
1123 if (irt == ict) {
1124 --nro;
1125 --nco;
1126 f11act_(&irstk[istk * irstk_dim1 + 1], i__, nro, &irstk[(istk + 1) *
1127 irstk_dim1 + 1]);
1128 f11act_(&icstk[istk * icstk_dim1 + 1], j, nco, &icstk[(istk + 1) *
1129 icstk_dim1 + 1]);
1130 } else if (irt > ict) {
1131 --nco;
1132 f11act_(&icstk[istk * icstk_dim1 + 1], j, nco, &icstk[(istk + 1) *
1133 icstk_dim1 + 1]);
1134 i__1 = irt - ict;
1135 f8xact_(&irstk[istk * irstk_dim1 + 1], &i__1, &i__, &nro, &irstk[(
1136 istk + 1) * irstk_dim1 + 1]);
1137 } else {
1138 --nro;
1139 f11act_(&irstk[istk * irstk_dim1 + 1], i__, nro, &irstk[(istk + 1) *
1140 irstk_dim1 + 1]);
1141 i__1 = ict - irt;
1142 f8xact_(&icstk[istk * icstk_dim1 + 1], &i__1, &j, &nco, &icstk[(istk
1143 + 1) * icstk_dim1 + 1]);
1144 }
1145
1146 if (nro == 1) {
1147 i__1 = nco;
1148 for (k = 1; k <= i__1; ++k) {
1149 y += fact[icstk[k + (istk + 1) * icstk_dim1]];
1150 /* L70: */
1151 }
1152 goto L90;
1153 }
1154
1155 if (nco == 1) {
1156 i__1 = nro;
1157 for (k = 1; k <= i__1; ++k) {
1158 y += fact[irstk[k + (istk + 1) * irstk_dim1]];
1159 /* L80: */
1160 }
1161 goto L90;
1162 }
1163
1164 lstk[istk] = l;
1165 mstk[istk] = m;
1166 nstk[istk] = n;
1167 ++istk;
1168 nrstk[istk] = nro;
1169 ncstk[istk] = nco;
1170 ystk[istk] = y;
1171 l = 1;
1172 goto L50;
1173
1174 L90:
1175 if (y > amx) {
1176 amx = y;
1177 if (*dsp - amx <= *tol) {
1178 *dsp = 0.f;
1179 goto L9000;
1180 }
1181 }
1182
1183 L100:
1184 --istk;
1185 if (istk == 0) {
1186 *dsp -= amx;
1187 if (*dsp - amx <= *tol) {
1188 *dsp = 0.f;
1189 }
1190 goto L9000;
1191 }
1192 l = lstk[istk] + 1;
1193
1194 L110:
1195 if (l > mstk[istk]) {
1196 goto L100;
1197 }
1198 n = nstk[istk];
1199 nro = nrstk[istk];
1200 nco = ncstk[istk];
1201 y = ystk[istk];
1202 if (n == 1) {
1203 if (irstk[l + istk * irstk_dim1] < irstk[l - 1 + istk * irstk_dim1]) {
1204 goto L60;
1205 }
1206 } else if (n == 2) {
1207 if (icstk[l + istk * icstk_dim1] < icstk[l - 1 + istk * icstk_dim1]) {
1208 goto L60;
1209 }
1210 }
1211
1212 ++l;
1213 goto L110;
1214 L9000:
1215 return 0;
1216 } /* f4xact_ */
1217
1218
1219 /* ----------------------------------------------------------------------- */
1220 /* Name: F3XACT */
1221
1222 /* Purpose: Computes the shortest path length for a given table. */
1223
1224 /* Usage: CALL F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, */
1225 /* IRO, IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, */
1226 /* TOL) */
1227
1228 /* Arguments: */
1229 /* NROW - The number of rows in the table. (Input) */
1230 /* IROW - Vector of length NROW containing the row sums for the */
1231 /* table. (Input) */
1232 /* NCOL - The number of columns in the table. (Input) */
1233 /* ICOL - Vector of length K containing the column sums for the */
1234 /* table. (Input) */
1235 /* DLP - The longest path for the table. (Output) */
1236 /* MM - The total count in the table. (Output) */
1237 /* FACT - Vector containing the logarithms of factorials. (Input) */
1238 /* ICO - Work vector of length MAX(NROW,NCOL). */
1239 /* IRO - Work vector of length MAX(NROW,NCOL). */
1240 /* IT - Work vector of length MAX(NROW,NCOL). */
1241 /* LB - Work vector of length MAX(NROW,NCOL). */
1242 /* NR - Work vector of length MAX(NROW,NCOL). */
1243 /* NT - Work vector of length MAX(NROW,NCOL). */
1244 /* NU - Work vector of length MAX(NROW,NCOL). */
1245 /* ITC - Work vector of length 400. */
1246 /* IST - Work vector of length 400. */
1247 /* STV - Work vector of length 400. */
1248 /* ALEN - Work vector of length MAX(NROW,NCOL). */
1249 /* TOL - Tolerance. (Input) */
1250 /* ----------------------------------------------------------------------- */
f3xact_(long * nrow,long * irow,long * ncol,long * icol,double * dlp,long * mm,double * fact,long * ico,long * iro,long * it,long * lb,long * nr,long * nt,long * nu,long * itc,long * ist,double * stv,double * alen,double * tol)1251 int f3xact_(long *nrow, long *irow, long *ncol,
1252 long *icol, double *dlp, long *mm, double *fact,
1253 long *ico, long *iro, long *it, long *lb, long *nr,
1254 long *nt, long *nu, long *itc, long *ist, double *stv,
1255 double *alen, double *tol)
1256 {
1257 /* Initialized data */
1258
1259 long ldst = 200;
1260 long nst = 0;
1261 long nitc = 0;
1262
1263 /* System generated locals */
1264 long i__1;
1265 double d__1, d__2;
1266
1267 /* Local variables */
1268 long i__, k;
1269 double v;
1270 long n11, n12, ii, nn, ks, ic1, ic2, nc1, nn1, nr1, nco;
1271 double val;
1272 long nct, ipn, irl, key, lev, itp, nro;
1273 double vmn;
1274 long nrt, kyy, nc1s;
1275 bool xmin;
1276 _String errMsg;
1277
1278 /* Parameter adjustments */
1279 --stv;
1280 --ist;
1281 --itc;
1282 --nu;
1283 --nt;
1284 --nr;
1285 --lb;
1286 --it;
1287 --iro;
1288 --ico;
1289 --icol;
1290 --irow;
1291
1292 /* Function Body */
1293
1294 i__1 = *ncol;
1295 for (i__ = 0; i__ <= i__1; ++i__) {
1296 alen[i__] = 0.f;
1297 /* L10: */
1298 }
1299 for (i__ = 1; i__ <= 400; ++i__) {
1300 ist[i__] = -1;
1301 /* L20: */
1302 }
1303 /* nrow is 1 */
1304 if (*nrow <= 1) {
1305 if (*nrow > 0) {
1306 *dlp -= fact[icol[1]];
1307 i__1 = *ncol;
1308 for (i__ = 2; i__ <= i__1; ++i__) {
1309 *dlp -= fact[icol[i__]];
1310 /* L30: */
1311 }
1312 }
1313 goto L9000;
1314 }
1315 /* ncol is 1 */
1316 if (*ncol <= 1) {
1317 if (*ncol > 0) {
1318 *dlp = *dlp - fact[irow[1]] - fact[irow[2]];
1319 i__1 = *nrow;
1320 for (i__ = 3; i__ <= i__1; ++i__) {
1321 *dlp -= fact[irow[i__]];
1322 /* L40: */
1323 }
1324 }
1325 goto L9000;
1326 }
1327 /* 2 by 2 table */
1328 if (*nrow * *ncol == 4) {
1329 n11 = (irow[1] + 1) * (icol[1] + 1) / (*mm + 2);
1330 n12 = irow[1] - n11;
1331 *dlp = *dlp - fact[n11] - fact[n12] - fact[icol[1] - n11] - fact[icol[
1332 2] - n12];
1333 goto L9000;
1334 }
1335 /* Test for optimal table */
1336 val = 0.f;
1337 xmin = false;
1338 if (irow[*nrow] <= irow[1] + *ncol) {
1339 f10act_(*nrow, &irow[1], *ncol, &icol[1], &val, &xmin, fact, &lb[1], &
1340 nu[1], &nr[1]);
1341 }
1342 if (! xmin) {
1343 if (icol[*ncol] <= icol[1] + *nrow) {
1344 f10act_(*ncol, &icol[1], *nrow, &irow[1], &val, &xmin, fact, &lb[1],
1345 &nu[1], &nr[1]);
1346 }
1347 }
1348
1349 if (xmin) {
1350 *dlp -= val;
1351 goto L9000;
1352 }
1353 /* Setup for dynamic programming */
1354 nn = *mm;
1355 /* Minimize ncol */
1356 if (*nrow >= *ncol) {
1357 nro = *nrow;
1358 nco = *ncol;
1359
1360 i__1 = *nrow;
1361 for (i__ = 1; i__ <= i__1; ++i__) {
1362 iro[i__] = irow[i__];
1363 /* L50: */
1364 }
1365
1366 ico[1] = icol[1];
1367 nt[1] = nn - ico[1];
1368 i__1 = *ncol;
1369 for (i__ = 2; i__ <= i__1; ++i__) {
1370 ico[i__] = icol[i__];
1371 nt[i__] = nt[i__ - 1] - ico[i__];
1372 /* L60: */
1373 }
1374 } else {
1375 nro = *ncol;
1376 nco = *nrow;
1377
1378 ico[1] = irow[1];
1379 nt[1] = nn - ico[1];
1380 i__1 = *nrow;
1381 for (i__ = 2; i__ <= i__1; ++i__) {
1382 ico[i__] = irow[i__];
1383 nt[i__] = nt[i__ - 1] - ico[i__];
1384 /* L70: */
1385 }
1386
1387 i__1 = *ncol;
1388 for (i__ = 1; i__ <= i__1; ++i__) {
1389 iro[i__] = icol[i__];
1390 /* L80: */
1391 }
1392 }
1393 /* Initialize pointers */
1394 vmn = 1e10;
1395 nc1s = nco - 1;
1396 irl = 1;
1397 ks = 0;
1398 k = ldst;
1399 kyy = ico[nco] + 1;
1400 goto L100;
1401 /* Test for optimality */
1402 L90:
1403 xmin = false;
1404 if (iro[nro] <= iro[irl] + nco) {
1405 f10act_(nro, &iro[irl], nco, &ico[1], &val, &xmin, fact, &lb[1], &
1406 nu[1], &nr[1]);
1407 }
1408 if (! xmin) {
1409 if (ico[nco] <= ico[1] + nro) {
1410 f10act_(nco, &ico[1], nro, &iro[irl], &val, &xmin, fact, &lb[1],
1411 &nu[1], &nr[1]);
1412 }
1413 }
1414
1415 if (xmin) {
1416 if (val < vmn) {
1417 vmn = val;
1418 }
1419 goto L200;
1420 }
1421 /* Setup to generate new node */
1422 L100:
1423 lev = 1;
1424 nr1 = nro - 1;
1425 nrt = iro[irl];
1426 nct = ico[1];
1427 lb[1] = (long) ((double) ((nrt + 1) * (nct + 1)) / (double) (
1428 nn + nr1 * nc1s + 1) - *tol) - 1;
1429 nu[1] = (long) ((double) ((nrt + nc1s) * (nct + nr1)) / (
1430 double) (nn + nr1 + nc1s)) - lb[1] + 1;
1431 nr[1] = nrt - lb[1];
1432 /* Generate a node */
1433 L110:
1434 --nu[lev];
1435 if (nu[lev] == 0) {
1436 if (lev == 1) {
1437 goto L200;
1438 }
1439 --lev;
1440 goto L110;
1441 }
1442 ++lb[lev];
1443 --nr[lev];
1444 L120:
1445 alen[lev] = alen[lev - 1] + fact[lb[lev]];
1446 if (lev < nc1s) {
1447 nn1 = nt[lev];
1448 nrt = nr[lev];
1449 ++lev;
1450 nc1 = nco - lev;
1451 nct = ico[lev];
1452 lb[lev] = (long) ((double) ((nrt + 1) * (nct + 1)) / (
1453 double) (nn1 + nr1 * nc1 + 1) - *tol);
1454 nu[lev] = (long) ((double) ((nrt + nc1) * (nct + nr1)) / (
1455 double) (nn1 + nr1 + nc1) - lb[lev] + 1);
1456 nr[lev] = nrt - lb[lev];
1457 goto L120;
1458 }
1459 alen[nco] = alen[lev] + fact[nr[lev]];
1460 lb[nco] = nr[lev];
1461
1462 v = val + alen[nco];
1463 if (nro == 2) {
1464 /* Only 1 row left */
1465 v = v + fact[ico[1] - lb[1]] + fact[ico[2] - lb[2]];
1466 i__1 = nco;
1467 for (i__ = 3; i__ <= i__1; ++i__) {
1468 v += fact[ico[i__] - lb[i__]];
1469 /* L130: */
1470 }
1471 if (v < vmn) {
1472 vmn = v;
1473 }
1474 } else if (nro == 3 && nco == 2) {
1475 /* 3 rows and 2 columns */
1476 nn1 = nn - iro[irl] + 2;
1477 ic1 = ico[1] - lb[1];
1478 ic2 = ico[2] - lb[2];
1479 n11 = (iro[irl + 1] + 1) * (ic1 + 1) / nn1;
1480 n12 = iro[irl + 1] - n11;
1481 v = v + fact[n11] + fact[n12] + fact[ic1 - n11] + fact[ic2 - n12];
1482 if (v < vmn) {
1483 vmn = v;
1484 }
1485 } else {
1486 /* Column marginals are new node */
1487 i__1 = nco;
1488 for (i__ = 1; i__ <= i__1; ++i__) {
1489 it[i__] = ico[i__] - lb[i__];
1490 /* L140: */
1491 }
1492 /* Sort column marginals */
1493 if (nco == 2) {
1494 if (it[1] > it[2]) {
1495 ii = it[1];
1496 it[1] = it[2];
1497 it[2] = ii;
1498 }
1499 } else if (nco == 3) {
1500 ii = it[1];
1501 if (ii > it[3]) {
1502 if (ii > it[2]) {
1503 if (it[2] > it[3]) {
1504 it[1] = it[3];
1505 it[3] = ii;
1506 } else {
1507 it[1] = it[2];
1508 it[2] = it[3];
1509 it[3] = ii;
1510 }
1511 } else {
1512 it[1] = it[3];
1513 it[3] = it[2];
1514 it[2] = ii;
1515 }
1516 } else if (ii > it[2]) {
1517 it[1] = it[2];
1518 it[2] = ii;
1519 } else if (it[2] > it[3]) {
1520 ii = it[2];
1521 it[2] = it[3];
1522 it[3] = ii;
1523 }
1524 } else {
1525 isort_(&nco, &it[1]);
1526 }
1527 /* Compute hash value */
1528 key = it[1] * kyy + it[2];
1529 for (i__ = 3; i__ <= nco; ++i__) {
1530 key = it[i__] + key * kyy;
1531 /* L150: */
1532 }
1533 /* Table index */
1534
1535 /*if(key < 0)
1536 {
1537 errMsg = "Negative key in f3xact";
1538 WarnError (errMsg);
1539 //return 0;
1540 }*/
1541
1542 if ((ipn = key % ldst + 1) < 1) {
1543 ipn += ldst;
1544 }
1545 /* Find empty position */
1546 ii = ks + ipn;
1547
1548 for (itp = ipn; itp <= ldst; ++itp) {
1549 if (ist[ii] < 0) {
1550 goto L180;
1551 } else if (ist[ii] == key) {
1552 goto L190;
1553 }
1554 ++ii;
1555 /* L160: */
1556 }
1557
1558 ii = ks + 1;
1559 i__1 = ipn - 1;
1560 for (itp = 1; itp <= i__1; ++itp) {
1561 if (ist[ii] < 0) {
1562 goto L180;
1563 } else if (ist[ii] == key) {
1564 goto L190;
1565 }
1566 ++ii;
1567 /* L170: */
1568 }
1569
1570 errMsg = ("Fisher Exact: Stack length exceeded in f3xact. This problem should not occur.");
1571 ReportWarning (errMsg);
1572 return 0;
1573 /* Push onto stack */
1574 L180:
1575 ist[ii] = key;
1576 stv[ii] = v;
1577 ++nst;
1578 ii = nst + ks;
1579 itc[ii] = itp;
1580 goto L110;
1581 /* Marginals already on stack */
1582 L190:
1583 /* Computing MIN */
1584 d__1 = v, d__2 = stv[ii];
1585 stv[ii] = Minimum(d__1,d__2);
1586 }
1587 goto L110;
1588 /* Pop item from stack */
1589 L200:
1590 if (nitc > 0) {
1591 /* Stack index */
1592 itp = itc[nitc + k] + k;
1593 --nitc;
1594 val = stv[itp];
1595 key = ist[itp];
1596 ist[itp] = -1;
1597 /* Compute marginals */
1598 for (i__ = nco; i__ >= 2; --i__) {
1599 ico[i__] = key % kyy;
1600 key /= kyy;
1601 /* L210: */
1602 }
1603 ico[1] = key;
1604 /* Set up nt array */
1605 nt[1] = nn - ico[1];
1606 i__1 = nco;
1607 for (i__ = 2; i__ <= i__1; ++i__) {
1608 nt[i__] = nt[i__ - 1] - ico[i__];
1609 /* L220: */
1610 }
1611 goto L90;
1612
1613 } else if (nro > 2 && nst > 0) {
1614 /* Go to next level */
1615 nitc = nst;
1616 nst = 0;
1617 k = ks;
1618 ks = ldst - ks;
1619 nn -= iro[irl];
1620 ++irl;
1621 --nro;
1622 goto L200;
1623 }
1624
1625 *dlp -= vmn;
1626 L9000:
1627 return 0;
1628 } /* f3xact_ */
1629
1630 /* ----------------------------------------------------------------------- */
1631 /* Name: F2XACT */
1632
1633 /* Purpose: Computes Fisher's exact test for a contingency table, */
1634 /* routine with workspace variables specified. */
1635
1636 /* Usage: CALL F2XACT (NROW, NCOL, TABLE, EXPECT, PERCNT, */
1637 /* EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF, */
1638 /* IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ, */
1639 /* DLP, DSP, TM, KEY2, IWK, RWK) */
1640 /* ----------------------------------------------------------------------- */
f2xact_(long * nrow,long * ncol,double * table,double * expect,double * percnt,double * emin,double * prt,double * pre,double * fact,long * ico,long * iro,long * kyy,long * idif,long * irn,long * key,long * ldkey,long * ipoin,double * stp,long * ldstp,long * ifrq,double * dlp,double * dsp,double * tm,long * key2,long * iwk,double * rwk)1641 int f2xact_(long *nrow, long *ncol, double *table,
1642 double *expect, double *percnt, double *
1643 emin, double *prt, double *pre, double *fact, long *
1644 ico, long *iro, long *kyy, long *idif, long *irn, long
1645 *key, long *ldkey, long *ipoin, double *stp, long *ldstp,
1646 long *ifrq, double *dlp, double *dsp, double *tm,
1647 long *key2, long *iwk, double *rwk)
1648 {
1649 /* Initialized data */
1650
1651 static long imax = 2147483647;
1652 static double amiss = -12345.f;
1653 static double tol = 3.45254e-7;
1654 static double emx = 1e30f;
1655
1656 /* System generated locals */
1657 long table_dim1, table_offset, i__1, i__2;
1658 double d__1, d__2, d__3, d__4;
1659
1660 /* Local variables */
1661 long i__, j, k, n, k1;
1662 double dd, df;
1663 long i31, i32, i33, i34, i35, i36, i37, i38, i39, i41, i42, i43,
1664 i44, i45, i46, i47, i48, ii, kb, kd, ks;
1665 double pv;
1666 long i310, i311;
1667 double ddf;
1668 long nco, nrb;
1669 double emn, drn, dro, obs;
1670 long ipn, ipo, itp, nro;
1671 double tmp, obs2, obs3;
1672 long nro2, kval, kmax, jkey, last;
1673 bool ipsh;
1674 long itmp;
1675 double dspt;
1676 long itop, jstp, ntot, jstp2, jstp3, jstp4, iflag, ncell, ifreq;
1677 bool chisq = false;
1678 long ikkey;
1679 double pastp;
1680 long ikstp;
1681 long ikstp2;
1682 long ifault;
1683
1684 _String errMsg ("Fisher Exact:");
1685
1686 /* SPECIFICATIONS FOR ARGUMENTS */
1687 /* SPECIFICATIONS FOR LOCAL VARIABLES */
1688 /* SPECIFICATIONS FOR INTRINSICS */
1689 /* SPECIFICATIONS FOR SUBROUTINES */
1690 /* SPECIFICATIONS FOR FUNCTIONS */
1691 /* *********************************************************************** */
1692 /* IMAX is the largest representable */
1693 /* long on the machine */
1694 /* *********************************************************************** */
1695 /* Parameter adjustments */
1696 table_dim1 = *nrow;
1697 table_offset = 1 + table_dim1;
1698 table -= table_offset;
1699 --ico;
1700 --iro;
1701 --kyy;
1702 --idif;
1703 --irn;
1704 --key;
1705 --ipoin;
1706 --stp;
1707 --ifrq;
1708 --dlp;
1709 --dsp;
1710 --tm;
1711 --key2;
1712 --iwk;
1713 --rwk;
1714
1715 /* Function Body */
1716 /* *********************************************************************** */
1717 /* AMISS is a missing value indicator */
1718 /* which is returned when the */
1719 /* probability is not defined. */
1720 /* *********************************************************************** */
1721 /* *********************************************************************** */
1722 /* TOL is chosen as the square root of */
1723 /* the smallest relative spacing */
1724 /* *********************************************************************** */
1725 /* *********************************************************************** */
1726 /* EMX is a large positive value used */
1727 /* in comparing expected values */
1728 /* *********************************************************************** */
1729 /* Initialize KEY array */
1730 i__1 = 2 * *ldkey;
1731 for (i__ = 1; i__ <= i__1; ++i__) {
1732 key[i__] = -9999;
1733 key2[i__] = -9999;
1734 /* L10: */
1735 }
1736 /* Initialize parameters */
1737 *pre = 0.f;
1738 itop = 0;
1739 if (*expect > 0.) {
1740 emn = *emin;
1741 } else {
1742 emn = emx;
1743 }
1744 /* Initialize pointers for workspace */
1745 k = Maximum(*nrow,*ncol);
1746 /* f3xact */
1747 i31 = 1;
1748 i32 = i31 + k;
1749 i33 = i32 + k;
1750 i34 = i33 + k;
1751 i35 = i34 + k;
1752 i36 = i35 + k;
1753 i37 = i36 + k;
1754 i38 = i37 + k;
1755 i39 = i38 + 400;
1756 i310 = 1;
1757 i311 = 401;
1758 /* f4xact */
1759 k = *nrow + *ncol + 1;
1760 i41 = 1;
1761 i42 = i41 + k;
1762 i43 = i42 + k;
1763 i44 = i43 + k;
1764 i45 = i44 + k;
1765 i46 = i45 + k;
1766 i47 = i46 + k * Maximum(*nrow,*ncol);
1767 i48 = 1;
1768 /* Check table dimensions */
1769 /* if (*nrow > *ldtabl) {
1770 errMsg = errMsg & "NROW must be less than or equal to LDTABL.";
1771 WarnError (errMsg);
1772 return 0;
1773 }*/
1774 if (*ncol <= 1) {
1775 HandleApplicationError( errMsg & "NCOL must be greater than 1.0");
1776 return 0;
1777 }
1778 /* Compute row marginals and total */
1779 ntot = 0;
1780 i__1 = *nrow;
1781 for (i__ = 1; i__ <= i__1; ++i__) {
1782 iro[i__] = 0;
1783 i__2 = *ncol;
1784 for (j = 1; j <= i__2; ++j) {
1785 if (table[i__ + j * table_dim1] < -1e-4) {
1786 HandleApplicationError (errMsg & "All elements of TABLE must be positive.");
1787 return 0;
1788 }
1789 iro[i__] += i_dnnt(&table[i__ + j * table_dim1]);
1790 ntot += i_dnnt(&table[i__ + j * table_dim1]);
1791 /* L20: */
1792 }
1793 /* L30: */
1794 }
1795
1796 if (ntot == 0) {
1797 errMsg = errMsg & "All elements of TABLE are zero. PRT and PRE are set to missing values (NaN, not a number)";
1798 ReportWarning (errMsg);
1799 return 0;
1800 *prt = amiss;
1801 *pre = amiss;
1802 goto L9000;
1803 }
1804 /* Column marginals */
1805 i__1 = *ncol;
1806 for (i__ = 1; i__ <= i__1; ++i__) {
1807 ico[i__] = 0;
1808 i__2 = *nrow;
1809 for (j = 1; j <= i__2; ++j) {
1810 ico[i__] += i_dnnt(&table[j + i__ * table_dim1]);
1811 /* L40: */
1812 }
1813 /* L50: */
1814 }
1815 /* sort */
1816 isort_(nrow, &iro[1]);
1817 isort_(ncol, &ico[1]);
1818 /* Determine row and column marginals */
1819
1820 if (*nrow > *ncol) {
1821 nro = *ncol;
1822 nco = *nrow;
1823 /* Interchange row and column marginals */
1824 i__1 = *nrow;
1825 for (i__ = 1; i__ <= i__1; ++i__) {
1826 itmp = iro[i__];
1827 if (i__ <= *ncol) {
1828 iro[i__] = ico[i__];
1829 }
1830 ico[i__] = itmp;
1831 /* L60: */
1832 }
1833 } else {
1834 nro = *nrow;
1835 nco = *ncol;
1836 }
1837
1838 /* Get multiplers for stack */
1839 kyy[1] = 1;
1840 i__1 = nro;
1841 for (i__ = 2; i__ <= i__1; ++i__) {
1842 /* Hash table multipliers */
1843 if (iro[i__ - 1] + 1 <= imax / kyy[i__ - 1]) {
1844 kyy[i__] = kyy[i__ - 1] * (iro[i__ - 1] + 1);
1845 j /= kyy[i__ - 1];
1846 } else {
1847 errMsg = errMsg & "The hash table key cannot be computed because the largest key is larger than the largest representable integer. The algorithm cannot proceed.";
1848 ReportWarning (errMsg);
1849 return 0;
1850 }
1851 /* L70: */
1852 }
1853 /* Maximum product */
1854 if (iro[nro - 1] + 1 <= imax / kyy[nro - 1]) {
1855 kmax = (iro[nro] + 1) * kyy[nro - 1];
1856 } else {
1857 errMsg = errMsg & "The hash table key cannot be computed because the largest key is larger than the largest representable integer. The algorithm cannot proceed.";
1858 ReportWarning (errMsg);
1859 goto L9000;
1860 }
1861 /* Compute log factorials */
1862 fact[0] = 0.;
1863 fact[1] = 0.;
1864 fact[2] = log(2.);
1865 i__1 = ntot;
1866 for (i__ = 3; i__ <= i__1; i__ += 2) {
1867 fact[i__] = fact[i__ - 1] + log((double) i__);
1868 j = i__ + 1;
1869 if (j <= ntot) {
1870 fact[j] = fact[i__] + fact[2] + fact[j / 2] - fact[j / 2 - 1];
1871 }
1872 /* L80: */
1873 }
1874 /* Compute observed path length: OBS */
1875 obs = tol;
1876 ntot = 0;
1877 i__1 = nco;
1878 for (j = 1; j <= i__1; ++j) {
1879 dd = 0.f;
1880 i__2 = nro;
1881 for (i__ = 1; i__ <= i__2; ++i__) {
1882 if (*nrow <= *ncol) {
1883 dd += fact[i_dnnt(&table[i__ + j * table_dim1])];
1884 ntot += i_dnnt(&table[i__ + j * table_dim1]);
1885 } else {
1886 dd += fact[i_dnnt(&table[j + i__ * table_dim1])];
1887 ntot += i_dnnt(&table[j + i__ * table_dim1]);
1888 }
1889 /* L90: */
1890 }
1891 obs = obs + fact[ico[j]] - dd;
1892 /* L100: */
1893 }
1894 /* Denominator of observed table: DRO */
1895 dro = f9xact_(&nro, &ntot, &iro[1], fact);
1896 *prt = exp(obs - dro);
1897 /* Initialize pointers */
1898 k = nco;
1899 last = *ldkey + 1;
1900 jkey = *ldkey + 1;
1901 jstp = *ldstp + 1;
1902 jstp2 = *ldstp * 3 + 1;
1903 jstp3 = (*ldstp << 2) + 1;
1904 jstp4 = *ldstp * 5 + 1;
1905 ikkey = 0;
1906 ikstp = 0;
1907 ikstp2 = *ldstp << 1;
1908 ipo = 1;
1909 ipoin[1] = 1;
1910 stp[1] = 0.f;
1911 ifrq[1] = 1;
1912 ifrq[ikstp2 + 1] = -1;
1913
1914 L110:
1915 kb = nco - k + 1;
1916 ks = 0;
1917 n = ico[kb];
1918 kd = nro + 1;
1919 kmax = nro;
1920 /* IDIF is the difference in going to th */
1921 /* daughter */
1922 i__1 = nro;
1923 for (i__ = 1; i__ <= i__1; ++i__) {
1924 idif[i__] = 0;
1925 /* L120: */
1926 }
1927 /* Generate the first daughter */
1928 L130:
1929 --kd;
1930 /* Computing MIN */
1931 i__1 = n, i__2 = iro[kd];
1932 ntot = Minimum(i__1,i__2);
1933 idif[kd] = ntot;
1934 if (idif[kmax] == 0) {
1935 --kmax;
1936 }
1937 n -= ntot;
1938 if (n > 0 && kd != 1) {
1939 goto L130;
1940 }
1941 if (n != 0) {
1942 goto L310;
1943 }
1944
1945 k1 = k - 1;
1946 n = ico[kb];
1947 ntot = 0;
1948 i__1 = nco;
1949 for (i__ = kb + 1; i__ <= i__1; ++i__) {
1950 ntot += ico[i__];
1951 /* L140: */
1952 }
1953 /* Arc to daughter length=ICO(KB) */
1954 L150:
1955 i__1 = nro;
1956 for (i__ = 1; i__ <= i__1; ++i__) {
1957 irn[i__] = iro[i__] - idif[i__];
1958 /* L160: */
1959 }
1960 /* Sort irn */
1961 if (k1 > 1) {
1962 if (nro == 2) {
1963 if (irn[1] > irn[2]) {
1964 ii = irn[1];
1965 irn[1] = irn[2];
1966 irn[2] = ii;
1967 }
1968 } else if (nro == 3) {
1969 ii = irn[1];
1970 if (ii > irn[3]) {
1971 if (ii > irn[2]) {
1972 if (irn[2] > irn[3]) {
1973 irn[1] = irn[3];
1974 irn[3] = ii;
1975 } else {
1976 irn[1] = irn[2];
1977 irn[2] = irn[3];
1978 irn[3] = ii;
1979 }
1980 } else {
1981 irn[1] = irn[3];
1982 irn[3] = irn[2];
1983 irn[2] = ii;
1984 }
1985 } else if (ii > irn[2]) {
1986 irn[1] = irn[2];
1987 irn[2] = ii;
1988 } else if (irn[2] > irn[3]) {
1989 ii = irn[2];
1990 irn[2] = irn[3];
1991 irn[3] = ii;
1992 }
1993 } else {
1994 i__1 = nro;
1995 for (j = 2; j <= i__1; ++j) {
1996 i__ = j - 1;
1997 ii = irn[j];
1998 L170:
1999 if (ii < irn[i__]) {
2000 irn[i__ + 1] = irn[i__];
2001 --i__;
2002 if (i__ > 0) {
2003 goto L170;
2004 }
2005 }
2006 irn[i__ + 1] = ii;
2007 /* L180: */
2008 }
2009 }
2010 /* Adjust start for zero */
2011 i__1 = nro;
2012 for (i__ = 1; i__ <= i__1; ++i__) {
2013 if (irn[i__] != 0) {
2014 goto L200;
2015 }
2016 /* L190: */
2017 }
2018 L200:
2019 nrb = i__;
2020 nro2 = nro - i__ + 1;
2021 } else {
2022 nrb = 1;
2023 nro2 = nro;
2024 }
2025 /* Some table values */
2026 ddf = f9xact_(&nro, &n, &idif[1], fact);
2027 drn = f9xact_(&nro2, &ntot, &irn[nrb], fact) - dro + ddf;
2028 /* Get hash value */
2029 if (k1 > 1) {
2030 kval = irn[1] + irn[2] * kyy[2];
2031 i__1 = nro;
2032 for (i__ = 3; i__ <= i__1; ++i__) {
2033 kval += irn[i__] * kyy[i__];
2034 /* L210: */
2035 }
2036 /* Get hash table entry */
2037 i__ = kval % (*ldkey << 1) + 1;
2038 /* Search for unused location */
2039 i__1 = *ldkey << 1;
2040 for (itp = i__; itp <= i__1; ++itp) {
2041 ii = key2[itp];
2042 if (ii == kval) {
2043 goto L240;
2044 } else if (ii < 0) {
2045 key2[itp] = kval;
2046 dlp[itp] = 1.;
2047 dsp[itp] = 1.;
2048 goto L240;
2049 }
2050 /* L220: */
2051 }
2052
2053 i__1 = i__ - 1;
2054 for (itp = 1; itp <= i__1; ++itp) {
2055 ii = key2[itp];
2056 if (ii == kval) {
2057 goto L240;
2058 } else if (ii < 0) {
2059 key2[itp] = kval;
2060 dlp[itp] = 1.f;
2061 goto L240;
2062 }
2063 /* L230: */
2064 }
2065
2066 errMsg = errMsg & "LDKEY is too small. It is not possible to give thevalue of LDKEY required, but you could try doubling LDKEY (and possibly LDSTP).";
2067 ReportWarning (errMsg);
2068 return 0;
2069 }
2070
2071 L240:
2072 ipsh = true;
2073 /* Recover pastp */
2074 ipn = ipoin[ipo + ikkey];
2075 pastp = stp[ipn + ikstp];
2076 ifreq = ifrq[ipn + ikstp];
2077 /* Compute shortest and longest path */
2078 if (k1 > 1) {
2079 obs2 = obs - fact[ico[kb + 1]] - fact[ico[kb + 2]] - ddf;
2080 i__1 = k1;
2081 for (i__ = 3; i__ <= i__1; ++i__) {
2082 obs2 -= fact[ico[kb + i__]];
2083 /* L250: */
2084 }
2085
2086 if (dlp[itp] > 0.) {
2087 dspt = obs - obs2 - ddf;
2088 /* Compute longest path */
2089 dlp[itp] = 0.;
2090 f3xact_(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dlp[itp], &ntot,
2091 fact, &iwk[i31], &iwk[i32], &iwk[i33], &iwk[i34], &iwk[
2092 i35], &iwk[i36], &iwk[i37], &iwk[i38], &iwk[i39], &rwk[
2093 i310], &rwk[i311], &tol);
2094 /* Computing Minimum */
2095 d__1 = 0., d__2 = dlp[itp];
2096 dlp[itp] = Minimum(d__1,d__2);
2097 /* Compute shortest path */
2098 dsp[itp] = dspt;
2099 f4xact_(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dsp[itp], fact, &
2100 iwk[i47], &iwk[i41], &iwk[i42], &iwk[i43], &iwk[i44], &
2101 iwk[i45], &iwk[i46], &rwk[i48], &tol);
2102 /* Computing Minimum */
2103 d__1 = 0., d__2 = dsp[itp] - dspt;
2104 dsp[itp] = Minimum(d__1,d__2);
2105 /* Use chi-squared approximation? */
2106 if ((double) (irn[nrb] * ico[kb + 1]) / (double) ntot >
2107 emn) {
2108 ncell = 0.f;
2109 i__1 = nro2;
2110 for (i__ = 1; i__ <= i__1; ++i__) {
2111 i__2 = k1;
2112 for (j = 1; j <= i__2; ++j) {
2113 if ((double) (irn[nrb + i__ - 1] * ico[kb + j]) >=
2114 ntot * *expect) {
2115 ++ncell;
2116 }
2117 /* L260: */
2118 }
2119 /* L270: */
2120 }
2121 if ((double) (ncell * 100) >= k1 * nro2 * *percnt) {
2122 tmp = 0.f;
2123 i__1 = nro2;
2124 for (i__ = 1; i__ <= i__1; ++i__) {
2125 tmp = tmp + fact[irn[nrb + i__ - 1]] - fact[irn[nrb +
2126 i__ - 1] - 1];
2127 /* L280: */
2128 }
2129 tmp *= k1 - 1;
2130 i__1 = k1;
2131 for (j = 1; j <= i__1; ++j) {
2132 tmp += (nro2 - 1) * (fact[ico[kb + j]] - fact[ico[kb
2133 + j] - 1]);
2134 /* L290: */
2135 }
2136 df = (double) ((nro2 - 1) * (k1 - 1));
2137 tmp += df * 1.83787706640934548356065947281;
2138 tmp -= (nro2 * k1 - 1) * (fact[ntot] - fact[ntot - 1]);
2139 tm[itp] = (obs - dro) * -2. - tmp;
2140 } else {
2141 /* tm(itp) set to a flag value */
2142 tm[itp] = -9876.;
2143 }
2144 } else {
2145 tm[itp] = -9876.;
2146 }
2147 }
2148 obs3 = obs2 - dlp[itp];
2149 obs2 -= dsp[itp];
2150 if (tm[itp] == -9876.) {
2151 chisq = false;
2152 } else {
2153 chisq = true;
2154 tmp = tm[itp];
2155 }
2156 } else {
2157 obs2 = obs - drn - dro;
2158 obs3 = obs2;
2159 }
2160 /* Process node with new PASTP */
2161 L300:
2162 if (pastp <= obs3) {
2163 /* Update pre */
2164 *pre += (double) ifreq * exp(pastp + drn);
2165
2166 } else if (pastp < obs2) {
2167 if (chisq) {
2168 df = (double) ((nro2 - 1) * (k1 - 1));
2169 /* Computing MAX */
2170 d__2 = 0., d__3 = tmp + (pastp + drn) * 2.;
2171 d__1 = Maximum(d__2,d__3) / 2.;
2172 d__4 = df / 2.;
2173 pv = 1.f - gammds_(&d__1, &d__4, &ifault);
2174 *pre += (double) ifreq * exp(pastp + drn) * pv;
2175 } else {
2176 /* Put daughter on queue */
2177 d__1 = pastp + ddf;
2178 f5xact_(&d__1, &tol, &kval, &key[jkey], ldkey, &ipoin[jkey], &stp[
2179 jstp], ldstp, &ifrq[jstp], &ifrq[jstp2], &ifrq[jstp3], &
2180 ifrq[jstp4], &ifreq, &itop, &ipsh);
2181 ipsh = false;
2182 }
2183 }
2184 /* Get next PASTP on chain */
2185 ipn = ifrq[ipn + ikstp2];
2186 if (ipn > 0) {
2187 pastp = stp[ipn + ikstp];
2188 ifreq = ifrq[ipn + ikstp];
2189 goto L300;
2190 }
2191 /* Generate a new daughter node */
2192 f7xact_(&kmax, &iro[1], &idif[1], &kd, &ks, &iflag);
2193 if (iflag != 1) {
2194 goto L150;
2195 }
2196 /* Go get a new mother from stage K */
2197 L310:
2198 iflag = 1;
2199 f6xact_(&nro, &iro[1], &iflag, &kyy[1], &key[ikkey + 1], ldkey, &last, &
2200 ipo);
2201 /* Update pointers */
2202 if (iflag == 3) {
2203 --k;
2204 itop = 0;
2205 ikkey = jkey - 1;
2206 ikstp = jstp - 1;
2207 ikstp2 = jstp2 - 1;
2208 jkey = *ldkey - jkey + 2;
2209 jstp = *ldstp - jstp + 2;
2210 jstp2 = (*ldstp << 1) + jstp;
2211 i__1 = *ldkey << 1;
2212 for (i__ = 1; i__ <= i__1; ++i__) {
2213 key2[i__] = -9999;
2214 /* L320: */
2215 }
2216 if (k >= 2) {
2217 goto L310;
2218 }
2219 } else {
2220 goto L110;
2221 }
2222
2223 L9000:
2224 return 0;
2225 } /* f2xact_ */
2226
2227 /* ----------------------------------------------------------------------- */
2228 /* Name: FEXACT */
2229
2230 /* Purpose: Computes Fisher's exact test probabilities and a hybrid */
2231 /* approximation to Fisher exact test probabilities for a */
2232 /* contingency table using the network algorithm. */
2233
2234 /* Usage: CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT, */
2235 /* EMIN, PRT, PRE) */
2236
2237 /* Arguments: */
2238 /* NROW - The number of rows in the table. (Input) */
2239 /* NCOL - The number of columns in the table. (Input) */
2240 /* TABLE - NROW by NCOL matrix containing the contingency table. */
2241 /* (Input) */
2242 /* LDTABL - Leading dimension of TABLE exactly as specified in the */
2243 /* dimension statement in the calling program. (Input) */
2244 /* EXPECT - Expected value used in the hybrid algorithm for */
2245 /* deciding when to use asymptotic theory probabilities. */
2246 /* (Input) */
2247 /* If EXPECT .LE. 0.0 then asymptotic theory probabilities */
2248 /* are not used and Fisher exact test probabilities are */
2249 /* computed. Otherwise, if PERCNT or more of the cells in */
2250 /* the remaining table have estimated expected values of */
2251 /* EXPECT or more, with no remaining cell having expected */
2252 /* value less than EMIN, then asymptotic chi-squared */
2253 /* probabilities are used. See the algorithm section of the */
2254 /* manual document for details. Use EXPECT = 5.0 to obtain */
2255 /* the 'Cochran' condition. */
2256 /* PERCNT - Percentage of remaining cells that must have estimated */
2257 /* expected values greater than EXPECT before asymptotic */
2258 /* probabilities can be used. (Input) */
2259 /* See argument EXPECT for details. Use PERCNT = 80.0 to */
2260 /* obtain the 'Cochran' condition. */
2261 /* EMIN - Minimum cell estimated expected value allowed for */
2262 /* asymptotic chi-squared probabilities to be used. (Input) */
2263 /* See argument EXPECT for details. Use EMIN = 1.0 to */
2264 /* obtain the 'Cochran' condition. */
2265 /* PRT - Probability of the observed table for fixed marginal */
2266 /* totals. (Output) */
2267 /* PRE - Table p-value. (Output) */
2268 /* PRE is the probability of a more extreme table, where */
2269 /* 'extreme' is in a probabilistic sense. */
2270 /* If EXPECT .LT. 0 then the Fisher exact probability */
2271 /* is returned. Otherwise, an approximation to the */
2272 /* Fisher exact probability is computed based upon */
2273 /* asymptotic chi-squared probabilities for ``large'' */
2274 /* table expected values. The user defines ``large'' */
2275 /* through the arguments EXPECT, PERCNT, and EMIN. */
2276
2277 /* Remarks: */
2278 /* 1. For many problems one megabyte or more of workspace can be */
2279 /* required. If the environment supports it, the user should begin */
2280 /* by increasing the workspace used to 200,000 units. */
2281
2282 /* 2. In FEXACT, LDSTP = 30*LDKEY. The proportion of table space used */
2283 /* by STP may be changed by changing the line MULT = 30 below to */
2284 /* another value. */
2285
2286 /* 3. FEXACT may be converted to single precision by setting IREAL = 3, */
2287 /* and converting all DOUBLE PRECISION specifications (except the */
2288 /* specifications for RWRK, IWRK, and DWRK) to REAL. This will */
2289 /* require changing the names and specifications of the intrinsic */
2290 /* functions ALOG, AMAX1, AMIN1, EXP, and REAL. In addition, the */
2291 /* machine specific constants will need to be changed, and the name */
2292 /* DWRK will need to be changed to RWRK in the call to F2XACT. */
2293
2294 /* 4. Machine specific constants are specified and documented in F2XACT. */
2295 /* A missing value code is specified in both FEXACT and F2XACT. */
2296
2297 /* 5. Although not a restriction, is is not generally practical to call */
2298 /* this routine with large tables which are not sparse and in */
2299 /* which the 'hybrid' algorithm has little effect. For example, */
2300 /* although it is feasible to compute exact probabilities for the */
2301 /* table */
2302 /* 1 8 5 4 4 2 2 */
2303 /* 5 3 3 4 3 1 0 */
2304 /* 10 1 4 0 0 0 0, */
2305 /* computing exact probabilities for a similar table which has been */
2306 /* enlarged by the addition of an extra row (or column) may not be */
2307 /* feasible. */
2308 /* ----------------------------------------------------------------------- */
fexact_(long nrow,long ncol,double * table,double expect,double percnt,double emin,double * prt,double * pre)2309 int fexact_(long nrow, long ncol, double *table, double expect, double percnt, double emin, double *prt, double *pre)
2310 {
2311 long ntot;
2312
2313 _String errMsg ("Fisher Exact:");
2314
2315 /* *********************************************************************** */
2316 /* Set IREAL = 4 for DOUBLE PRECISION */
2317 /* Set IREAL = 3 for SINGLE PRECISION */
2318 /* *********************************************************************** */
2319
2320 /* if (*nrow > *ldtabl) {
2321 errMsg = errMsg & "NROW must be less than or equal to LDTABL.";
2322 WarnError(errMsg);
2323 free (equiv_1);
2324 return 0;
2325 }*/
2326
2327 ntot = 0;
2328 for (long i = 0; i<ncol*nrow; i++) {
2329 if (table[i] < 0.) {
2330 HandleApplicationError (errMsg & "All elements of TABLE must be non-negative.");
2331 return 0;
2332 }
2333 ntot += (long)(table[i]+0.5);
2334 }
2335
2336
2337 if (ntot == 0) {
2338 ReportWarning(errMsg & "All elements of TABLE are zero. PRT and PRE are set to missing values (NaN, not a number).");
2339 *prt = -1.;
2340 *pre = -1.;
2341 return 0;
2342 }
2343
2344 long k = nrow + ncol + 1,
2345 kk = k * ncol;
2346
2347
2348 double * i1 = (double*)MemAllocate((ntot + 1)*sizeof(double)),
2349 * irwk = (double*)MemAllocate(Maximum(ncol + 401,k)*sizeof(double));
2350
2351 bool didAlloc = false;
2352
2353 long *i2 = (long*) MemAllocate (ncol*sizeof (long)),
2354 *i3 = (long*) MemAllocate (ncol*sizeof (long)),
2355 *i3a = (long*) MemAllocate (ncol*sizeof (long)),
2356 *i3b = (long*) MemAllocate (nrow*sizeof (long)),
2357 *i3c = (long*) MemAllocate (nrow*sizeof (long)),
2358 *iiwk = (long*) MemAllocate (Maximum(k * 5 + (kk << 1),ncol * 7 + 800)*sizeof (long));
2359
2360
2361 if (!fexact_i4) {
2362 didAlloc = true;
2363 allocate_fexact_keys(4096,30);
2364 }
2365
2366
2367 f2xact_ (&nrow, &ncol, table, &expect, &percnt, &emin, prt, pre, i1, i2, i3, i3a,
2368 i3b, i3c, fexact_i4, &fexact_ldkey, fexact_i5, fexact_i6, &fexact_ldstp, fexact_i7, fexact_i8,
2369 fexact_i9, fexact_i9a, fexact_i10, iiwk, irwk);
2370
2371 free (i1);
2372 free (i2);
2373 free (i3);
2374 free (i3a);
2375 free (i3b);
2376 free (i3c);
2377
2378 free (irwk);
2379 free (iiwk);
2380
2381 if (didAlloc) {
2382 free_fexact_keys();
2383 }
2384 return 0;
2385 } /* fexact_ */
2386
2387
2388
2389