sharedbdiversity.cpp

   1 /*
   2  *  bdiversity.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 3/13/09.
   6  *  Copyright 2009 __MyCompanyName__. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedbdiversity.h"
  11 #include "calculator.h"
  12
  13 /***********************************************************************/
  14
  15 double SharedBDiversity::whitt(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
  16         double nz1 = (double)shared1->numNZ();
  17         double nz2 = (double)shared2->numNZ();
  18         double sum = nz1;
  19         for(int i = 1; i < shared1->size(); i++)
  20         {
  21                 if(shared1->get(i).abundance == 0 && shared2->get(i).abundance != 0)
  22                         sum++;
  23                         }
  24         return 2*sum/(nz1+nz2)-1;
  25 }
  26
  27 /***********************************************************************/
  28
  29 double SharedBDiversity::ms(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
  30         double a = 0;
  31         double b = 0;
  32         double c = 0;
  33         for(int i = 1; i < shared1->size(); i++)
  34         {
  35                 int abund1 = shared1->get(i).abundance;
  36                 int abund2 = shared2->get(i).abundance;
  37
  38                 if(abund1 > 0 && abund2 > 0)
  39                         a++;
  40                 else if(abund1 > 0 && abund2 == 0)
  41                         b++;
  42                 else if(abund1 == 0 && abund2 > 0)
  43                         c++;
  44         }
  45         return (b+c)/(a+b+c);
  46 }
  47
  48 /***********************************************************************/
  49
  50 double SharedBDiversity::sor(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
  51         double sum = 0;
  52         double asum = 0;
  53         double bsum = 0;
  54         for(int i = 1; i < shared1->size(); i++)
  55         {
  56                 int abund1 = shared1->get(i).abundance;
  57                 int abund2 = shared2->get(i).abundance;
  58
  59                 asum += abund1;
  60                 bsum += abund2;
  61                 if(abund1 >= abund2)
  62                         sum += abund2;
  63                 else
  64                         sum += abund1;
  65         }
  66         return 2*sum/(asum+bsum);
  67 }
  68
  69 /***********************************************************************/
  70
  71 double SharedBDiversity::mor(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
  72         double multSum = 0;
  73         double powSum1 = 0;
  74         double powSum2 = 0;
  75         double ind1 = 0;
  76         double ind2 = 0;
  77         for(int i = 1; i < shared1->size(); i++)
  78         {
  79                 double abund1 = (double)shared1->get(i).abundance;
  80                 double abund2 = (double)shared2->get(i).abundance;
  81                 multSum += abund1*abund2;
  82                 powSum1 += pow(abund1, 2);
  83                 powSum2 += pow(abund2, 2);
  84                 ind1 += abund1;
  85                 ind2 += abund2;
  86         }
  87         return 2*multSum / ((powSum1/pow(ind1, 2) + powSum2/pow(ind2, 2)) * (ind1*ind2));
  88 }
  89 /***********************************************************************/
  90
  91 EstOutput SharedBDiversity::getValues(SharedRAbundVector* shared1, SharedRAbundVector* shared2){
  92         try {
  93                 data.resize(4,0);
  94
  95                 double whitt1 = whitt(shared1, shared2);
  96                 double jac1 = 1-ms(shared1, shared2);
  97                 double sor1 = sor(shared1, shared2);
  98                 double mor1 = mor(shared1, shared2);
  99                 /*cout << "Whittaker's Measure: " << whitt1 << "\n";
 100                 cout << "Marczewski-Steinhaus distance: " << jac1 << "\n";
 101                 cout << "Sorensen index: " << sor1 << "\n";
 102                 cout << "Morisita-Horn index: " << mor1 << "\n";*/
 103
 104                 data[0] = whitt1;
 105                 data[1] = jac1;
 106                 data[2] = sor1;
 107                 data[3] = mor1;
 108                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
 109                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 110                 if (isnan(data[2]) || isinf(data[0])) { data[2] = 0; }
 111                 if (isnan(data[3]) || isinf(data[1])) { data[3] = 0; }
 112
 113                 return data;
 114         }
 115
 116         catch(exception& e) {
 117                 cout << "Standard Error: " << e.what() << " has occurred in the NPShannon class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 118                 exit(1);
 119         }
 120         catch(...) {
 121                 cout << "An unknown error has occurred in the NPShannon class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 122                 exit(1);
 123         }
 124 }
 125
 126 /***********************************************************************/