--- /dev/null
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+//#include "ctest.h"
+
+
+#include <limits.h>
+#define SINT_MAX INT_MAX
+
+#define max(a, b) ((a) < (b) ? (b) : (a))
+#define min(a, b) ((a) > (b) ? (b) : (a))
+
+static void f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
+ double *expect, double *percnt, double *emin, double
+ *prt, double *pre, double *fact, int *ico, int
+ *iro, int *kyy, int *idif, int *irn, int *key,
+ int *ldkey, int *ipoin, double *stp, int *ldstp,
+ int *ifrq, double *dlp, double *dsp, double *tm,
+ int *key2, int *iwk, double *rwk);
+static void f3xact(int *nrow, int *irow, int *ncol, int *icol,
+ double *dlp, int *mm, double *fact, int *ico, int
+ *iro, int *it, int *lb, int *nr, int *nt, int
+ *nu, int *itc, int *ist, double *stv, double *alen,
+ const double *tol);
+static void f4xact(int *nrow, int *irow, int *ncol, int *icol,
+ double *dsp, double *fact, int *icstk, int *ncstk,
+ int *lstk, int *mstk, int *nstk, int *nrstk, int
+ *irstk, double *ystk, const double *tol);
+static void f5xact(double *pastp, const double *tol, int *kval, int *key,
+ int *ldkey, int *ipoin, double *stp, int *ldstp,
+ int *ifrq, int *npoin, int *nr, int *nl, int
+ *ifreq, int *itop, int *ipsh);
+static void f6xact(int *nrow, int *irow, int *iflag, int *kyy,
+ int *key, int *ldkey, int *last, int *ipn);
+static void f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
+ int *iflag);
+static void f8xact(int *irow, int *is, int *i1, int *izero, int *new);
+static double f9xact(int *n, int *mm, int *ir, double *fact);
+static void f10act(int *nrow, int *irow, int *ncol, int *icol,
+ double *val, int *xmin, double *fact, int *nd,
+ int *ne, int *m);
+static void f11act(int *irow, int *i1, int *i2, int *new);
+static void prterr(int icode, char *mes);
+static int iwork(int iwkmax, int *iwkpt, int number, int itype);
+// void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
+// double *expect, double *percnt, double *emin, double *prt,
+// double *pre, /* new in C : */ int *workspace);
+ static void isort(int *n, int *ix);
+ static double gammds(double *y, double *p, int *ifault);
+ static double alogam(double *x, int *ifault);
+
+
+
+
+/* The only public function : */
+void
+fexact(int *nrow, int *ncol, double *table, int *ldtabl,
+ double *expect, double *percnt, double *emin, double *prt,
+ double *pre, /* new in C : */ int *workspace)
+{
+
+/*
+ ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
+ THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
+ VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
+ -----------------------------------------------------------------------
+ Name: FEXACT
+ Purpose: Computes Fisher's exact test probabilities and a hybrid
+ approximation to Fisher exact test probabilities for a
+ contingency table using the network algorithm.
+ Usage: CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
+ EMIN, PRT, PRE)
+ Arguments:
+ NROW - The number of rows in the table. (Input)
+ NCOL - The number of columns in the table. (Input)
+ TABLE - NROW by NCOL matrix containing the contingency
+ table. (Input)
+ LDTABL - Leading dimension of TABLE exactly as specified
+ in the dimension statement in the calling
+ program. (Input)
+ EXPECT - Expected value used in the hybrid algorithm for
+ deciding when to use asymptotic theory
+ probabilities. (Input)
+ If EXPECT <= 0.0 then asymptotic theory probabilities
+ are not used and Fisher exact test probabilities are
+ computed. Otherwise, if PERCNT or more of the cells in
+ the remaining table have estimated expected values of
+ EXPECT or more, with no remaining cell having expected
+ value less than EMIN, then asymptotic chi-squared
+ probabilities are used. See the algorithm section of the
+ manual document for details.
+ Use EXPECT = 5.0 to obtain the 'Cochran' condition.
+ PERCNT - Percentage of remaining cells that must have
+ estimated expected values greater than EXPECT
+ before asymptotic probabilities can be used. (Input)
+ See argument EXPECT for details.
+ Use PERCNT = 80.0 to obtain the 'Cochran' condition.
+ EMIN - Minimum cell estimated expected value allowed for
+ asymptotic chi-squared probabilities to be used. (Input)
+ See argument EXPECT for details.
+ Use EMIN = 1.0 to obtain the 'Cochran' condition.
+ PRT - Probability of the observed table for fixed
+ marginal totals. (Output)
+ PRE - Table p-value. (Output)
+ PRE is the probability of a more extreme table,
+ where `extreme' is in a probabilistic sense.
+ If EXPECT < 0 then the Fisher exact probability
+ is returned. Otherwise, an approximation to the
+ Fisher exact probability is computed based upon
+ asymptotic chi-squared probabilities for ``large''
+ table expected values. The user defines ``large''
+ through the arguments EXPECT, PERCNT, and EMIN.
+
+ Remarks:
+ 1. For many problems one megabyte or more of workspace can be
+ required. If the environment supports it, the user should begin
+ by increasing the workspace used to 200,000 units.
+ 2. In FEXACT, LDSTP = 30*LDKEY. The proportion of table space used
+ by STP may be changed by changing the line MULT = 30 below to
+ another value.
+ 3. FEXACT may be converted to single precision by setting IREAL = 3,
+ and converting all DOUBLE PRECISION specifications (except the
+ specifications for RWRK, IWRK, and DWRK) to REAL. This will
+ require changing the names and specifications of the intrinsic
+ functions ALOG, AMAX1, AMIN1, EXP, and REAL. In addition, the
+ machine specific constants will need to be changed, and the name
+ DWRK will need to be changed to RWRK in the call to F2XACT.
+ 4. Machine specific constants are specified and documented in F2XACT.
+ A missing value code is specified in both FEXACT and F2XACT.
+ 5. Although not a restriction, is is not generally practical to call
+ this routine with large tables which are not sparse and in
+ which the 'hybrid' algorithm has little effect. For example,
+ although it is feasible to compute exact probabilities for the
+ table
+ 1 8 5 4 4 2 2
+ 5 3 3 4 3 1 0
+ 10 1 4 0 0 0 0,
+ computing exact probabilities for a similar table which has been
+ enlarged by the addition of an extra row (or column) may not be
+ feasible.
+ -----------------------------------------------------------------------
+ */
+
+ /* CONSTANT Parameters : */
+
+ /* To increase the length of the table of paste path lengths relative
+ to the length of the hash table, increase MULT.
+ */
+ const int mult = 30;
+ /* AMISS is a missing value indicator which is returned when the
+ probability is not defined.
+ */
+ const double amiss = -12345.;
+ /*
+ Set IREAL = 4 for DOUBLE PRECISION
+ Set IREAL = 3 for SINGLE PRECISION
+ */
+#define i_real 4
+#define i_int 2
+
+ /* System generated locals */
+ int ikh;
+ /* Local variables */
+ int nco, nro, ntot, numb, iiwk, irwk;
+ int i, j, k, kk, ldkey, ldstp, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
+ int i3a, i3b, i3c, i9a, iwkmax, iwkpt;
+
+ /* Workspace Allocation (freed at end) */
+ double *equiv;
+ iwkmax = 2 * (int) (*workspace / 2);
+// equiv = (double *) R_alloc(iwkmax / 2, sizeof(double));
+ equiv = (double *) calloc(iwkmax / 2, sizeof(double));
+
+ /* The check could never happen with Calloc!
+ equiv = Calloc(iwkmax / 2, double);
+ if (!equiv) {
+ prterr(0, "Can not allocate specified workspace");
+ } */
+
+#define dwrk (equiv)
+#define iwrk ((int *)equiv)
+#define rwrk ((float *)equiv)
+
+ /* Parameter adjustments */
+ table -= *ldtabl + 1;
+
+ /* Function Body */
+ iwkpt = 0;
+
+ if (*nrow > *ldtabl)
+ prterr(1, "NROW must be less than or equal to LDTABL.");
+
+ ntot = 0;
+ for (i = 1; i <= *nrow; ++i) {
+ for (j = 1; j <= *ncol; ++j) {
+ if (table[i + j * *ldtabl] < 0.)
+ prterr(2, "All elements of TABLE must be positive.");
+ ntot = (int) (ntot + table[i + j * *ldtabl]);
+ }
+ }
+ if (ntot == 0) {
+ prterr(3, "All elements of TABLE are zero.\n"
+ "PRT and PRE are set to missing values.");
+ *prt = amiss;
+ *pre = amiss;
+ goto L_End;
+ }
+
+ nco = max(*nrow, *ncol);
+ nro = *nrow + *ncol - nco;/* = min(*nrow, *ncol) */
+ k = *nrow + *ncol + 1;
+ kk = k * nco;
+
+ ikh = ntot + 1;
+ i1 = iwork(iwkmax, &iwkpt, ikh, i_real);
+ i2 = iwork(iwkmax, &iwkpt, nco, i_int);
+ i3 = iwork(iwkmax, &iwkpt, nco, i_int);
+ i3a = iwork(iwkmax, &iwkpt, nco, i_int);
+ i3b = iwork(iwkmax, &iwkpt, nro, i_int);
+ i3c = iwork(iwkmax, &iwkpt, nro, i_int);
+ ikh = max(k * 5 + (kk << 1), nco * 7 + 800);
+ iiwk= iwork(iwkmax, &iwkpt, ikh, i_int);
+ ikh = max(nco + 401, k);
+ irwk= iwork(iwkmax, &iwkpt, ikh, i_real);
+
+ /* NOTE:
+ What follows below splits the remaining amount iwkmax - iwkpt of
+ (int) workspace into hash tables as follows.
+ type size index
+ INT 2 * ldkey i4 i5 i11
+ REAL 2 * ldkey i8 i9 i10
+ REAL 2 * ldstp i6
+ INT 6 * ldstp i7
+ Hence, we need ldkey times
+ 3 * 2 + 3 * 2 * s + 2 * mult * s + 6 * mult
+ chunks of integer memory, where s = sizeof(REAL) / sizeof(INT).
+ If doubles are used and are twice as long as ints, this gives
+ 18 + 10 * mult
+ so that the value of ldkey can be obtained by dividing available
+ (int) workspace by this number.
+
+ In fact, because iwork() can actually s * n + s - 1 int chunks
+ when allocating a REAL, we use ldkey = available / numb - 1.
+
+ FIXME:
+ Can we always assume that sizeof(double) / sizeof(int) is 2?
+ */
+
+ if (i_real == 4) { /* Double precision reals */
+ numb = 18 + 10 * mult;
+ } else { /* Single precision reals */
+ numb = (mult << 3) + 12;
+ }
+ ldkey = (iwkmax - iwkpt) / numb - 1;
+ ldstp = mult * ldkey;
+ ikh = ldkey << 1; i4 = iwork(iwkmax, &iwkpt, ikh, i_int);
+ ikh = ldkey << 1; i5 = iwork(iwkmax, &iwkpt, ikh, i_int);
+ ikh = ldstp << 1; i6 = iwork(iwkmax, &iwkpt, ikh, i_real);
+ ikh = ldstp * 6; i7 = iwork(iwkmax, &iwkpt, ikh, i_int);
+ ikh = ldkey << 1; i8 = iwork(iwkmax, &iwkpt, ikh, i_real);
+ ikh = ldkey << 1; i9 = iwork(iwkmax, &iwkpt, ikh, i_real);
+ ikh = ldkey << 1; i9a = iwork(iwkmax, &iwkpt, ikh, i_real);
+ ikh = ldkey << 1; i10 = iwork(iwkmax, &iwkpt, ikh, i_int);
+
+ /* To convert to double precision, change RWRK to DWRK in the next CALL.
+ */
+ f2xact(nrow,
+ ncol,
+ &table[*ldtabl + 1],
+ ldtabl,
+ expect,
+ percnt,
+ emin,
+ prt,
+ pre,
+ dwrk + i1,
+ iwrk + i2,
+ iwrk + i3,
+ iwrk + i3a,
+ iwrk + i3b,
+ iwrk + i3c,
+ iwrk + i4,
+ &ldkey,
+ iwrk + i5,
+ dwrk + i6,
+ &ldstp,
+ iwrk + i7,
+ dwrk + i8,
+ dwrk + i9,
+ dwrk + i9a,
+ iwrk + i10,
+ iwrk + iiwk,
+ dwrk + irwk);
+
+L_End:
+ /* Free(equiv); */
+ free(equiv);
+ return;
+}
+
+#undef rwrk
+#undef iwrk
+#undef dwrk
+
+
+/*
+ -----------------------------------------------------------------------
+ Name: F2XACT
+ Purpose: Computes Fisher's exact test for a contingency table,
+ routine with workspace variables specified.
+ Usage: F2XACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
+ EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF,
+ IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ,
+ DLP, DSP, TM, KEY2, IWK, RWK)
+ -----------------------------------------------------------------------
+ */
+void
+f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
+ double *expect, double *percnt, double *emin, double *prt,
+ double *pre, double *fact, int *ico, int *iro, int *kyy,
+ int *idif, int *irn, int *key, int *ldkey, int *ipoin,
+ double *stp, int *ldstp, int *ifrq, double *dlp, double *dsp,
+ double *tm, int *key2, int *iwk, double *rwk)
+{
+ /* IMAX is the largest representable int on the machine. */
+ const int imax = SINT_MAX;
+// const int imax = 2147483647; //xx: I DON´T like this, and
+// thanks to the hint from Jason Turner I don't do it anymore. (R.D-U).
+
+ /* AMISS is a missing value indicator which is returned when the
+ probability is not defined. */
+ const double amiss = -12345.;
+
+ /* TOL is chosen as the square root of the smallest relative spacing. */
+#ifndef Macintosh
+ const static double tol = 3.45254e-7;
+#else
+ static double tol = 3.45254e-7;
+#endif
+ /* EMX is a large positive value used in comparing expected values. */
+ const static double emx = 1e30;
+
+ /* Local variables {{any really need to be static ???}} */
+ static int kval, kmax, jkey, last, ipsh, itmp, itop, jstp, ntot,
+ jstp2, jstp3, jstp4, i, ii, j, k, n, iflag, ncell, ifreq, chisq,
+ ikkey, ikstp, ikstp2, k1, kb, kd, ks,
+ i31, i32, i33, i34, i35, i36, i37, i38, i39,
+ i41, i42, i43, i44, i45, i46, i47, i48, i310, i311,
+ nco, nrb, ipn, ipo, itp, nro, nro2;
+ static double dspt, dd, df,ddf, drn,dro, emn, obs, obs2, obs3,
+ pastp, pv, tmp;
+ double d1;
+#ifndef USING_R
+ double d2;
+ static int ifault;
+#endif
+ int nr_gt_nc=0;
+
+ /* Parameter adjustments */
+ table -= *ldtabl + 1;
+ --ico;
+ --iro;
+ --kyy;
+ --idif;
+ --irn;
+ --key;
+ --ipoin;
+ --stp;
+ --ifrq;
+ --dlp;
+ --dsp;
+ --tm;
+ --key2;
+ --iwk;
+ --rwk;
+
+
+ /* Check table dimensions */
+ if (*nrow > *ldtabl)
+ prterr(1, "NROW must be less than or equal to LDTABL.");
+ if (*ncol <= 1)
+ prterr(4, "NCOL must be at least 2");
+
+ /* Initialize KEY array */
+ for (i = 1; i <= *ldkey << 1; ++i) {
+ key[i] = -9999;
+ key2[i] = -9999;
+ }
+ /* Initialize parameters */
+ *pre = 0.;
+ itop = 0;
+ if (*expect > 0.)
+ emn = *emin;
+ else
+ emn = emx;
+ if (*nrow > *ncol){
+ nr_gt_nc = 1;
+}
+else{
+ nr_gt_nc = 0;
+}
+ /* nco := max(nrow, ncol) : */
+ if(nr_gt_nc)
+ nco = *nrow;
+ else
+ nco = *ncol;
+ /* Initialize pointers for workspace */
+ /* f3xact */
+ i31 = 1;
+ i32 = i31 + nco;
+ i33 = i32 + nco;
+ i34 = i33 + nco;
+ i35 = i34 + nco;
+ i36 = i35 + nco;
+ i37 = i36 + nco;
+ i38 = i37 + nco;
+ i39 = i38 + 400;
+ i310 = 1;
+ i311 = 401;
+ /* f4xact */
+ k = *nrow + *ncol + 1;
+ i41 = 1;
+ i42 = i41 + k;
+ i43 = i42 + k;
+ i44 = i43 + k;
+ i45 = i44 + k;
+ i46 = i45 + k;
+ i47 = i46 + k * nco;
+ i48 = 1;
+
+ /* Compute row marginals and total */
+ ntot = 0;
+ for (i = 1; i <= *nrow; ++i) {
+ iro[i] = 0;
+ for (j = 1; j <= *ncol; ++j) {
+ if (table[i + j * *ldtabl] < -1e-4)
+ prterr(2, "All elements of TABLE must be positive.");
+ iro[i] += (int) table[i + j * *ldtabl];
+ }
+ ntot += iro[i];
+ }
+
+ if (ntot == 0) {
+ prterr(3, "All elements of TABLE are zero.\n"
+ "PRT and PRE are set to missing values.");
+ *pre = *prt = amiss;
+ return;
+ }
+
+ /* Column marginals */
+ for (i = 1; i <= *ncol; ++i) {
+ ico[i] = 0;
+ for (j = 1; j <= *nrow; ++j)
+ ico[i] += (int) table[j + i * *ldtabl];
+ }
+
+ /* sort marginals */
+ isort(nrow, &iro[1]);
+ isort(ncol, &ico[1]);
+
+ /* Determine row and column marginals.
+ Define max(nrow,ncol) =: nco >= nro := min(nrow,ncol)
+ nco is defined above
+
+ Swap marginals if necessary to ico[1:nco] & iro[1:nro]
+ */
+ if (nr_gt_nc) {
+ nro = *ncol;
+ /* Swap marginals */
+ for (i = 1; i <= nco; ++i) {
+ itmp = iro[i];
+ if (i <= nro)
+ iro[i] = ico[i];
+ ico[i] = itmp;
+ }
+ } else
+ nro = *nrow;
+
+
+ /* Get multiplers for stack */
+ kyy[1] = 1;
+ for (i = 2; i <= nro; ++i) {
+ /* Hash table multipliers */
+ if (iro[i - 1] + 1 <= imax / kyy[i - 1]) {
+ kyy[i] = kyy[i - 1] * (iro[i - 1] + 1);
+ j /= kyy[i - 1];
+ } else
+ goto L_ERR_5;
+ }
+ /* Maximum product */
+ if (iro[nro - 1] + 1 <= imax / kyy[nro - 1]) {
+ kmax = (iro[nro] + 1) * kyy[nro - 1];
+ } else {
+ L_ERR_5:
+ prterr(5, "The hash table key cannot be computed because "
+ "the largest key\n"
+ "is larger than the largest representable int.\n"
+ "The algorithm cannot proceed.\n"
+ "Reduce the workspace size, or use `exact = FALSE'.");
+ return;
+ }
+
+ /* Compute log factorials */
+ fact[0] = 0.;
+ fact[1] = 0.;
+ if(ntot >= 2) fact[2] = log(2.);
+/* MM: old code assuming log() to be SLOW */
+ for (i = 3; i <= ntot; i += 2) {
+ fact[i] = fact[i - 1] + log((double) i);
+ j = i + 1;
+ if (j <= ntot)
+ fact[j] = fact[i] + fact[2] + fact[j / 2] - fact[j / 2 - 1];
+ }
+ /* Compute obs := observed path length */
+ obs = tol;
+ ntot = 0;
+ for (j = 1; j <= nco; ++j) {
+ dd = 0.;
+ for (i = 1; i <= nro; ++i) {
+ if (nr_gt_nc) {
+ dd += fact[(int) table[j + i * *ldtabl]];
+ ntot += (int) table[j + i * *ldtabl];
+ } else {
+ dd += fact[(int) table[i + j * *ldtabl]];
+ ntot += (int) table[i + j * *ldtabl];
+ }
+ }
+ obs += fact[ico[j]] - dd;
+ }
+ /* Denominator of observed table: DRO */
+ dro = f9xact(&nro, &ntot, &iro[1], fact);
+ *prt = exp(obs - dro);
+ /* Initialize pointers */
+ k = nco;
+ last = *ldkey + 1;
+ jkey = *ldkey + 1;
+ jstp = *ldstp + 1;
+ jstp2 = *ldstp * 3 + 1;
+ jstp3 = (*ldstp << 2) + 1;
+ jstp4 = *ldstp * 5 + 1;
+ ikkey = 0;
+ ikstp = 0;
+ ikstp2 = *ldstp << 1;
+ ipo = 1;
+ ipoin[1] = 1;
+ stp[1] = 0.;
+ ifrq[1] = 1;
+ ifrq[ikstp2 + 1] = -1;
+
+Outer_Loop:
+ kb = nco - k + 1;
+ ks = 0;
+ n = ico[kb];
+ kd = nro + 1;
+ kmax = nro;
+ /* IDIF is the difference in going to the daughter */
+ for (i = 1; i <= nro; ++i)
+ idif[i] = 0;
+
+ /* Generate the first daughter */
+ do {
+ --kd;
+ ntot = min(n, iro[kd]);
+ idif[kd] = ntot;
+ if (idif[kmax] == 0)
+ --kmax;
+ n -= ntot;
+ }
+ while (n > 0 && kd != 1);
+
+ if (n != 0) {
+ goto L310;
+ }
+
+ k1 = k - 1;
+ n = ico[kb];
+ ntot = 0;
+ for (i = kb + 1; i <= nco; ++i)
+ ntot += ico[i];
+
+
+L150:
+ /* Arc to daughter length=ICO(KB) */
+ for (i = 1; i <= nro; ++i)
+ irn[i] = iro[i] - idif[i];
+
+ /* Sort irn */
+ if (k1 > 1) {
+ if (nro == 2) {
+ if (irn[1] > irn[2]) {
+ ii = irn[1];
+ irn[1] = irn[2];
+ irn[2] = ii;
+ }
+ } else if (nro == 3) {
+ ii = irn[1];
+ if (ii > irn[3]) {
+ if (ii > irn[2]) {
+ if (irn[2] > irn[3]) {
+ irn[1] = irn[3];
+ irn[3] = ii;
+ } else {
+ irn[1] = irn[2];
+ irn[2] = irn[3];
+ irn[3] = ii;
+ }
+ } else {
+ irn[1] = irn[3];
+ irn[3] = irn[2];
+ irn[2] = ii;
+ }
+ } else if (ii > irn[2]) {
+ irn[1] = irn[2];
+ irn[2] = ii;
+ } else if (irn[2] > irn[3]) {
+ ii = irn[2];
+ irn[2] = irn[3];
+ irn[3] = ii;
+ }
+ } else {
+ for (j = 2; j <= nro; ++j) {
+ i = j - 1;
+ ii = irn[j];
+
+ while (ii < irn[i]) {
+ irn[i + 1] = irn[i];
+ --i;
+ if (i == 0)
+ break;
+ }
+ irn[i + 1] = ii;
+ }
+ }
+ /* Adjust start for zero */
+ for (i = 1; i <= nro; ++i) {
+ if (irn[i] != 0)
+ break;
+ }
+
+ nrb = i;
+ nro2 = nro - i + 1;
+ } else {
+ nrb = 1;
+ nro2 = nro;
+ }
+ /* Some table values */
+ ddf = f9xact(&nro, &n, &idif[1], fact);
+ drn = f9xact(&nro2, &ntot, &irn[nrb], fact) - dro + ddf;
+ /* Get hash value */
+ if (k1 > 1) {
+ kval = irn[1] + irn[2] * kyy[2];
+ for (i = 3; i <= nro; ++i) {
+ kval += irn[i] * kyy[i];
+ }
+ /* Get hash table entry */
+ i = kval % (*ldkey << 1) + 1;
+ /* Search for unused location */
+ for (itp = i; itp <= *ldkey << 1; ++itp) {
+ ii = key2[itp];
+ if (ii == kval) {
+ goto L240;
+ } else if (ii < 0) {
+ key2[itp] = kval;
+ dlp[itp] = 1.;
+ dsp[itp] = 1.;
+ goto L240;
+ }
+ }
+
+ for (itp = 1; itp <= i - 1; ++itp) {
+ ii = key2[itp];
+ if (ii == kval) {
+ goto L240;
+ } else if (ii < 0) {
+ key2[itp] = kval;
+ dlp[itp] = 1.;
+ goto L240;
+ }
+ }
+
+ /* KH
+ prterr(6, "LDKEY is too small.\n"
+ "It is not possible to give the value of LDKEY required,\n"
+ "but you could try doubling LDKEY (and possibly LDSTP).");
+ */
+ prterr(6, "LDKEY is too small for this problem.\n"
+ "Try increasing the size of the workspace.");
+ }
+
+L240:
+ ipsh = (1);
+ /* Recover pastp */
+ ipn = ipoin[ipo + ikkey];
+ pastp = stp[ipn + ikstp];
+ ifreq = ifrq[ipn + ikstp];
+ /* Compute shortest and longest path */
+ if (k1 > 1) {
+ obs2 = obs - fact[ico[kb + 1]] - fact[ico[kb + 2]] - ddf;
+ for (i = 3; i <= k1; ++i) {
+ obs2 -= fact[ico[kb + i]];
+ }
+ if (dlp[itp] > 0.) {
+ dspt = obs - obs2 - ddf;
+ /* Compute longest path */
+ dlp[itp] = 0.;
+ f3xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dlp[itp],
+ &ntot, fact, &iwk[i31], &iwk[i32], &iwk[i33],
+ &iwk[i34], &iwk[i35], &iwk[i36], &iwk[i37],
+ &iwk[i38], &iwk[i39], &rwk[i310], &rwk[i311], &tol);
+ dlp[itp] = min(0., dlp[itp]);
+ /* Compute shortest path */
+ dsp[itp] = dspt;
+ f4xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dsp[itp], fact,
+ &iwk[i47], &iwk[i41], &iwk[i42], &iwk[i43],
+ &iwk[i44], &iwk[i45], &iwk[i46], &rwk[i48], &tol);
+ dsp[itp] = min(0., dsp[itp] - dspt);
+ /* Use chi-squared approximation? */
+ if ((irn[nrb] * ico[kb + 1]) > ntot * emn) {
+ ncell = 0.;
+ for (i = 0; i < nro2; ++i)
+ for (j = 1; j <= k1; ++j)
+ if (irn[nrb + i] * ico[kb + j] >= ntot * *expect)
+ ncell++;
+
+ if (ncell * 100 >= k1 * nro2 * *percnt) {
+ tmp = 0.;
+ for (i = 0; i < nro2; ++i)
+ tmp += (fact[irn[nrb + i]] -
+ fact[irn[nrb + i] - 1]);
+ tmp *= k1 - 1;
+ for (j = 1; j <= k1; ++j)
+ tmp += (nro2 - 1) * (fact[ico[kb + j]] -
+ fact[ico[kb + j] - 1]);
+ df = (double) ((nro2 - 1) * (k1 - 1));
+ tmp += df * 1.83787706640934548356065947281;
+ tmp -= (nro2 * k1 - 1) * (fact[ntot] - fact[ntot - 1]);
+ tm[itp] = (obs - dro) * -2. - tmp;
+ } else {
+ /* tm(itp) set to a flag value */
+ tm[itp] = -9876.;
+ }
+ } else {
+ tm[itp] = -9876.;
+ }
+ }
+ obs3 = obs2 - dlp[itp];
+ obs2 -= dsp[itp];
+ if (tm[itp] == -9876.) {
+ chisq = (0);
+ } else {
+ chisq = (1);
+ tmp = tm[itp];
+ }
+ } else {
+ obs2 = obs - drn - dro;
+ obs3 = obs2;
+ }
+
+L300:
+ /* Process node with new PASTP */
+ if (pastp <= obs3) {
+ /* Update pre */
+ *pre += (double) ifreq * exp(pastp + drn);
+ } else if (pastp < obs2) {
+ if (chisq) {
+ df = (double) ((nro2 - 1) * (k1 - 1));
+#ifdef USING_R
+ pv = pgamma(fmax2(0., tmp + (pastp + drn) * 2.) / 2.,
+ df / 2., /*scale = */ 1.,
+ /*lower_tail = */FALSE, /*log_p = */ FALSE);
+#else
+ d1 = max(0., tmp + (pastp + drn) * 2.) / 2.;
+ d2 = df / 2.;
+ pv = 1. - gammds(&d1, &d2, &ifault);
+#endif
+ *pre += (double) ifreq * exp(pastp + drn) * pv;
+ } else {
+ /* Put daughter on queue */
+ d1 = pastp + ddf;
+ f5xact(&d1, &tol, &kval, &key[jkey], ldkey, &ipoin[jkey],
+ &stp[jstp], ldstp, &ifrq[jstp], &ifrq[jstp2],
+ &ifrq[jstp3], &ifrq[jstp4], &ifreq, &itop, &ipsh);
+ ipsh = (0);
+ }
+ }
+ /* Get next PASTP on chain */
+ ipn = ifrq[ipn + ikstp2];
+ if (ipn > 0) {
+ pastp = stp[ipn + ikstp];
+ ifreq = ifrq[ipn + ikstp];
+ goto L300;
+ }
+ /* Generate a new daughter node */
+ f7xact(&kmax, &iro[1], &idif[1], &kd, &ks, &iflag);
+ if (iflag != 1) {
+ goto L150;
+ }
+
+L310:
+ /* Go get a new mother from stage K */
+ do {
+ iflag = 1;
+ f6xact(&nro, &iro[1], &iflag, &kyy[1], &key[ikkey + 1], ldkey,
+ &last, &ipo);
+ /* Update pointers */
+ if (iflag != 3)
+ goto Outer_Loop;
+ /* else iflag == 3 : no additional nodes to process */
+ --k;
+ itop = 0;
+ ikkey = jkey - 1;
+ ikstp = jstp - 1;
+ ikstp2 = jstp2 - 1;
+ jkey = *ldkey - jkey + 2;
+ jstp = *ldstp - jstp + 2;
+ jstp2 = (*ldstp << 1) + jstp;
+ for (i = 1; i <= *ldkey << 1; ++i)
+ key2[i] = -9999;
+
+ } while (k >= 2);
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F3XACT
+ Purpose: Computes the shortest path length for a given table.
+ Usage: F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, IRO,
+ IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, TOL)
+ Arguments:
+ NROW - The number of rows in the table. (Input)
+ IROW - Vector of length NROW containing the row sums
+ for the table. (Input)
+ NCOL - The number of columns in the table. (Input)
+ ICOL - Vector of length K containing the column sums
+ for the table. (Input)
+ DLP - The longest path for the table. (Output)
+ MM - The total count in the table. (Output)
+ FACT - Vector containing the logarithms of factorials. (Input)
+ ICO - Work vector of length MAX(NROW,NCOL).
+ IRO - Work vector of length MAX(NROW,NCOL).
+ IT - Work vector of length MAX(NROW,NCOL).
+ LB - Work vector of length MAX(NROW,NCOL).
+ NR - Work vector of length MAX(NROW,NCOL).
+ NT - Work vector of length MAX(NROW,NCOL).
+ NU - Work vector of length MAX(NROW,NCOL).
+ ITC - Work vector of length 400.
+ IST - Work vector of length 400.
+ STV - Work vector of length 400.
+ ALEN - Work vector of length MAX(NROW,NCOL).
+ TOL - Tolerance. (Input)
+ -----------------------------------------------------------------------
+ */
+
+void
+f3xact(int *nrow, int *irow, int *ncol, int *icol, double *dlp,
+ int *mm, double *fact, int *ico, int *iro, int *it,
+ int *lb, int *nr, int *nt, int *nu, int *itc, int *ist,
+ double *stv, double *alen, const double *tol)
+{
+ /* Initialized data */
+ static int ldst = 200;
+ static int nst = 0;
+ static int nitc = 0;
+
+ /* Local variables */
+ static int xmin;
+ static int i, k;
+ static double v;
+ static int n11, n12, ii, nn, ks, ic1, ic2, nc1, nn1;
+ static int nr1, nco;
+ static double val;
+ static int nct, ipn, irl, key, lev, itp, nro;
+ static double vmn;
+ static int nrt, kyy, nc1s;
+
+ /* Parameter adjustments */
+ --stv;
+ --ist;
+ --itc;
+ --nu;
+ --nt;
+ --nr;
+ --lb;
+ --it;
+ --iro;
+ --ico;
+ --icol;
+ --irow;
+
+ /* Function Body */
+ for (i = 0; i <= *ncol; ++i) {
+ alen[i] = 0.;
+ }
+ for (i = 1; i <= 400; ++i) {
+ ist[i] = -1;
+ }
+ /* nrow is 1 */
+ if (*nrow <= 1) {
+ if (*nrow > 0) {
+ *dlp -= fact[icol[1]];
+ for (i = 2; i <= *ncol; ++i) {
+ *dlp -= fact[icol[i]];
+ }
+ }
+ return;
+ }
+ /* ncol is 1 */
+ if (*ncol <= 1) {
+ if (*ncol > 0) {
+ *dlp = *dlp - fact[irow[1]] - fact[irow[2]];
+ for (i = 3; i <= *nrow; ++i) {
+ *dlp -= fact[irow[i]];
+ }
+ }
+ return;
+ }
+ /* 2 by 2 table */
+ if (*nrow * *ncol == 4) {
+ n11 = (irow[1] + 1) * (icol[1] + 1) / (*mm + 2);
+ n12 = irow[1] - n11;
+ *dlp = *dlp - fact[n11] - fact[n12] - fact[icol[1] - n11]
+ - fact[icol[2] - n12];
+ return;
+ }
+ /* Test for optimal table */
+ val = 0.;
+ xmin = (0);
+ if (irow[*nrow] <= irow[1] + *ncol) {
+ f10act(nrow, &irow[1], ncol, &icol[1], &val, &xmin, fact,
+ &lb[1], &nu[1], &nr[1]);
+ }
+ if (! xmin) {
+ if (icol[*ncol] <= icol[1] + *nrow) {
+ f10act(ncol, &icol[1], nrow, &irow[1], &val, &xmin, fact,
+ &lb[1], &nu[1], &nr[1]);
+ }
+ }
+
+ if (xmin) {
+ *dlp -= val;
+ return;
+ }
+ /* Setup for dynamic programming */
+ nn = *mm;
+ /* Minimize ncol */
+ if (*nrow >= *ncol) {
+ nro = *nrow;
+ nco = *ncol;
+ for (i = 1; i <= *nrow; ++i) {
+ iro[i] = irow[i];
+ }
+ ico[1] = icol[1];
+ nt[1] = nn - ico[1];
+ for (i = 2; i <= *ncol; ++i) {
+ ico[i] = icol[i];
+ nt[i] = nt[i - 1] - ico[i];
+ }
+ } else {
+ nro = *ncol;
+ nco = *nrow;
+ ico[1] = irow[1];
+ nt[1] = nn - ico[1];
+ for (i = 2; i <= *nrow; ++i) {
+ ico[i] = irow[i];
+ nt[i] = nt[i - 1] - ico[i];
+ }
+ for (i = 1; i <= *ncol; ++i)
+ iro[i] = icol[i];
+ }
+ /* Initialize pointers */
+ vmn = 1e10;
+ nc1s = nco - 1;
+ irl = 1;
+ ks = 0;
+ k = ldst;
+ kyy = ico[nco] + 1;
+
+LnewNode: /* Setup to generate new node */
+
+ lev = 1;
+ nr1 = nro - 1;
+ nrt = iro[irl];
+ nct = ico[1];
+ lb[1] = (int) ((double) ((nrt + 1) * (nct + 1)) /
+ (double) (nn + nr1 * nc1s + 1) - *tol) - 1;
+ nu[1] = (int) ((double) ((nrt + nc1s) * (nct + nr1)) /
+ (double) (nn + nr1 + nc1s)) - lb[1] + 1;
+ nr[1] = nrt - lb[1];
+
+LoopNode: /* Generate a node */
+ --nu[lev];
+ if (nu[lev] == 0) {
+ if (lev == 1)
+ goto L200;
+
+ --lev;
+ goto LoopNode;
+ }
+ ++lb[lev];
+ --nr[lev];
+L120:
+ alen[lev] = alen[lev - 1] + fact[lb[lev]];
+ if (lev < nc1s) {
+ nn1 = nt[lev];
+ nrt = nr[lev];
+ ++lev;
+ nc1 = nco - lev;
+ nct = ico[lev];
+ lb[lev] = (int) ((double) ((nrt + 1) * (nct + 1)) /
+ (double) (nn1 + nr1 * nc1 + 1) - *tol);
+ nu[lev] = (int) ((double) ((nrt + nc1) * (nct + nr1)) /
+ (double) (nn1 + nr1 + nc1) - lb[lev] + 1);
+ nr[lev] = nrt - lb[lev];
+ goto L120;
+ }
+ alen[nco] = alen[lev] + fact[nr[lev]];
+ lb[nco] = nr[lev];
+
+ v = val + alen[nco];
+ if (nro == 2) {
+ /* Only 1 row left */
+ v = v + fact[ico[1] - lb[1]] + fact[ico[2] - lb[2]];
+ for (i = 3; i <= nco; ++i) {
+ v += fact[ico[i] - lb[i]];
+ }
+ if (v < vmn) {
+ vmn = v;
+ }
+ } else if (nro == 3 && nco == 2) {
+ /* 3 rows and 2 columns */
+ nn1 = nn - iro[irl] + 2;
+ ic1 = ico[1] - lb[1];
+ ic2 = ico[2] - lb[2];
+ n11 = (iro[irl + 1] + 1) * (ic1 + 1) / nn1;
+ n12 = iro[irl + 1] - n11;
+ v = v + fact[n11] + fact[n12] + fact[ic1 - n11]
+ + fact[ic2 - n12];
+ if (v < vmn) {
+ vmn = v;
+ }
+ } else {
+ /* Column marginals are new node */
+ for (i = 1; i <= nco; ++i) {
+ it[i] = ico[i] - lb[i];
+ }
+ /* Sort column marginals */
+ if (nco == 2) {
+ if (it[1] > it[2]) {
+ ii = it[1];
+ it[1] = it[2];
+ it[2] = ii;
+ }
+ } else if (nco == 3) {
+ ii = it[1];
+ if (ii > it[3]) {
+ if (ii > it[2]) {
+ if (it[2] > it[3]) {
+ it[1] = it[3];
+ it[3] = ii;
+ } else {
+ it[1] = it[2];
+ it[2] = it[3];
+ it[3] = ii;
+ }
+ } else {
+ it[1] = it[3];
+ it[3] = it[2];
+ it[2] = ii;
+ }
+ } else if (ii > it[2]) {
+ it[1] = it[2];
+ it[2] = ii;
+ } else if (it[2] > it[3]) {
+ ii = it[2];
+ it[2] = it[3];
+ it[3] = ii;
+ }
+ } else {
+ isort(&nco, &it[1]);
+ }
+ /* Compute hash value */
+ key = it[1] * kyy + it[2];
+ for (i = 3; i <= nco; ++i) {
+ key = it[i] + key * kyy;
+ }
+ if(key < 0)
+ //PROBLEM "Bug in FEXACT: gave negative key" RECOVER(NULL_ENTRY);
+ printf("Bug in FEXACT: gave negative key \n"); //xx:another one of my ugly kluges (R.D-U)
+
+ /* Table index */
+ ipn = key % ldst + 1;
+
+ /* Find empty position */
+ for (itp = ipn, ii = ks + ipn; itp <= ldst; ++itp, ++ii) {
+ if (ist[ii] < 0) {
+ goto L180;
+ } else if (ist[ii] == key) {
+ goto L190;
+ }
+ }
+
+ for (itp = 1, ii = ks + 1; itp <= ipn - 1; ++itp, ++ii) {
+ if (ist[ii] < 0) {
+ goto L180;
+ } else if (ist[ii] == key) {
+ goto L190;
+ }
+ }
+
+ prterr(30, "Stack length exceeded in f3xact.\n"
+ "This problem should not occur.");
+
+L180: /* Push onto stack */
+ ist[ii] = key;
+ stv[ii] = v;
+ ++nst;
+ ii = nst + ks;
+ itc[ii] = itp;
+ goto LoopNode;
+
+L190: /* Marginals already on stack */
+ stv[ii] = min(v, stv[ii]);
+ }
+ goto LoopNode;
+
+
+L200: /* Pop item from stack */
+ if (nitc > 0) {
+ /* Stack index */
+ itp = itc[nitc + k] + k;
+ --nitc;
+ val = stv[itp];
+ key = ist[itp];
+ ist[itp] = -1;
+ /* Compute marginals */
+ for (i = nco; i >= 2; --i) {
+ ico[i] = key % kyy;
+ key /= kyy;
+ }
+ ico[1] = key;
+ /* Set up nt array */
+ nt[1] = nn - ico[1];
+ for (i = 2; i <= nco; ++i)
+ nt[i] = nt[i - 1] - ico[i];
+
+ /* Test for optimality (L90) */
+ xmin = (0);
+ if (iro[nro] <= iro[irl] + nco) {
+ f10act(&nro, &iro[irl], &nco, &ico[1], &val, &xmin, fact,
+ &lb[1], &nu[1], &nr[1]);
+ }
+ if (!xmin && ico[nco] <= ico[1] + nro)
+ f10act(&nco, &ico[1], &nro, &iro[irl], &val, &xmin, fact,
+ &lb[1], &nu[1], &nr[1]);
+ if (xmin) {
+ if (vmn > val)
+ vmn = val;
+ goto L200;
+ }
+ else goto LnewNode;
+
+ } else if (nro > 2 && nst > 0) {
+ /* Go to next level */
+ nitc = nst;
+ nst = 0;
+ k = ks;
+ ks = ldst - ks;
+ nn -= iro[irl];
+ ++irl;
+ --nro;
+ goto L200;
+ }
+
+ *dlp -= vmn;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F4XACT
+ Purpose: Computes the longest path length for a given table.
+ Usage: CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK,
+ NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK,
+ TOL)
+ Arguments:
+ NROW - The number of rows in the table. (Input)
+ IROW - Vector of length NROW containing the row sums for the
+ table. (Input)
+ NCOL - The number of columns in the table. (Input)
+ ICOL - Vector of length K containing the column sums for the
+ table. (Input)
+ DSP - The shortest path for the table. (Output)
+ FACT - Vector containing the logarithms of factorials. (Input)
+ ICSTK - NCOL by NROW+NCOL+1 work array.
+ NCSTK - Work vector of length NROW+NCOL+1.
+ LSTK - Work vector of length NROW+NCOL+1.
+ MSTK - Work vector of length NROW+NCOL+1.
+ NSTK - Work vector of length NROW+NCOL+1.
+ NRSTK - Work vector of length NROW+NCOL+1.
+ IRSTK - NROW by MAX(NROW,NCOL) work array.
+ YSTK - Work vector of length NROW+NCOL+1.
+ TOL - Tolerance. (Input)
+ -----------------------------------------------------------------------
+ */
+
+void
+f4xact(int *nrow, int *irow, int *ncol, int *icol, double *dsp,
+ double *fact, int *icstk, int *ncstk, int *lstk, int *mstk,
+ int *nstk, int *nrstk, int *irstk, double *ystk, const double *tol)
+{
+ /* System generated locals */
+ int ikh;
+
+ /* Local variables */
+ int i, j, k, l, m, n, mn, ic1, ir1, ict, irt, istk, nco, nro;
+ double y, amx;
+
+ /* Parameter adjustments */
+ irstk -= *nrow + 1;
+ --irow;
+ icstk -= *ncol + 1;
+ --icol;
+ --ncstk;
+ --lstk;
+ --mstk;
+ --nstk;
+ --nrstk;
+ --ystk;
+
+ /* Function Body */
+ /* Take care of the easy cases first */
+ if (*nrow == 1) {
+ for (i = 1; i <= *ncol; ++i) {
+ *dsp -= fact[icol[i]];
+ }
+ return;
+ }
+ if (*ncol == 1) {
+ for (i = 1; i <= *nrow; ++i) {
+ *dsp -= fact[irow[i]];
+ }
+ return;
+ }
+ if (*nrow * *ncol == 4) {
+ if (irow[2] <= icol[2]) {
+ *dsp = *dsp - fact[irow[2]] - fact[icol[1]]
+ - fact[icol[2] - irow[2]];
+ } else {
+ *dsp = *dsp - fact[icol[2]] - fact[irow[1]]
+ - fact[irow[2] - icol[2]];
+ }
+ return;
+ }
+ /* initialization before loop */
+ for (i = 1; i <= *nrow; ++i) {
+ irstk[i + *nrow] = irow[*nrow - i + 1];
+ }
+ for (j = 1; j <= *ncol; ++j) {
+ icstk[j + *ncol] = icol[*ncol - j + 1];
+ }
+
+ nro = *nrow;
+ nco = *ncol;
+ nrstk[1] = nro;
+ ncstk[1] = nco;
+ ystk[1] = 0.;
+ y = 0.;
+ istk = 1;
+ l = 1;
+ amx = 0.;
+
+ /* First LOOP */
+ do {
+ ir1 = irstk[istk * *nrow + 1];
+ ic1 = icstk[istk * *ncol + 1];
+ if (ir1 > ic1) {
+ if (nro >= nco) {
+ m = nco - 1; n = 2;
+ } else {
+ m = nro; n = 1;
+ }
+ } else if (ir1 < ic1) {
+ if (nro <= nco) {
+ m = nro - 1; n = 1;
+ } else {
+ m = nco; n = 2;
+ }
+ } else {
+ if (nro <= nco) {
+ m = nro - 1; n = 1;
+ } else {
+ m = nco - 1; n = 2;
+ }
+ }
+
+ L60:
+ if (n == 1) {
+ i = l; j = 1;
+ } else {
+ i = 1; j = l;
+ }
+
+ irt = irstk[i + istk * *nrow];
+ ict = icstk[j + istk * *ncol];
+ mn = irt;
+ if (mn > ict) {
+ mn = ict;
+ }
+ y += fact[mn];
+ if (irt == ict) {
+ --nro;
+ --nco;
+ f11act(&irstk[istk * *nrow + 1], &i, &nro,
+ &irstk[(istk + 1) * *nrow + 1]);
+ f11act(&icstk[istk * *ncol + 1], &j, &nco,
+ &icstk[(istk + 1) * *ncol + 1]);
+ } else if (irt > ict) {
+ --nco;
+ f11act(&icstk[istk * *ncol + 1], &j, &nco,
+ &icstk[(istk + 1) * *ncol + 1]);
+ ikh = irt - ict;
+ f8xact(&irstk[istk * *nrow + 1], &ikh, &i,
+ &nro, &irstk[(istk + 1) * *nrow + 1]);
+ } else {
+ --nro;
+ f11act(&irstk[istk * *nrow + 1], &i, &nro,
+ &irstk[(istk + 1) * *nrow + 1]);
+ ikh = ict - irt;
+ f8xact(&icstk[istk * *ncol + 1], &ikh, &j,
+ &nco, &icstk[(istk + 1) * *ncol + 1]);
+ }
+
+ if (nro == 1) {
+ for (k = 1; k <= nco; ++k) {
+ y += fact[icstk[k + (istk + 1) * *ncol]];
+ }
+ break;
+ }
+ if (nco == 1) {
+ for (k = 1; k <= nro; ++k) {
+ y += fact[irstk[k + (istk + 1) * *nrow]];
+ }
+ break;
+ }
+
+ lstk[istk] = l;
+ mstk[istk] = m;
+ nstk[istk] = n;
+ ++istk;
+ nrstk[istk] = nro;
+ ncstk[istk] = nco;
+ ystk[istk] = y;
+ l = 1;
+ } while(1);/* end do */
+
+/* L90:*/
+ if (y > amx) {
+ amx = y;
+ if (*dsp - amx <= *tol) {
+ *dsp = 0.;
+ return;
+ }
+ }
+
+L100:
+ --istk;
+ if (istk == 0) {
+ *dsp -= amx;
+ if (*dsp - amx <= *tol) {
+ *dsp = 0.;
+ }
+ return;
+ }
+ l = lstk[istk] + 1;
+
+/* L110: */
+ for(;; ++l) {
+ if (l > mstk[istk]) goto L100;
+
+ n = nstk[istk];
+ nro = nrstk[istk];
+ nco = ncstk[istk];
+ y = ystk[istk];
+ if (n == 1) {
+ if (irstk[l + istk * *nrow] <
+ irstk[l - 1 + istk * *nrow]) goto L60;
+ }
+ else if (n == 2) {
+ if (icstk[l + istk * *ncol] <
+ icstk[l - 1 + istk * *ncol]) goto L60;
+ }
+ }
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F5XACT
+ Purpose: Put node on stack in network algorithm.
+ Usage: CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP,
+ LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP,
+ IPSH)
+ Arguments:
+ PASTP - The past path length. (Input)
+ TOL - Tolerance for equivalence of past path lengths. (Input)
+ KVAL - Key value. (Input)
+ KEY - Vector of length LDKEY containing the key values. (in/out)
+ LDKEY - Length of vector KEY. (Input)
+ IPOIN - Vector of length LDKEY pointing to the
+ linked list of past path lengths. (in/out)
+ STP - Vector of length LSDTP containing the
+ linked lists of past path lengths. (in/out)
+ LDSTP - Length of vector STP. (Input)
+ IFRQ - Vector of length LDSTP containing the past path
+ frequencies. (in/out)
+ NPOIN - Vector of length LDSTP containing the pointers to
+ the next past path length. (in/out)
+ NR - Vector of length LDSTP containing the right object
+ pointers in the tree of past path lengths. (in/out)
+ NL - Vector of length LDSTP containing the left object
+ pointers in the tree of past path lengths. (in/out)
+ IFREQ - Frequency of the current path length. (Input)
+ ITOP - Pointer to the top of STP. (Input)
+ IPSH - Option parameter. (Input)
+ If IPSH is true, the past path length is found in the
+ table KEY. Otherwise the location of the past path
+ length is assumed known and to have been found in
+ a previous call. ==>>>>> USING "static" variables
+ -----------------------------------------------------------------------
+ */
+
+void
+f5xact(double *pastp, const double *tol, int *kval, int *key, int *ldkey,
+ int *ipoin, double *stp, int *ldstp, int *ifrq, int *npoin,
+ int *nr, int *nl, int *ifreq, int *itop, int *ipsh)
+{
+ /* Local variables */
+ static int itmp, ird, ipn, itp;
+ double test1, test2;
+
+ /* Parameter adjustments */
+ --nl;
+ --nr;
+ --npoin;
+ --ifrq;
+ --stp;
+ --ipoin;
+ --key;
+
+ /* Function Body */
+ if (*ipsh) {
+ /* Convert KVAL to int in range 1, ..., LDKEY. */
+ ird = *kval % *ldkey + 1;
+ /* Search for an unused location */
+ for (itp = ird; itp <= *ldkey; ++itp) {
+ if (key[itp] == *kval) {
+ goto L40;
+ }
+ if (key[itp] < 0) {
+ goto L30;
+ }
+ }
+ for (itp = 1; itp <= ird - 1; ++itp) {
+ if (key[itp] == *kval) {
+ goto L40;
+ }
+ if (key[itp] < 0) {
+ goto L30;
+ }
+ }
+ /* Return if KEY array is full */
+ /* KH
+ prterr(6, "LDKEY is too small for this problem.\n"
+ "It is not possible to estimate the value of LDKEY "
+ "required,\n"
+ "but twice the current value may be sufficient.");
+ */
+ prterr(6, "LDKEY is too small for this problem.\n"
+ "Try increasing the size of the workspace.");
+
+ /* Update KEY */
+L30:
+ key[itp] = *kval;
+ ++(*itop);
+ ipoin[itp] = *itop;
+ /* Return if STP array full */
+ if (*itop > *ldstp) {
+ /* KH
+ prterr(7, "LDSTP is too small for this problem.\n"
+ "It is not possible to estimate the value of LDSTP "
+ "required,\n"
+ "but twice the current value may be sufficient.");
+ */
+ prterr(7, "LDSTP is too small for this problem.\n"
+ "Try increasing the size of the workspace.");
+ }
+ /* Update STP, etc. */
+ npoin[*itop] = -1;
+ nr[*itop] = -1;
+ nl[*itop] = -1;
+ stp[*itop] = *pastp;
+ ifrq[*itop] = *ifreq;
+ return;
+ }
+
+ /* Find location, if any, of pastp */
+L40:
+ ipn = ipoin[itp];
+ test1 = *pastp - *tol;
+ test2 = *pastp + *tol;
+
+L50:
+ if (stp[ipn] < test1) {
+ ipn = nl[ipn];
+ if (ipn > 0) {
+ goto L50;
+ }
+ } else if (stp[ipn] > test2) {
+ ipn = nr[ipn];
+ if (ipn > 0) {
+ goto L50;
+ }
+ } else {
+ ifrq[ipn] += *ifreq;
+ return;
+ }
+ /* Return if STP array full */
+ ++(*itop);
+ if (*itop > *ldstp) {
+ /*
+ prterr(7, "LDSTP is too small for this problem.\n"
+ "It is not possible to estimate the value of LDSTP "
+ "required,\n"
+ "but twice the current value may be sufficient.");
+ */
+ prterr(7, "LDSTP is too small for this problem.\n"
+ "Try increasing the size of the workspace.");
+ return;
+ }
+ /* Find location to add value */
+ ipn = ipoin[itp];
+ itmp = ipn;
+
+L60:
+ if (stp[ipn] < test1) {
+ itmp = ipn;
+ ipn = nl[ipn];
+ if (ipn > 0) {
+ goto L60;
+ } else {
+ nl[itmp] = *itop;
+ }
+ } else if (stp[ipn] > test2) {
+ itmp = ipn;
+ ipn = nr[ipn];
+ if (ipn > 0) {
+ goto L60;
+ } else {
+ nr[itmp] = *itop;
+ }
+ }
+ /* Update STP, etc. */
+ npoin[*itop] = npoin[itmp];
+ npoin[itmp] = *itop;
+ stp[*itop] = *pastp;
+ ifrq[*itop] = *ifreq;
+ nl[*itop] = -1;
+ nr[*itop] = -1;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F6XACT
+ Purpose: Pop a node off the stack.
+ Usage: CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, IPN)
+ Arguments:
+ NROW - The number of rows in the table. (Input)
+ IROW - Vector of length nrow containing the row sums on
+ output. (Output)
+ IFLAG - Set to 3 if there are no additional nodes to process.
+ (Output)
+ KYY - Constant mutlipliers used in forming the hash
+ table key. (Input)
+ KEY - Vector of length LDKEY containing the hash table
+ keys. (In/out)
+ LDKEY - Length of vector KEY. (Input)
+ LAST - Index of the last key popped off the stack. (In/out)
+ IPN - Pointer to the linked list of past path lengths. (Output)
+ -----------------------------------------------------------------------
+ */
+void
+f6xact(int *nrow, int *irow, int *iflag, int *kyy, int *key, int
+ *ldkey, int *last, int *ipn)
+{
+ int kval, j;
+
+ /* Parameter adjustments */
+ --key;
+ --kyy;
+ --irow;
+
+ /* Function Body */
+L10:
+ ++(*last);
+ if (*last <= *ldkey) {
+ if (key[*last] < 0) {
+ goto L10;
+ }
+ /* Get KVAL from the stack */
+ kval = key[*last];
+ key[*last] = -9999;
+ for (j = *nrow; j >= 2; --j) {
+ irow[j] = kval / kyy[j];
+ kval -= irow[j] * kyy[j];
+ }
+ irow[1] = kval;
+ *ipn = *last;
+ } else {
+ *last = 0;
+ *iflag = 3;
+ }
+ return;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F7XACT
+ Purpose: Generate the new nodes for given marinal totals.
+ Usage: CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG)
+ Arguments:
+ NROW - The number of rows in the table. (Input)
+ IMAX - The row marginal totals. (Input)
+ IDIF - The column counts for the new column. (in/out)
+ K - Indicator for the row to decrement. (in/out)
+ KS - Indicator for the row to increment. (in/out)
+ IFLAG - Status indicator. (Output)
+ If IFLAG is zero, a new table was generated. For
+ IFLAG = 1, no additional tables could be generated.
+ -----------------------------------------------------------------------
+ */
+
+void
+f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
+ int *iflag)
+
+{
+ int i, m, k1, mm;
+
+ /* Parameter adjustments */
+ --idif;
+ --imax;
+
+ /* Function Body */
+ *iflag = 0;
+ /* Find node which can be incremented, ks */
+ if (*ks == 0)
+ do {
+ ++(*ks);
+ } while (idif[*ks] == imax[*ks]);
+
+ /* Find node to decrement (>ks) */
+ if (idif[*k] > 0 && *k > *ks) {
+ --idif[*k];
+ do {
+ --(*k);
+ } while(imax[*k] == 0);
+
+ m = *k;
+
+ /* Find node to increment (>=ks) */
+ while (idif[m] >= imax[m]) {
+ --m;
+ }
+ ++idif[m];
+ /* Change ks */
+ if (m == *ks) {
+ if (idif[m] == imax[m]) {
+ *ks = *k;
+ }
+ }
+ }
+ else {
+ Loop:
+ /* Check for finish */
+ for (k1 = *k + 1; k1 <= *nrow; ++k1) {
+ if (idif[k1] > 0) {
+ goto L70;
+ }
+ }
+ *iflag = 1;
+ return;
+
+ L70:
+ /* Reallocate counts */
+ mm = 1;
+ for (i = 1; i <= *k; ++i) {
+ mm += idif[i];
+ idif[i] = 0;
+ }
+ *k = k1;
+
+ do {
+ --(*k);
+ m = min(mm, imax[*k]);
+ idif[*k] = m;
+ mm -= m;
+ } while (mm > 0 && *k != 1);
+
+ /* Check that all counts reallocated */
+ if (mm > 0) {
+ if (k1 != *nrow) {
+ *k = k1;
+ goto Loop;
+ }
+ *iflag = 1;
+ return;
+ }
+ /* Get ks */
+ --idif[k1];
+ *ks = 0;
+ do {
+ ++(*ks);
+ if (*ks > *k) {
+ return;
+ }
+ } while (idif[*ks] >= imax[*ks]);
+ }
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F8XACT
+ Purpose: Routine for reducing a vector when there is a zero
+ element.
+ Usage: CALL F8XACT (IROW, IS, I1, IZERO, NEW)
+ Arguments:
+ IROW - Vector containing the row counts. (Input)
+ IS - Indicator. (Input)
+ I1 - Indicator. (Input)
+ IZERO - Position of the zero. (Input)
+ NEW - Vector of new row counts. (Output)
+ -----------------------------------------------------------------------
+ */
+
+void
+f8xact(int *irow, int *is, int *i1, int *izero, int *new)
+{
+ int i;
+
+ /* Parameter adjustments */
+ --new;
+ --irow;
+
+ /* Function Body */
+ for (i = 1; i < *i1; ++i)
+ new[i] = irow[i];
+
+ for (i = *i1; i <= *izero - 1; ++i) {
+ if (*is >= irow[i + 1])
+ break;
+ new[i] = irow[i + 1];
+ }
+
+ new[i] = *is;
+
+ for(;;) {
+ ++i;
+ if (i > *izero)
+ return;
+ new[i] = irow[i];
+ }
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F9XACT
+ Purpose: Computes the log of a multinomial coefficient.
+ Usage: F9XACT(N, MM, IR, FACT)
+ Arguments:
+ N - Length of IR. (Input)
+ MM - Number for factorial in numerator. (Input)
+ IR - Vector of length N containing the numbers for
+ the denominator of the factorial. (Input)
+ FACT - Table of log factorials. (Input)
+ F9XACT - The log of the multinomal coefficient. (Output)
+ -----------------------------------------------------------------------
+ */
+
+double
+f9xact(int *n, int *mm, int *ir, double *fact)
+{
+ double d;
+ int k;
+
+ d = fact[*mm];
+ for (k = 0; k < *n; k++)
+ d -= fact[ir[k]];
+ return d;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F10ACT
+ Purpose: Computes the shortest path length for special tables.
+ Usage: F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, NE, M)
+ Arguments:
+ NROW - The number of rows in the table. (Input)
+ IROW - Vector of length NROW containing the row totals. (Input)
+ NCOL - The number of columns in the table. (Input)
+ ICO - Vector of length NCOL containing the column totals.(Input)
+ VAL - The shortest path. (Output)
+ XMIN - Set to true if shortest path obtained. (Output)
+ FACT - Vector containing the logarithms of factorials. (Input)
+ ND - Workspace vector of length NROW. (Input)
+ NE - Workspace vector of length NCOL. (Input)
+ M - Workspace vector of length NCOL. (Input)
+
+ Chapter: STAT/LIBRARY Categorical and Discrete Data Analysis
+ -----------------------------------------------------------------------
+ */
+
+void
+f10act(int *nrow, int *irow, int *ncol, int *icol, double *val,
+ int *xmin, double *fact, int *nd, int *ne, int *m)
+{
+ /* Local variables */
+ int i, is, ix, nrw1;
+
+ /* Parameter adjustments */
+ --m;
+ --ne;
+ --nd;
+ --icol;
+ --irow;
+
+ /* Function Body */
+ for (i = 1; i <= *nrow - 1; ++i)
+ nd[i] = 0;
+
+ is = icol[1] / *nrow;
+ ix = icol[1] - *nrow * is;
+ ne[1] = is;
+ m[1] = ix;
+ if (ix != 0)
+ ++nd[ix];
+
+ for (i = 2; i <= *ncol; ++i) {
+ ix = icol[i] / *nrow;
+ ne[i] = ix;
+ is += ix;
+ ix = icol[i] - *nrow * ix;
+ m[i] = ix;
+ if (ix != 0)
+ ++nd[ix];
+ }
+
+ for (i = *nrow - 2; i >= 1; --i)
+ nd[i] += nd[i + 1];
+
+ ix = 0;
+ nrw1 = *nrow + 1;
+ for (i = *nrow; i >= 2; --i) {
+ ix = ix + is + nd[nrw1 - i] - irow[i];
+ if (ix < 0)
+ return;
+ }
+
+ for (i = 1; i <= *ncol; ++i) {
+ ix = ne[i];
+ is = m[i];
+ *val = *val + is * fact[ix + 1] + (*nrow - is) * fact[ix];
+ }
+ *xmin = (1);
+
+ return;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: F11ACT
+ Purpose: Routine for revising row totals.
+ Usage: CALL F11ACT (IROW, I1, I2, NEW)
+ Arguments:
+ IROW - Vector containing the row totals. (Input)
+ I1 - Indicator. (Input)
+ I2 - Indicator. (Input)
+ NEW - Vector containing the row totals. (Output)
+ -----------------------------------------------------------------------
+ */
+void
+f11act(int *irow, int *i1, int *i2, int *new)
+{
+ int i;
+
+ /* Parameter adjustments */
+ --new;
+ --irow;
+
+ for (i = 1; i <= (*i1 - 1); ++i) new[i] = irow[i];
+ for (i = *i1; i <= *i2; ++i) new[i] = irow[i + 1];
+
+ return;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: prterr
+ Purpose: Print an error message and stop.
+ Usage: prterr(icode, mes)
+ Arguments:
+ icode - Integer code for the error message. (Input)
+ mes - Character string containing the error message. (Input)
+ -----------------------------------------------------------------------
+ */
+void
+prterr(int icode, char *mes)
+{
+// PROBLEM "FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY);
+// printf("FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY));
+ printf("FEXACT error %d.\n", icode); //xx:another one of my ugly kluges
+ return;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: iwork
+ Purpose: Routine for allocating workspace.
+ Usage: iwork (iwkmax, iwkpt, number, itype)
+ Arguments:
+ iwkmax - Maximum (int) amount of workspace. (Input)
+ iwkpt - Amount of (int) workspace currently allocated. (in/out)
+ number - Number of elements of workspace desired. (Input)
+ itype - Workspace type. (Input)
+ ITYPE TYPE
+ 2 integer
+ 3 float
+ 4 double
+ iwork(): Index in rwrk, dwrk, or iwrk of the beginning of
+ the first free element in the workspace array. (Output)
+ -----------------------------------------------------------------------
+ */
+int
+iwork(int iwkmax, int *iwkpt, int number, int itype)
+{
+ int i;
+
+ i = *iwkpt;
+ if (itype == 2 || itype == 3)
+ *iwkpt += number;
+ else { /* double */
+ if (i % 2 != 0)
+ ++i;
+ *iwkpt += (number << 1);
+ i /= 2;
+ }
+ if (*iwkpt > iwkmax)
+ prterr(40, "Out of workspace.");
+
+ return i;
+}
+
+#ifndef USING_R
+
+void isort(int *n, int *ix)
+{
+/*
+ -----------------------------------------------------------------------
+ Name: ISORT
+ Purpose: Shell sort for an int vector.
+ Usage: CALL ISORT (N, IX)
+ Arguments:
+ N - Lenth of vector IX. (Input)
+ IX - Vector to be sorted. (in/out)
+ -----------------------------------------------------------------------
+ */
+ static int ikey, i, j, m, il[10], kl, it, iu[10], ku;
+
+ /* Parameter adjustments */
+ --ix;
+
+ /* Function Body */
+ m = 1;
+ i = 1;
+ j = *n;
+
+L10:
+ if (i >= j) {
+ goto L40;
+ }
+ kl = i;
+ ku = j;
+ ikey = i;
+ ++j;
+ /* Find element in first half */
+L20:
+ ++i;
+ if (i < j) {
+ if (ix[ikey] > ix[i]) {
+ goto L20;
+ }
+ }
+ /* Find element in second half */
+L30:
+ --j;
+ if (ix[j] > ix[ikey]) {
+ goto L30;
+ }
+ /* Interchange */
+ if (i < j) {
+ it = ix[i];
+ ix[i] = ix[j];
+ ix[j] = it;
+ goto L20;
+ }
+ it = ix[ikey];
+ ix[ikey] = ix[j];
+ ix[j] = it;
+ /* Save upper and lower subscripts of the array yet to be sorted */
+ if (m < 11) {
+ if (j - kl < ku - j) {
+ il[m - 1] = j + 1;
+ iu[m - 1] = ku;
+ i = kl;
+ --j;
+ } else {
+ il[m - 1] = kl;
+ iu[m - 1] = j - 1;
+ i = j + 1;
+ j = ku;
+ }
+ ++m;
+ goto L10;
+ } else {
+ prterr(20, "This should never occur.");
+ }
+ /* Use another segment */
+L40:
+ --m;
+ if (m == 0) {
+ return;
+ }
+ i = il[m - 1];
+ j = iu[m - 1];
+ goto L10;
+}
+
+double gammds(double *y, double *p, int *ifault)
+{
+/*
+ -----------------------------------------------------------------------
+ Name: GAMMDS
+ Purpose: Cumulative distribution for the gamma distribution.
+ Usage: PGAMMA (Q, ALPHA,IFAULT)
+ Arguments:
+ Q - Value at which the distribution is desired. (Input)
+ ALPHA - Parameter in the gamma distribution. (Input)
+ IFAULT - Error indicator. (Output)
+ IFAULT DEFINITION
+ 0 No error
+ 1 An argument is misspecified.
+ 2 A numerical error has occurred.
+ PGAMMA - The cdf for the gamma distribution with parameter alpha
+ evaluated at Q. (Output)
+ -----------------------------------------------------------------------
+
+ Algorithm AS 147 APPL. Statist. (1980) VOL. 29, P. 113
+
+ Computes the incomplete gamma integral for positive parameters Y, P
+ using and infinite series.
+ */
+
+ static double a, c, f, g;
+ static int ifail;
+
+ /* Checks for the admissibility of arguments and value of F */
+ *ifault = 1;
+ g = 0.;
+ if (*y <= 0. || *p <= 0.) {
+ return g;
+ }
+ *ifault = 2;
+
+ /*
+ ALOGAM is natural log of gamma function no need to test ifail as
+ an error is impossible
+ */
+
+ a = *p + 1.;
+ f = exp(*p * log(*y) - alogam(&a, &ifail) - *y);
+ if (f == 0.) {
+ return g;
+ }
+ *ifault = 0;
+
+ /* Series begins */
+ c = 1.;
+ g = 1.;
+ a = *p;
+L10:
+ a += 1.;
+ c = c * *y / a;
+ g += c;
+ if (c / g > 1e-6) {
+ goto L10;
+ }
+ g *= f;
+ return g;
+}
+
+/*
+ -----------------------------------------------------------------------
+ Name: ALOGAM
+ Purpose: Value of the log-gamma function.
+ Usage: ALOGAM (X, IFAULT)
+ Arguments:
+ X - Value at which the log-gamma function is to be evaluated.
+ (Input)
+ IFAULT - Error indicator. (Output)
+ IFAULT DEFINITION
+ 0 No error
+ 1 X < 0
+ ALGAMA - The value of the log-gamma function at XX. (Output)
+ -----------------------------------------------------------------------
+
+ Algorithm ACM 291, Comm. ACM. (1966) Vol. 9, P. 684
+
+ Evaluates natural logarithm of gamma(x) for X greater than zero.
+ */
+
+double alogam(double *x, int *ifault)
+{
+ /* Initialized data */
+
+ static double a1 = .918938533204673;
+ static double a2 = 5.95238095238e-4;
+ static double a3 = 7.93650793651e-4;
+ static double a4 = .002777777777778;
+ static double a5 = .083333333333333;
+
+ /* Local variables */
+ static double f, y, z;
+
+ *ifault = 1;
+ if (*x < 0.) {
+ return(0.);
+ }
+ *ifault = 0;
+ y = *x;
+ f = 0.;
+ if (y >= 7.) {
+ goto L30;
+ }
+ f = y;
+L10:
+ y += 1.;
+ if (y >= 7.) {
+ goto L20;
+ }
+ f *= y;
+ goto L10;
+L20:
+ f = -log(f);
+L30:
+ z = 1. / (y * y);
+ return(f + (y - .5) * log(y) - y + a1 +
+ (((-a2 * z + a3) * z - a4) * z + a5) / y);
+}
+
+#endif /* not USING_R */
+
--- /dev/null
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <time.h>
+#include <math.h>
+#include "fisher2.h"
+
+void testp(double *permuted_ttests,int *B,double *permuted,double
+ *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double *ps);
+void permute_matrix(double *Imatrix,int *nc,int *nr,double
+ *permuted,int *g,double *trial_ts,double *Tinitial,double
+ *counter);
+void permute_array(int *array, int n);
+void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,double
+ *Ts,double *Tinitial,double *counter1);
+void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *storage);
+void start(double *Imatrix,int *g,int *nr,int *nc,double *testing,
+ double storage[][9]);
+
+int main (int argc, char *argv[]){
+
+ int col=-1, row=-1, size,g=0,c=0,i=0,j=0,k,counter=0,lines=0, B=10000;
+ double placeholder=0,min=0, thresh=0.05;
+
+ char arr[10000], str[51], a;
+ char location[41]="jobj.txt", output[41]="out.txt";
+
+ for(i=0;i<10000;i++){
+ arr[i]='q';
+ }
+
+int u,rflag=0,cflag=0, bflag=0;
+char *filename;
+int numbers;
+double numb;
+extern char *optarg;
+extern int optind, optopt, opterr;
+
+while ((u = getopt(argc, argv, ":r:c:g:b:t:f:o:")) != -1) {
+ switch(u) {
+ case 'r':
+ numbers = atoi(optarg);
+ printf("The number of features/rows is %d.\n", numbers);
+ row = numbers;
+ rflag = 1;
+ break;
+ case 'c':
+ numbers = atoi(optarg);
+ printf("The number of samples/columns is %d.\n", numbers);
+ col = numbers;
+ cflag = 1;
+ break;
+ case 'g':
+ numbers = atoi(optarg);
+ printf("Your g-value is %d.\n", numbers);
+ g = numbers;
+ break;
+ case 'b':
+ numbers = atoi(optarg);
+ printf("The number of permutations is %d\n", numbers);
+ B = numbers;
+ break;
+ case 't':
+ numb = atof(optarg);
+ printf("Threshold is is %lf\n", numb);
+ thresh = numb;
+ break;
+ case 'f':
+ filename = optarg;
+ printf("filename input is %s\n", filename);
+ strcpy(location,filename);
+ break;
+ case 'o':
+ filename = optarg;
+ printf("filename output %s\n", filename);
+ strcpy(output,filename);
+ break;
+ case ':':
+ printf("-%c without filename\n", optopt);
+ break;
+ case '?':
+ printf("unknown arg %c\n", optopt);
+ break;
+ }
+}
+
+ FILE *jobj, *out;
+ jobj=fopen(location,"r");
+
+ if(jobj == NULL){
+ printf("Please don't forget to save your matrix in the active");
+ printf(" directory as \"%s\".\n",location);
+ return 0;
+ }
+
+ // Gets the first line of samples names and checks for user error.
+ fgets(arr,10000,jobj);
+
+ for(i=0;i<10000;i++){
+ if(isspace(arr[i])){
+ counter++; }
+ }
+
+ if (cflag == 0) {
+ printf("You didn't tell us how many subjects there are!\n");
+ printf("But we'll still do the analysis as if there are %d subjects.\n\n",col=counter-1);
+ }
+ if (cflag == 1) {
+ if (col != counter-1){
+ printf("We would expect %d subjects, but you said %d.\n",counter-1,col);
+ }
+ }
+
+ while((a = fgetc(jobj)) != EOF){
+ if(a == '\n'){
+ lines++; }
+ }
+
+ if (rflag == 0) {
+ printf("You didn't specify the number of features!\n");
+ printf("But we'll still do the analysis assuming %d features.\n\n", row=lines-1);
+ }
+ if (rflag == 1) {
+ if ( lines-1 != row ){
+ printf("We would expect %d features, but you said %d.\n",lines-1,row);
+ }
+ }
+
+ if (g>=col || g<=0){
+ printf("Check your g value\n");
+ }
+
+ // Initialize the matrices
+ size = row*col;
+ double matrix[row][col];
+ double pmatrix[size],pmatrix2[size],permuted[size];
+ double storage[row][9];
+
+ for (i=0;i<row;i++){
+ for (j =0;j<9;j++){
+ storage[i][j]=0;
+ }
+ }
+ // Reset file below and create a separate matrix.
+ rewind(jobj);
+ fgets(arr,10000,jobj);
+
+ for(i=0; i<row; i++){
+ fscanf(jobj,"%s",str);
+ for(j=0; j<col;j++){
+ fscanf(jobj,"%lf",&placeholder);
+ matrix[i][j]=placeholder;
+ if(isalnum(placeholder)!=0){ // check for ""
+ printf("Your matrix isn't set up properly!\n");
+ return 0;
+ }
+ pmatrix[c]=0; // initializing to zero
+ permuted[c]=0;
+ c++;
+ }
+ }
+
+ fclose(jobj);
+
+
+ // Produces the sum of each column
+ double total[col],total1=0,total2=0;
+ double ratio[col];
+
+ for(i=0;i<col;i++){
+ total[i]=0;
+ ratio[i]=0; }
+
+ for(i=0; i<col; i++){
+ for(j=0;j<row;j++){
+ total[i]=total[i]+matrix[j][i];
+ }
+ }
+
+ for(i=0;i<g-1;i++){
+ total1=total1+total[i];}
+
+ for(i=g-1;i<col;i++){
+ total2=total2+total[i];}
+
+
+ // Creates the ratios by first finding the minimum of totals
+ min = total[0];
+ if (col==2){
+ if (total[0]<total[1]){
+ min = total[1];}
+ }
+ if (col >2){
+ for(i=1;i<col;i++){
+ if (min > total[i]){
+ min = total[i];}
+ }
+ }
+ if (min<=0){
+ printf("Error, the sum of one of the columns <= 0.");
+ return 0;
+ }
+
+
+ // Ratio time...
+ for(i=0;i<col;i++){
+ ratio[i]=total[i]/min;
+ }
+
+ // Change matrix into an array as received by R for compatibility.
+
+ c=0;
+ for(i=0;i<col;i++){
+ for(j=0;j<row;j++){
+ pmatrix[c]=matrix[j][i];
+ c++;
+ }
+ }
+
+ if(row == 1){
+ for (i =0; i<col;i++){
+ pmatrix[i]=pmatrix[i]/ratio[i];
+ }
+ }
+ else {
+ counter = 0;
+ j=-1;
+ for (i=0; i<size; i++) {
+ if (counter % row == 0) {
+ j++;
+ }
+ pmatrix[i]=pmatrix[i]/ratio[j];
+ counter++;
+ }
+ }
+ // pass everything to the rest of the code using pointers. then
+ // write to output file. below pointers for most of the values are
+ // created to send everything by reference.
+
+ int ptt_size, *permutes,*nc,*nr,*gvalue;
+
+ nc = &col;
+ nr = &row;
+ gvalue = &g;
+
+ permutes = &B;
+ ptt_size = B*row;
+
+ //changing ptt_size to row
+ double permuted_ttests[row], pvalues[row], tinitial[row];
+
+ for(i=0;i<row;i++){
+ permuted_ttests[i]=0;}
+
+ for(i=0;i<row;i++){
+ pvalues[i]=0;
+ tinitial[i]=0; }
+
+ // Find the initial values for the matrix.
+ start(pmatrix,gvalue,nr,nc,tinitial,storage);
+
+ // Start the calculations.
+
+ if ( (col==2) || ((g-1)<8) || ((col-g+1) < 8) ){
+
+ double fish[row], fish2[row];
+ for(i=0;i<row;i++){
+ fish[i]=0;
+ fish2[i]=0;}
+
+ for(i=0;i<row;i++){
+
+ for(j=0;j<g-1;j++){
+ fish[i]=fish[i]+matrix[i][j];
+ }
+
+ for(j=g-1;j<col;j++){
+ fish2[i]=fish2[i]+matrix[i][j];
+ }
+
+ double f11,f12,f21,f22;
+
+ f11=fish[i];
+ f12=fish2[i];
+
+ f21=total1-f11;
+ f22=total2-f12;
+
+ double data[] = {f11, f12, f21, f22};
+
+ // CONTINGENGCY TABLE:
+ // f11 f12
+ // f21 f22
+
+ int *nr, *nc, *ldtabl, *work;
+ int nrow=2, ncol=2, ldtable=2, workspace=100000;
+ double *expect, *prc, *emin,*prt,*pre;
+ double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
+
+ nr = &nrow;
+ nc = &ncol;
+ ldtabl=&ldtable;
+ work = &workspace;
+
+ expect = &e;
+ prc = &prc1;
+ emin=&emin1;
+ prt=&prt1;
+ pre=&pre1;
+
+ fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
+
+ if (*pre>.999999999){
+ *pre=1;
+ }
+ storage[i][8] = *pre;
+ pvalues[i]=*pre;
+ }
+ }
+ else{
+
+ testp(permuted_ttests, permutes, permuted,pmatrix, nc, nr, gvalue,tinitial,pvalues);
+
+ // Checks to make sure the matrix isn't sparse.
+ double sparse[row], sparse2[row];
+ for(i=0;i<row;i++){
+ sparse[i]=0;
+ sparse2[i]=0;}
+
+ c=0;
+ for(i=0;i<row;i++){
+
+ for(j=0;j<g-1;j++){
+ sparse[i]=sparse[i]+matrix[i][j];
+ }
+
+ if(sparse[i] < (double)(g-1)){
+ c++;
+ }
+ for(j=g-1;j<col;j++){ // ?<= for col
+ sparse2[i]=sparse2[i]+matrix[i][j];
+ }
+
+ if( (sparse2[i] <(double)(col-g+1))) {
+ c++;
+ }
+
+ if (c==2){
+ c=0;
+
+ double f11,f12,f21,f22;
+
+ f11=sparse[i];
+ sparse[i]=0;
+
+ f12=sparse2[i];
+ sparse2[i]=0;
+
+ f21=total1-f11;
+ f22=total2-f12;
+
+ double data[] = {f11, f12, f21, f22};
+
+ int *nr, *nc, *ldtabl, *work;
+ int nrow=2, ncol=2, ldtable=2, workspace=10000000; // I added two zeros for larger data sets
+ double *expect, *prc, *emin,*prt,*pre;
+ double e=0, prc1=0, emin1=0, prt1=0, pre1=0;
+
+ nr = &nrow;
+ nc = &ncol;
+ ldtabl=&ldtable;
+ work = &workspace;
+
+ expect = &e;
+ prc = &prc1;
+ emin=&emin1;
+ prt=&prt1;
+ pre=&pre1;
+
+ fexact(nr,nc,data, ldtabl,expect,prc,emin,prt,pre,work);
+
+ if (*pre>.999999999){
+ *pre=1;
+ }
+ storage[i][8] = *pre;
+ pvalues[i]=*pre;
+ }
+ }
+ // End of else statement
+ bflag = 1;
+ }
+
+ // Calculates the mean of counts (not normalized)
+ double temp[row][2];
+
+ for(j=0;j<row;j++){
+ for(i=0;i<2;i++){
+ temp[j][i]=0;
+ }
+ }
+
+ for (j=0;j<row;j++){
+ for (i=1; i<=(g-1); i++){
+ temp[j][0]=temp[j][0]+matrix[j][i-1];
+ }
+ temp[j][0]= (double) temp[j][0]/(g-1);
+ for(i=g;i<=col;i++){
+ temp[j][1]=temp[j][1]+matrix[j][i-1];
+ }
+ temp[j][1]= (double) temp[j][1]/(col-g+1);
+ }
+
+ for(i=0;i<row;i++){
+ storage[i][3]=temp[i][0];
+ storage[i][7]=temp[i][1];
+ storage[i][8]=pvalues[i];
+ }
+
+// BACKUP checks
+
+ for (i=0;i<row;i++){
+ if(pvalues[i]<thresh){
+ printf("Feature %d is significant, p = %.10lf \n",i+1,pvalues[i]);
+ }
+ }
+
+ // And now we write the files to a text file.
+ struct tm *local;
+ time_t t;
+ t = time(NULL);
+ local = localtime(&t);
+
+ jobj= fopen(location,"r");
+ fgets(arr,10000,jobj);
+
+ out = fopen(output,"a+");
+
+ fprintf(out,"Local time and date of test: %s\n", asctime(local));
+ fprintf(out,"# rows = %d, # col = %d, g = %d\n\n",row,col,g);
+ if (bflag == 1){
+ fprintf(out,"%d permutations\n\n",B);
+ }
+ for(i=0; i<row; i++){
+ fscanf(jobj,"%s",str);
+ fprintf(out,"%s",str);
+ for(k=0;k<col;k++){
+ fscanf(jobj,"%*lf",&placeholder);
+ }
+ for(j=0; j<9;j++){
+ fprintf(out,"\t%.12lf",storage[i][j]);
+ }
+ fprintf(out,"\n");
+ }
+
+ fprintf(out,"\n \n");
+
+ fclose(jobj);
+ fclose(out);
+
+ return 0;
+}
+
+void testp(double *permuted_ttests,int *B,double *permuted,
+ double *Imatrix,int *nc,int *nr,int *g,double *Tinitial,double
+ *ps) {
+
+ double Tvalues[*nr];
+ int a, b, n, i, j,k=0;
+
+ a = *B;
+ b = *nr;
+ n = a*b;
+
+ double counter[b];
+
+ for(j=0;j<b;j++){
+ counter[j]=0;
+ }
+
+ for (j=1; j<=*B; j++){
+ permute_matrix(Imatrix,nc,nr,permuted,g,Tvalues,Tinitial,counter);
+ // for(i=0;i<*nr;i++){
+ // permuted_ttests[k]=fabs(Tvalues[i]);
+ // k++;
+ }
+
+
+ for(j=0;j<*nr;j++){
+ ps[j]=((counter[j]+1)/(double)(a+1));
+ }
+}
+
+void permute_matrix(double *Imatrix,int *nc,int *nr,double *permuted,
+ int *g,double *trial_ts,double *Tinitial,double *counter1){
+
+ int i=0,j=0,n=0,a=0,b=0,f=0,c=0,k=0;
+
+ a = *nr; // number of rows
+ b = *nc;
+ n = a*b;
+
+ int y[b];
+
+ for (i=1; i<=*nc; i++){
+ y[i-1] = i;
+ }
+
+ permute_array(y, b);
+
+ for (i=0; i<*nc; i++){
+ f = y[i]; //column number
+ c=1;
+ c*=(f-1);
+ c*=a;
+ if (f == 1){
+ c = 0;
+ } // starting value position in the Imatrix
+ for(j=1; j<=*nr; j++){
+ permuted[k] = Imatrix[c];
+ c++;
+ k++;
+ }
+ }
+
+ calc_twosample_ts(permuted,g,nc,nr,trial_ts,Tinitial,counter1);
+}
+
+void permute_array(int *array, int n) {
+ static int seeded = 0;
+ int i;
+
+ if (! seeded) {
+ seeded = 1;
+ srandom(time(NULL));
+ }
+
+ for (i = 0; i < n; i++) {
+ int selection = random() % (n - i);
+ int tmp = array[i + selection];
+ array[i + selection] = array[i];
+ array[i] = tmp;
+ }
+}
+
+void calc_twosample_ts(double *Pmatrix,int *g,int *nc,int *nr,
+ double *Ts,double *Tinitial,double *counter) {
+ int i,a;
+ a = *nr;
+ a*=4;
+
+ double C1[*nr][3], C2[*nr][3], storage[a],tool[a];
+ double nrows,ncols,gvalue, xbardiff=0, denom=0;
+
+ nrows = (double) *nr;
+ ncols = (double) *nc;
+ gvalue= (double) *g;
+
+ meanvar(Pmatrix,g,nr,nc,storage);
+ for(i=0;i<=a-1;i++){
+ tool[i]=storage[i];
+ }
+ for (i=0; i<*nr;i++){
+ C1[i][0]=tool[i];
+ C1[i][1]=tool[i+*nr+*nr];
+ C1[i][2]=C1[i][1]/(gvalue-1);
+
+ C2[i][0]=tool[i+*nr];
+ C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
+ C2[i][2]=C2[i][1]/(ncols-gvalue+1);
+ }
+
+ for (i=0; i<*nr; i++){
+ xbardiff = C1[i][0]-C2[i][0];
+ denom = sqrt(C1[i][2]+C2[i][2]);
+ Ts[i]=fabs(xbardiff/denom);
+ if (fabs(Ts[i])>(fabs(Tinitial[i])+.0000000000001)){ //13th place
+ counter[i]++;
+ }
+ }
+}
+
+void meanvar(double *pmatrix,int *g,int *nr,int *nc,double *store){
+ double temp[*nr], temp2[*nr],var[*nr],var2[*nr],a,b;
+
+ int i,m,k,l,n;
+
+ a = (double) *g-1;
+ b = (double) (*nc-a);
+
+ for (i = 0; i<*nr; i++){
+ temp[i]=0;
+ temp2[i]=0;
+ var[i]=0;
+ var2[i]=0;
+ }
+
+ k = *nr; // number of rows
+ l = *nc;
+ n = k*l;
+
+ m=0;
+ m=*g-1;
+ k=*nr;
+ m*=k; // m = g * nr now
+ for (i=0;i<m;i++){
+ temp[i%k]=temp[i%k]+pmatrix[i];
+ }
+ for (i=0;i<n;i++){
+ temp2[i%k]=temp2[i%k]+pmatrix[i];
+ }
+ for (i=0;i<*nr;i++){
+ temp2[i]=temp2[i]-temp[i];
+ }
+ for (i=0;i<=*nr-1;i++){
+ store[i]=temp[i]/a;
+ store[i+*nr]=temp2[i]/b;
+ }
+
+ // That completes the mean calculations.
+
+ for (i=0;i<m;i++){
+ var[i%k]=var[i%k]+pow((pmatrix[i]-store[i%k]),2);
+ }
+ for (i=m;i<n;i++){
+ var2[i%k]=var2[i%k]+pow((pmatrix[i]-store[(i%k)+*nr]),2);
+ }
+
+ for (i=0;i<=*nr-1;i++){
+ store[i+2*k]=var[i]/(a-1);
+ store[i+3*k]=var2[i]/(b-1);
+ }
+ // That completes var calculations.
+}
+
+void start(double *Imatrix,int *g,int *nr,int *nc,double *initial,
+ double storage[][9]){
+ int i, a = *nr;
+ a*=4;
+
+ double store[a], tool[a], C1[*nr][3], C2[*nr][3];
+ double nrows,ncols,gvalue, xbardiff=0, denom=0;
+
+ nrows = (double) *nr;
+ ncols = (double) *nc;
+ gvalue= (double) *g;
+
+ meanvar(Imatrix,g,nr,nc,store);
+
+ for(i=0;i<=a-1;i++){
+ tool[i]=store[i];
+ }
+ for (i=0; i<*nr;i++){
+ C1[i][0]=tool[i]; //mean group 1
+ storage[i][0]=C1[i][0];
+ C1[i][1]=tool[i+*nr+*nr]; // var group 1
+ storage[i][1]=C1[i][1];
+ C1[i][2]=C1[i][1]/(gvalue-1);
+ storage[i][2]=sqrt(C1[i][2]);
+
+ C2[i][0]=tool[i+*nr]; // mean group 2
+ storage[i][4]=C2[i][0];
+ C2[i][1]=tool[i+*nr+*nr+*nr]; // var group 2
+ storage[i][5]=C2[i][1];
+ C2[i][2]=C2[i][1]/(ncols-gvalue+1);
+ storage[i][6]=sqrt(C2[i][2]);
+ }
+ for (i=0; i<*nr; i++){
+ xbardiff = C1[i][0]-C2[i][0];
+ denom = sqrt(C1[i][2]+C2[i][2]);
+ initial[i]=fabs(xbardiff/denom);
+ }
+}
projects / mothur.git / commitdiff