sharedthetan.cpp

   1 /*
   2  *  sharedthetan.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/8/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedthetan.h"
  11
  12 /***********************************************************************/
  13 EstOutput ThetaN::getValues(vector<SharedRAbundVector*> shared) {
  14         try {
  15                 data.resize(1,0);
  16
  17                 double Atotal, Btotal, tempA, tempB;
  18                 Atotal = 0; Btotal = 0;
  19                 double numerator, denominator, thetaN, sumSharedA, sumSharedB, a, b, d;
  20                 numerator = 0.0; denominator = 0.0; thetaN = 0.0; sumSharedA = 0.0; sumSharedB = 0.0; a = 0.0; b = 0.0; d = 0.0;
  21
  22                 //get the total values we need to calculate the theta denominator sums
  23                 for (int i = 0; i < shared[0]->getNumBins(); i++) {
  24                         //store in temps to avoid multiple repetitive function calls
  25                         Atotal += shared[0]->getAbundance(i);
  26                         Btotal += shared[1]->getAbundance(i);
  27                 }
  28
  29                 //calculate the theta denominator sums
  30                 for (int j = 0; j < shared[0]->getNumBins(); j++) {
  31                         //store in temps to avoid multiple repetitive function calls
  32                         tempA = shared[0]->getAbundance(j);
  33                         tempB = shared[1]->getAbundance(j);
  34
  35                         //they are shared
  36                         if ((tempA != 0) && (tempB != 0)) {
  37                                 if (Atotal != 0)        { sumSharedA = (tempA / (float)Atotal); }
  38                                 if (Btotal != 0)        { sumSharedB = (tempB / (float)Btotal); }
  39
  40                                 a += sumSharedA;
  41                                 b += sumSharedB;
  42                         }
  43                 }
  44
  45                 thetaN = (a * b) / (a + b - (a * b));
  46
  47                 if (isnan(thetaN) || isinf(thetaN)) { thetaN = 0; }
  48
  49                 data[0] = 1.0 - thetaN;
  50
  51                 return data;
  52         }
  53         catch(exception& e) {
  54                 m->errorOut(e, "ThetaN", "getValues");
  55                 exit(1);
  56         }
  57 }
  58
  59 /***********************************************************************/