X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=linearalgebra.cpp;h=bacf2c1c07bb420430d07416b6ca7354c650a79b;hp=fdd95ec9602cb3b44cf3e1aa54934fc3045982c8;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=3a5dd9e428ab93a6dcdce7912e8ebb977be0b893

diff --git a/linearalgebra.cpp b/linearalgebra.cpp
index fdd95ec..bacf2c1 100644
--- a/linearalgebra.cpp
+++ b/linearalgebra.cpp
@@ -10,6 +10,8 @@
 #include "linearalgebra.h"
 #include "wilcox.h"
 
+#define PI 3.1415926535897932384626433832795
+
 // This class references functions used from "Numerical Recipes in C++" //
 
 /*********************************************************************************************************************************/
@@ -112,7 +114,7 @@ double LinearAlgebra::gammp(const double a, const double x) {
 }
 /*********************************************************************************************************************************/
 //Numerical Recipes pg. 223
-double LinearAlgebra::gammq(const double a, const double x) {
+/*double LinearAlgebra::gammq(const double a, const double x) {
     try {
         double gamser,gammcf,gln;
         
@@ -128,11 +130,11 @@ double LinearAlgebra::gammq(const double a, const double x) {
         return 0;
     }
 	catch(exception& e) {
-		m->errorOut(e, "LinearAlgebra", "gammp");
+		m->errorOut(e, "LinearAlgebra", "gammq");
 		exit(1);
 	}
 }
-/*********************************************************************************************************************************/
+*********************************************************************************************************************************/
 //Numerical Recipes pg. 224
 double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double& gln){
     try {
@@ -159,7 +161,7 @@ double LinearAlgebra::gcf(double& gammcf, const double a, const double x, double
             h *= del;
             if (fabs(del-1.0) <= EPS) break;
         }
-        if (i > ITMAX)  { m->mothurOut("[ERROR]: a too large, ITMAX too small in gcf\n"); m->control_pressed = true; }
+        if (i > ITMAX)  { m->mothurOut("[ERROR]: " + toString(a) + " too large, ITMAX=100 too small in gcf\n"); m->control_pressed = true; }
         gammcf=exp(-x+a*log(x)-gln)*h;
         
         return 0.0;
@@ -1300,7 +1302,7 @@ double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pV
         double lastTerm = 3 * (values.size()+1);
         
         H = firstTerm * middleTerm - lastTerm;
-        
+       
         //adjust for ties
         if (TIES.size() != 0) {
             double sum = 0.0;
@@ -1309,6 +1311,8 @@ double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pV
             H /= result;
         }
         
+        if (isnan(H) || isinf(H)) { H = 0; }
+        
         //Numerical Recipes pg221
         pValue = 1.0 - (gammp(((treatments.size()-1)/(double)2.0), H/2.0));
         
@@ -1320,24 +1324,14 @@ double LinearAlgebra::calcKruskalWallis(vector<spearmanRank>& values, double& pV
 	}
 }
 /*********************************************************************************************************************************/
-//python random.normalvariate - thanks http://madssj.com/blog/2008/05/07/porting-normalvariate-from-python-to-c/
-double LinearAlgebra::normalvariate(double mu, double sigma) {
+double LinearAlgebra::normalvariate(double mean, double standardDeviation) {
     try {
-        double NV_MAGICCONST = 1.7155277699214135; /* (4 * exp(-0.5) / sqrt(2.0)); */
-        unsigned long int MAX_RANDOM = 2147483647; /* (2 ** 31) - 1; */
-        
-        double u1, u2, z, zz;
-        for (;;) {
-            if (m->control_pressed) { break; }
-            u1 = ((float)random()) / MAX_RANDOM;
-            u2 = 1.0 - (((float)random()) / MAX_RANDOM);
-            z = NV_MAGICCONST * (u1 - 0.5) / u2;
-            zz = z * z / 4.0;
-            if (zz <= -(log(u2))) {
-                break;
-            }
-        }
-        return mu + z * sigma;
+        double u1 = ((double)(rand()) + 1.0 )/( (double)(RAND_MAX) + 1.0);
+        double u2 = ((double)(rand()) + 1.0 )/( (double)(RAND_MAX) + 1.0);
+        //double r = sqrt( -2.0*log(u1) );
+        //double theta = 2.0*PI*u2;
+        //cout << cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1)) << endl;
+        return cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1));
     }
 	catch(exception& e) {
 		m->errorOut(e, "LinearAlgebra", "normalvariate");
@@ -1367,6 +1361,7 @@ double LinearAlgebra::pnorm(double x){
         double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
         
         return 0.5*(1.0 + sign*y);
+        
     }
 	catch(exception& e) {
 		m->errorOut(e, "LinearAlgebra", "pnorm");