5 * Created by westcott on 12/17/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "mempearson.h"
12 /***********************************************************************/
13 EstOutput MemPearson::getValues(vector<SharedRAbundVector*> shared) {
19 int numOTUS = shared[0]->getNumBins();
22 for (int i = 0; i < shared[0]->getNumBins(); i++) {
23 if (shared[0]->getAbundance(i) != 0) { nonZeroA++; }
24 if (shared[1]->getAbundance(i) != 0) { nonZeroB++; }
27 double numTerm1 = 0.0;
28 double numTerm2 = 0.0;
29 double denomTerm1 = 0.0;
30 double denomTerm2 = 0.0;
32 for (int i = 0; i < shared[0]->getNumBins(); i++) {
33 int Aij = shared[0]->getAbundance(i);
34 int Bij = shared[1]->getAbundance(i);
36 if (Aij > 0) { Aij = 1; }
37 if (Bij > 0) { Bij = 1; }
39 numTerm1 += (Aij - (nonZeroA / (float) numOTUS));
40 numTerm2 += (Bij - (nonZeroB / (float) numOTUS));
42 denomTerm1 += ((Aij - (nonZeroA / (float) numOTUS)) * (Aij - (nonZeroA / (float) numOTUS)));
43 denomTerm2 += ((Bij - (nonZeroB / (float) numOTUS)) * (Bij - (nonZeroB / (float) numOTUS)));
46 denomTerm1 = sqrt(denomTerm1);
47 denomTerm2 = sqrt(denomTerm2);
49 double denom = denomTerm1 * denomTerm2;
50 double numerator = numTerm1 * numTerm2;
53 data[0] = 1.0 - (numerator / denom);
58 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
63 m->errorOut(e, "MemPearson", "getValues");
67 /***********************************************************************/