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