wilcox.cpp

   1 //
   2 //  wilcox.cpp
   3 //  Mothur
   4 //
   5 //  Created by SarahsWork on 8/6/13.
   6 //  Copyright (c) 2013 Schloss Lab. All rights reserved.
   7 //
   8 //  Modified to avoid using Rmath.h
   9
  10 #include "mothur.h"
  11
  12 /*
  13  Mathlib : A C Library of Special Functions
  14  Copyright (C) 1999-2012  The R Core Team
  15
  16  This program is free software; you can redistribute it and/or modify
  17  it under the terms of the GNU General Public License as published by
  18  the Free Software Foundation; either version 2 of the License, or (at
  19  your option) any later version.
  20
  21  This program is distributed in the hope that it will be useful, but
  22  WITHOUT ANY WARRANTY; without even the implied warranty of
  23  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.   See the GNU
  24  General Public License for more details.
  25
  26  You should have received a copy of the GNU General Public License
  27  along with this program; if not, a copy is available at
  28  http://www.r-project.org/Licenses/
  29
  30  SYNOPSIS
  31
  32  #include <Rmath.h>
  33  double dwilcox(double x, double m, double n, int give_log)
  34  double pwilcox(double x, double m, double n, int lower_tail, int log_p)
  35  double qwilcox(double x, double m, double n, int lower_tail, int log_p);
  36  double rwilcox(double m, double n)
  37
  38  DESCRIPTION
  39
  40  dwilcox        The density of the Wilcoxon distribution.
  41  pwilcox        The distribution function of the Wilcoxon distribution.
  42  qwilcox        The quantile function of the Wilcoxon distribution.
  43  rwilcox        Random variates from the Wilcoxon distribution.
  44
  45  */
  46
  47 /*
  48  Note: the checks here for R_CheckInterrupt also do stack checking.
  49
  50  calloc/free are remapped for use in R, so allocation checks are done there.
  51  freeing is completed by an on.exit action in the R wrappers.
  52  */
  53 #include "wilcox.h"
  54
  55 #define WILCOX_MAX 50
  56 #define give_log log_p
  57 #define R_D__0  (log_p ? 1 : 0.)                /* 0 */
  58 #define R_D__1  (log_p ? 0. : 1.)                       /* 1 */
  59 #define R_DT_0  (lower_tail ? R_D__0 : R_D__1)          /* 0 */
  60 #define R_DT_1  (lower_tail ? R_D__1 : R_D__0)          /* 1 */
  61 #define R_D_val(x)      (log_p  ? log(x) : (x))         /*  x  in pF(x,..) */
  62 #define R_D_Clog(p)     (log_p  ? log1p(-(p)) : (0.5 - (p) + 0.5)) /* [log](1-p) */
  63 #define R_DT_val(x)     (lower_tail ? R_D_val(x)  : R_D_Clog(x))
  64
  65 /*********************************************************************************************************************************/
  66 //Numerical Recipes pg. 219
  67 double PWilcox::gammln(const double xx) {
  68     try {
  69         int j;
  70         double x,y,tmp,ser;
  71         static const double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,
  72             -1.231739572450155,0.120858003e-2,-0.536382e-5};
  73
  74         y=x=xx;
  75         tmp=x+5.5;
  76         tmp -= (x+0.5)*log(tmp);
  77         ser=1.0;
  78         for (j=0;j<6;j++) {
  79             ser += cof[j]/++y;
  80         }
  81         return -tmp+log(2.5066282746310005*ser/x);
  82     }
  83     catch(exception& e) {
  84         mout->errorOut(e, "PWilcox", "gammln");
  85         exit(1);
  86     }
  87 }
  88 /*********************************************************************************************************************************/
  89 double PWilcox::choose(double n, double k){
  90     try {
  91         n = floor(n + 0.5);
  92         k = floor(k + 0.5);
  93
  94         double lchoose = gammln(n + 1.0) - gammln(k + 1.0) - gammln(n - k + 1.0);
  95
  96         return (floor(exp(lchoose) + 0.5));
  97     }
  98     catch(exception& e) {
  99         mout->errorOut(e, "PWilcox", "choose");
 100         exit(1);
 101     }
 102 }
 103 /*********************************************************************************************************************************/
 104
 105 /* This counts the number of choices with statistic = k */
 106 double PWilcox::cwilcox(int k, int m, int n, double*** w) {
 107     try {
 108         int c, u, i, j, l;
 109
 110         if (mout->control_pressed) { return 0; }
 111
 112         u = m * n;
 113         if (k < 0 || k > u)
 114             return(0);
 115         c = (int)(u / 2);
 116         if (k > c)
 117             k = u - k; /* hence  k <= floor(u / 2) */
 118         if (m < n) {
 119             i = m; j = n;
 120         } else {
 121             i = n; j = m;
 122         } /* hence  i <= j */
 123
 124         if (j == 0) /* and hence i == 0 */
 125             return (k == 0);
 126
 127
 128         /* We can simplify things if k is small.  Consider the Mann-Whitney
 129          definition, and sort y.  Then if the statistic is k, no more
 130          than k of the y's can be <= any x[i], and since they are sorted
 131          these can only be in the first k.  So the count is the same as
 132          if there were just k y's.
 133          */
 134         if (j > 0 && k < j) return cwilcox(k, i, k, w);
 135
 136         if (w[i][j] == 0) {
 137             w[i][j] = (double *) calloc((size_t) c + 1, sizeof(double));
 138
 139             for (l = 0; l <= c; l++)
 140                 w[i][j][l] = -1;
 141         }
 142         if (w[i][j][k] < 0) {
 143             if (j == 0) /* and hence i == 0 */
 144                 w[i][j][k] = (k == 0);
 145             else
 146                 w[i][j][k] = cwilcox(k - j, i - 1, j, w) + cwilcox(k, i, j - 1, w);
 147
 148         }
 149         return(w[i][j][k]);
 150     }
 151     catch(exception& e) {
 152         mout->errorOut(e, "PWilcox", "cwilcox");
 153         exit(1);
 154     }
 155 }
 156 /*********************************************************************************************************************************/
 157
 158 /* args have the same meaning as R function pwilcox */
 159 double PWilcox::pwilcox(double q, double m, double n, bool lower_tail){
 160     try {
 161         int i;
 162         double c, p;
 163         bool log_p = false;
 164         double*** w;
 165
 166
 167         if (isnan(m) || isnan(n))
 168         {return 0;}
 169         m = floor(m + 0.5);
 170         n = floor(n + 0.5);
 171         if (m <= 0 || n <= 0) {return 0;}
 172
 173         q = floor(q + 1e-7);
 174
 175         if (q < 0.0)
 176             return(R_DT_0);
 177         if (q >= m * n)
 178             return(R_DT_1);
 179
 180         int mm = (int) m, nn = (int) n;
 181
 182         if (mout->control_pressed) { return 0; }
 183         //w_init_maybe(mm, nn);
 184         /********************************************/
 185         int thisi;
 186         if (mm > nn) { thisi = nn; nn = mm; mm = thisi; }
 187
 188         mm = max(mm, 50);
 189         nn = max(nn, 50);
 190         w = (double ***) calloc((size_t) mm + 1, sizeof(double **));
 191
 192         for (thisi = 0; thisi <= mm; thisi++) {
 193             w[thisi] = (double **) calloc((size_t) nn + 1, sizeof(double *));
 194         }
 195         allocated_m = m; allocated_n = n;
 196         /********************************************/
 197
 198         c = choose(m + n, n);
 199         p = 0;
 200         /* Use summation of probs over the shorter range */
 201         if (q <= (m * n / 2)) {
 202             for (i = 0; i <= q; i++)
 203                 p += cwilcox(i, m, n, w) / c;
 204         }
 205         else {
 206             q = m * n - q;
 207             for (i = 0; i < q; i++) {
 208                 p += cwilcox(i, m, n, w) / c;
 209             }
 210             lower_tail = !lower_tail; /* p = 1 - p; */
 211         }
 212
 213         //free w
 214         /********************************************/
 215         for (int i = allocated_m; i >= 0; i--) {
 216             for (int j = allocated_n; j >= 0; j--) {
 217                 if (w[i][j] != 0)
 218                     free((void *) w[i][j]);
 219             }
 220             free((void *) w[i]);
 221         }
 222         free((void *) w);
 223         w = 0; allocated_m = allocated_n = 0;
 224         /********************************************/
 225
 226         return(R_DT_val(p));
 227     }
 228     catch(exception& e) {
 229         mout->errorOut(e, "PWilcox", "pwilcox");
 230         exit(1);
 231     }
 232 } /* pwilcox */
 233