]> git.donarmstrong.com Git - mothur.git/blob - memchi2.cpp
added odum, canberra, structchi2, structchord, structeuclidean, gower, hellinger...
[mothur.git] / memchi2.cpp
1 /*
2  *  memchi2.cpp
3  *  Mothur
4  *
5  *  Created by westcott on 12/17/10.
6  *  Copyright 2010 Schloss Lab. All rights reserved.
7  *
8  */
9
10 #include "memchi2.h"
11
12 /***********************************************************************/
13 EstOutput MemChi2::getValues(vector<SharedRAbundVector*> shared) {
14         try {
15                 data.resize(1,0);
16                 
17                 int nonZeroA = 0;
18                 int nonZeroB = 0;
19                 int totalOtus = shared[0]->getNumBins();
20                 int totalGroups = shared.size();
21                 
22                 //for each otu
23                 for (int i = 0; i < shared[0]->getNumBins(); i++) {
24                         if (shared[0]->getAbundance(i) != 0) { nonZeroA++; }
25                         if (shared[1]->getAbundance(i) != 0) { nonZeroB++; }
26                 }
27                 
28                 double totalTerm = 1 / (float) totalGroups;
29                 double sum = 0.0;
30                 for (int i = 0; i < shared[0]->getNumBins(); i++) {
31                         int A = shared[0]->getAbundance(i);
32                         int B = shared[1]->getAbundance(i);
33                         
34                         if (A > 0) { A = 1; }
35                         if (B > 0) { B = 1; }
36                         
37                         double Aterm = A / (float) nonZeroA;
38                         double Bterm = B / (float) nonZeroB;
39                         
40                         sum += (totalTerm * ((Aterm-Bterm)*(Aterm-Bterm)));
41                 }
42                 
43                 data[0] = sqrt((totalOtus * sum));
44                 
45                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
46                 
47                 return data;
48         }
49         catch(exception& e) {
50                 m->errorOut(e, "MemChi2", "getValues");
51                 exit(1);
52         }
53 }
54 /***********************************************************************/
55