From: westcott Date: Thu, 16 Sep 2010 10:34:29 +0000 (+0000) Subject: initial add of metastats X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=2bf3df7736ef2a17286d99394e211f51751d6829 initial add of metastats --- diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 875a3c5..0594ef6 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -486,6 +486,9 @@ A7DF0AE2121EBB14004A03EA /* prng.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = prng.h; sourceTree = ""; }; A7E8338B115BBDAA00739EC4 /* parsesffcommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = parsesffcommand.cpp; sourceTree = ""; }; A7E8338C115BBDAA00739EC4 /* parsesffcommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = parsesffcommand.h; sourceTree = ""; }; + A7F6C8E1124229F900299875 /* fisher2.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = fisher2.c; sourceTree = ""; }; + A7F6C8E2124229F900299875 /* fisher2.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = fisher2.h; sourceTree = ""; }; + A7F6C8E3124229F900299875 /* metastats2.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = metastats2.c; sourceTree = ""; }; /* End PBXFileReference section */ /* Begin PBXGroup section */ @@ -537,6 +540,7 @@ A7DA2082113FECD400BF472F /* inputdata.h */, A7DA208B113FECD400BF472F /* libshuff.cpp */, A7DA208C113FECD400BF472F /* libshuff.h */, + A7F6C8DC1242299300299875 /* metastatsource */, A7DA209E113FECD400BF472F /* mothur.cpp */, A7DA209F113FECD400BF472F /* mothur.h */, A7DA20A0113FECD400BF472F /* mothurout.cpp */, @@ -1031,6 +1035,16 @@ name = clearcutsource; sourceTree = ""; }; + A7F6C8DC1242299300299875 /* metastatsource */ = { + isa = PBXGroup; + children = ( + A7F6C8E2124229F900299875 /* fisher2.h */, + A7F6C8E1124229F900299875 /* fisher2.c */, + A7F6C8E3124229F900299875 /* metastats2.c */, + ); + name = metastatsource; + sourceTree = ""; + }; /* End PBXGroup section */ /* Begin PBXLegacyTarget section */ diff --git a/clearcutcommand.cpp b/clearcutcommand.cpp index 1a7ab1f..c76de38 100644 --- a/clearcutcommand.cpp +++ b/clearcutcommand.cpp @@ -180,7 +180,9 @@ int ClearcutCommand::execute() { vector cPara; - char* tempClearcut = new char[8]; strcpy(tempClearcut, "clearcut"); cPara.push_back(tempClearcut); + char* tempClearcut = new char[8]; + strcpy(tempClearcut, "clearcut"); + cPara.push_back(tempClearcut); //you gave us a distance matrix if (phylipfile != "") { char* temp = new char[10]; strcpy(temp, "--distance"); cPara.push_back(temp); } diff --git a/fisher2.c b/fisher2.c new file mode 100644 index 0000000..9e62302 --- /dev/null +++ b/fisher2.c @@ -0,0 +1,2168 @@ +#include +#include +#include +//#include "ctest.h" + + +#include +#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 */ + diff --git a/fisher2.h b/fisher2.h new file mode 100644 index 0000000..cba5966 --- /dev/null +++ b/fisher2.h @@ -0,0 +1,9 @@ +//#ifndef GUARD_fisher2 +//#define GUARD_fisher2 + +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); + +//#endif GUARD_fisher2 + diff --git a/heatmapcommand.cpp b/heatmapcommand.cpp index fb8e021..bb82702 100644 --- a/heatmapcommand.cpp +++ b/heatmapcommand.cpp @@ -184,8 +184,9 @@ int HeatMapCommand::execute(){ if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) { string saveLabel = lookup[0]->getLabel(); - + for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } + lookup = input->getSharedRAbundVectors(lastLabel); m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); diff --git a/inputdata.cpp b/inputdata.cpp index 5d71450..dc99a68 100644 --- a/inputdata.cpp +++ b/inputdata.cpp @@ -16,6 +16,7 @@ InputData::InputData(string fName, string f) : format(f){ m = MothurOut::getInstance(); + globaldata = GlobalData::getInstance(); m->openInputFile(fName, fileHandle); filename = fName; @@ -37,6 +38,7 @@ InputData::~InputData(){ InputData::InputData(string fName, string orderFileName, string f) : format(f){ try { m = MothurOut::getInstance(); + globaldata = GlobalData::getInstance(); ifstream ofHandle; m->openInputFile(orderFileName, ofHandle); string name; @@ -453,7 +455,8 @@ vector InputData::getSharedRAbundVectors(string label){ string thisLabel; m->openInputFile(filename, in); - + globaldata->saveNextLabel = ""; + if(in){ if (format == "sharedfile") { while (in.eof() != true) { @@ -461,6 +464,7 @@ vector InputData::getSharedRAbundVectors(string label){ SharedRAbundVector* SharedRAbund = new SharedRAbundVector(in); if (SharedRAbund != NULL) { thisLabel = SharedRAbund->getLabel(); + //if you are at the last label if (thisLabel == label) { in.close(); return SharedRAbund->getSharedRAbundVectors(); } else { @@ -536,6 +540,7 @@ vector InputData::getSharedRAbundFloatVectors(string l string thisLabel; m->openInputFile(filename, in); + globaldata->saveNextLabel = ""; if(in){ if (format == "sharedfile") { diff --git a/metastats2.c b/metastats2.c new file mode 100644 index 0000000..c4d1f46 --- /dev/null +++ b/metastats2.c @@ -0,0 +1,672 @@ +#include +#include +#include +#include +#include +#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;i2){ + for(i=1;i 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.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.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(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;icontrol_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; } for (int i = 1; i < processors; i++) { @@ -253,10 +253,10 @@ int ScreenSeqsCommand::execute(){ numSeqsPerProcessor = numFastaSeqs / processors; int startIndex = pid * numSeqsPerProcessor; if(pid == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor; } - cout << pid << '\t' << numSeqsPerProcessor << '\t' << startIndex << endl; + //cout << pid << '\t' << numSeqsPerProcessor << '\t' << startIndex << endl; //align your part driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPIGood, outMPIBadAccnos, MPIPos, badSeqNames); -cout << pid << " done" << endl; +//cout << pid << " done" << endl; if (m->control_pressed) { MPI_File_close(&inMPI); MPI_File_close(&outMPIGood); MPI_File_close(&outMPIBadAccnos); return 0; } //send bad list diff --git a/seqsummarycommand.cpp b/seqsummarycommand.cpp index 5dbcdec..abef126 100644 --- a/seqsummarycommand.cpp +++ b/seqsummarycommand.cpp @@ -347,10 +347,10 @@ int SeqSummaryCommand::driverCreateSummary(vector& startPosition, vector= filePos->end)) { break; } //report progress - if((count) % 100 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } + //if((count) % 100 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } } //report progress - if((count) % 100 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } + //if((count) % 100 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); } in.close();