From 4bc720215844701594bcc61067a2cf79a232ff96 Mon Sep 17 00:00:00 2001 From: Pat Schloss Date: Thu, 9 Jan 2014 10:34:03 -0500 Subject: [PATCH] converted jsd calculation to be based on relative abundance --- sharedjsd.cpp | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/sharedjsd.cpp b/sharedjsd.cpp index 46cb977..914c875 100644 --- a/sharedjsd.cpp +++ b/sharedjsd.cpp @@ -19,21 +19,35 @@ EstOutput JSD::getValues(vector shared) { double KLD1 = 0.0; double KLD2 = 0.0; + vector countA = shared[0]->getAbundances(); + vector countB = shared[1]->getAbundances(); + double totalA = 0; + double totalB = 0; + for (int i = 0; i < shared[0]->getNumBins(); i++) { - //store in temps to avoid multiple repetitive function calls - double tempA = shared[0]->getAbundance(i); - double tempB = shared[1]->getAbundance(i); + totalA += countA[i]; + totalB += countB[i]; + } + + for (int i = 0; i < shared[0]->getNumBins(); i++) { + double tempA = countA[i] / totalA; + double tempB = countB[i] / totalB; + tempA = countA[i] / totalA; + tempB = countB[i] / totalB; + if (tempA == 0) { tempA = 0.000001; } if (tempB == 0) { tempB = 0.000001; } - + double denom = (tempA+tempB)/(double)2.0; - + if (tempA != 0) { KLD1 += tempA * log(tempA/denom); } //KLD(x,m) if (tempB != 0) { KLD2 += tempB * log(tempB/denom); } //KLD(y,m) + } - - data[0] = sqrt((0.5*KLD1) + (0.5*KLD2)); + + + data[0] = ((0.5*KLD1) + (0.5*KLD2)); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } -- 2.39.2