From 0adea182fc7adb2d83492ac2c065fb30e0054e41 Mon Sep 17 00:00:00 2001 From: westcott Date: Wed, 26 Jan 2011 18:45:51 +0000 Subject: [PATCH] fixed mempearson and structpearson --- mempearson.cpp | 16 +++++++--------- structpearson.cpp | 22 +++++++--------------- 2 files changed, 14 insertions(+), 24 deletions(-) diff --git a/mempearson.cpp b/mempearson.cpp index b39afd8..934bd01 100644 --- a/mempearson.cpp +++ b/mempearson.cpp @@ -24,10 +24,11 @@ EstOutput MemPearson::getValues(vector shared) { if (shared[1]->getAbundance(i) != 0) { nonZeroB++; } } - double numTerm1 = 0.0; - double numTerm2 = 0.0; + double numTerm = 0.0; double denomTerm1 = 0.0; double denomTerm2 = 0.0; + double averageA = nonZeroA / (float) numOTUS; + double averageB = nonZeroB / (float) numOTUS; for (int i = 0; i < shared[0]->getNumBins(); i++) { int Aij = shared[0]->getAbundance(i); @@ -36,21 +37,18 @@ EstOutput MemPearson::getValues(vector shared) { if (Aij > 0) { Aij = 1; } if (Bij > 0) { Bij = 1; } - numTerm1 += (Aij - (nonZeroA / (float) numOTUS)); - numTerm2 += (Bij - (nonZeroB / (float) numOTUS)); - - denomTerm1 += ((Aij - (nonZeroA / (float) numOTUS)) * (Aij - (nonZeroA / (float) numOTUS))); - denomTerm2 += ((Bij - (nonZeroB / (float) numOTUS)) * (Bij - (nonZeroB / (float) numOTUS))); + numTerm += ((Aij - averageA) * (Bij - averageB)); + denomTerm1 += ((Aij - averageA) * (Aij - averageA)); + denomTerm2 += ((Bij - averageB) * (Bij - averageB)); } denomTerm1 = sqrt(denomTerm1); denomTerm2 = sqrt(denomTerm2); double denom = denomTerm1 * denomTerm2; - double numerator = numTerm1 * numTerm2; if (denom != 0) { - data[0] = 1.0 - (numerator / denom); + data[0] = 1.0 - (numTerm / denom); }else { data[0] = 1.0; } diff --git a/structpearson.cpp b/structpearson.cpp index 9d030ef..bc2ac76 100644 --- a/structpearson.cpp +++ b/structpearson.cpp @@ -14,17 +14,11 @@ EstOutput StructPearson::getValues(vector shared) { try { data.resize(1,0); - double sumA = 0.0; - double sumB = 0.0; int numOTUS = shared[0]->getNumBins(); + double averageA = shared[0]->getNumSeqs() / (float) numOTUS; + double averageB = shared[1]->getNumSeqs() / (float) numOTUS; - for (int i = 0; i < shared[0]->getNumBins(); i++) { - sumA += shared[0]->getAbundance(i); - sumB += shared[1]->getAbundance(i); - } - - double numTerm1 = 0.0; - double numTerm2 = 0.0; + double numTerm = 0.0; double denomTerm1 = 0.0; double denomTerm2 = 0.0; @@ -32,20 +26,18 @@ EstOutput StructPearson::getValues(vector shared) { int Aij = shared[0]->getAbundance(i); int Bij = shared[1]->getAbundance(i); - numTerm1 += (Aij - (sumA / (float) numOTUS)); - numTerm2 += (Bij - (sumB / (float) numOTUS)); - denomTerm1 += ((Aij - (sumA / (float) numOTUS)) * (Aij - (sumA / (float) numOTUS))); - denomTerm2 += ((Bij - (sumB / (float) numOTUS)) * (Bij - (sumB / (float) numOTUS))); + numTerm += ((Aij - averageA) * (Bij - averageB)); + denomTerm1 += ((Aij - averageA) * (Aij - averageA)); + denomTerm2 += ((Bij - averageB) * (Bij - averageB)); } denomTerm1 = sqrt(denomTerm1); denomTerm2 = sqrt(denomTerm2); double denom = denomTerm1 * denomTerm2; - double numerator = numTerm1 * numTerm2; - data[0] = 1.0 - (numerator / denom); + data[0] = 1.0 - (numTerm / denom); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } -- 2.39.2