calculator.cpp

   1 /*
   2  *  calculator.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 11/18/08.
   6  *  Copyright 2008 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "calculator.h"
  11
  12 /***********************************************************************/
  13 void VecCalc::printElements(vector<double> vec){
  14         try {
  15                 for(int i = 0; i < vec.size(); i++)
  16                         cout << vec.at(i) << " ";
  17                 cout << "\n\n";
  18         }
  19         catch(exception& e) {
  20                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  21                 exit(1);
  22         }
  23         catch(...) {
  24                 cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  25                 exit(1);
  26         }
  27 }
  28
  29 /***********************************************************************/
  30
  31 void VecCalc::printElements(vector<string> vec){
  32         try {
  33                 for(int i = 0; i < vec.size(); i++)
  34                         cout << vec.at(i) << " ";
  35                 cout << "\n\n";
  36         }
  37         catch(exception& e) {
  38                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  39                 exit(1);
  40         }
  41         catch(...) {
  42                 cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  43                 exit(1);
  44         }
  45 }
  46 /***********************************************************************/
  47 int VecCalc::findString(vector<string> vec, string str){
  48         try {
  49                 for(int i = 0; i < vec.size(); i++)
  50                         if(str.compare(vec.at(i)) == 0)
  51                                 return i;
  52                 return -1;
  53         }
  54         catch(exception& e) {
  55                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  56                 exit(1);
  57         }
  58         catch(...) {
  59                 cout << "An unknown error has occurred in the VecCalc class function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  60                 exit(1);
  61         }
  62 }
  63 /***********************************************************************/
  64 double VecCalc::mean(vector<double> vec){
  65         return sumElements(vec)/vec.size();
  66 }
  67 /***********************************************************************/
  68 double VecCalc::stError(vector<double> vec){
  69         try {
  70                 double sum = 0;
  71                 double m = mean(vec);
  72                 for(int i = 0; i < vec.size(); i++)
  73                         sum += pow(vec.at(i)-m,2);
  74                 return pow(sum/(vec.size()*(vec.size()-1)), .5);
  75         }
  76         catch(exception& e) {
  77                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  78                 exit(1);
  79         }
  80         catch(...) {
  81                 cout << "An unknown error has occurred in the VecCalc class function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  82                 exit(1);
  83         }
  84 }
  85 /***********************************************************************/
  86 double VecCalc::sumElements(vector<double> vec){
  87         try {
  88                 double sum = 0;
  89                 for(int i = 0; i < vec.size(); i++)
  90                         sum += vec.at(i);
  91                 return sum;
  92         }
  93         catch(exception& e) {
  94                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  95                 exit(1);
  96         }
  97         catch(...) {
  98                 cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  99                 exit(1);
 100         }
 101 }
 102 /***********************************************************************/
 103 double VecCalc::sumElements(vector<double> vec, int index){
 104         try {
 105                 double sum = 0;
 106                 for(int i = index; i < vec.size(); i++)
 107                         sum += vec.at(i);
 108                 return sum;
 109         }
 110         catch(exception& e) {
 111                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 112                 exit(1);
 113         }
 114         catch(...) {
 115                 cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 116                 exit(1);
 117         }
 118 }
 119 /***********************************************************************/
 120 double VecCalc::findMax(vector<double> vec){
 121         try {
 122                 double max = -1000000.0;
 123                 for(int i = 0; i < vec.size(); i++)
 124                         if(vec.at(i) > max)
 125                                 max = vec.at(i);
 126                 return max;
 127         }
 128         catch(exception& e) {
 129                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 130                 exit(1);
 131         }
 132         catch(...) {
 133                 cout << "An unknown error has occurred in the VecCalc class function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 134                 exit(1);
 135         }
 136 }
 137 /***********************************************************************/
 138 double VecCalc::numNZ(vector<double> vec){
 139         try {
 140                 double numNZ = 0;
 141                 for(int i = 0; i < vec.size(); i++)
 142                         if(vec.at(i) != 0)
 143                                 numNZ++;
 144                 return numNZ;
 145         }
 146         catch(exception& e) {
 147                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 148                 exit(1);
 149         }
 150         catch(...) {
 151                 cout << "An unknown error has occurred in the VecCalc class function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 152                 exit(1);
 153         }
 154 }
 155 /***********************************************************************/
 156 double VecCalc::numPos(vector<double> vec){
 157         try {
 158                 double numPos = 0;
 159                 for(int i = 0 ; i < vec.size(); i++)
 160                         if(vec.at(i) > 0)
 161                                 numPos++;
 162                 return numPos;
 163         }
 164         catch(exception& e) {
 165                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 166                 exit(1);
 167         }
 168         catch(...) {
 169                 cout << "An unknown error has occurred in the VecCalc class function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 170                 exit(1);
 171         }
 172 }
 173 /***********************************************************************/
 174 double VecCalc::findMaxDiff(vector<double> vec1, vector<double> vec2){
 175         try {
 176                 double maxDiff = -10000000000.0;
 177                 for(int i = 0; i < vec1.size(); i++)
 178                         if(fabs(vec1.at(i)-vec2.at(i)) > maxDiff)
 179                                 maxDiff = fabs(vec1.at(i)-vec2.at(i));
 180                 return maxDiff;
 181         }
 182         catch(exception& e) {
 183                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 184                 exit(1);
 185         }
 186         catch(...) {
 187                 cout << "An unknown error has occurred in the VecCalc class function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 188                 exit(1);
 189         }
 190 }
 191 /***********************************************************************/
 192 double VecCalc::findDStat(vector<double> vec1, vector<double> vec2, double num){
 193         try {
 194                 double mDiff = findMaxDiff(vec1, vec2);
 195                 return (mDiff+.5)/num;
 196         }
 197         catch(exception& e) {
 198                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 199                 exit(1);
 200         }
 201         catch(...) {
 202                 cout << "An unknown error has occurred in the VecCalc class function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 203                 exit(1);
 204         }
 205 }
 206 /***********************************************************************/
 207 vector<int> VecCalc::findQuartiles(vector<double> vec){
 208         try {
 209                 vector<int> quartiles;
 210                 double max = vec.at(vec.size()-1);
 211                 double r1 = max/4;
 212                 double r2 = max*3/4;
 213                 bool r1found = false;
 214                 bool r2found = false;
 215                 for(int i = 0; i < vec.size(); i++)
 216                 {
 217                         if(vec.at(i) > r1 && !r1found)
 218                         {
 219                                 quartiles.push_back(i);
 220                                 r1found = true;
 221                         }
 222                         if(vec.at(i) > r2 && !r2found)
 223                         {
 224                                 quartiles.push_back(i);
 225                                 i = vec.size();
 226                         }
 227                 }
 228                 return quartiles;
 229         }
 230         catch(exception& e) {
 231                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 232                 exit(1);
 233         }
 234         catch(...) {
 235                 cout << "An unknown error has occurred in the VecCalc class function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 236                 exit(1);
 237         }
 238 }
 239 /***********************************************************************/
 240 vector<double> VecCalc::add(vector<double> vec, double x){
 241         try {
 242                 vector<double> newVec;
 243                 for(int i = 0; i < vec.size(); i++)
 244                         newVec.push_back(vec.at(i)+x);
 245                 return newVec;
 246         }
 247         catch(exception& e) {
 248                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 249                 exit(1);
 250         }
 251         catch(...) {
 252                 cout << "An unknown error has occurred in the VecCalc class function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 253                 exit(1);
 254         }
 255 }
 256 /***********************************************************************/
 257 vector<double> VecCalc::multiply(vector<double> vec, double x){
 258         try {
 259                 vector<double> newVec;
 260                 for(int i = 0; i < vec.size(); i++)
 261                         newVec.push_back(vec.at(i)*x);
 262                 return newVec;
 263         }
 264         catch(exception& e) {
 265                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 266                 exit(1);
 267         }
 268         catch(...) {
 269                 cout << "An unknown error has occurred in the VecCalc class function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 270                 exit(1);
 271         }
 272 }
 273 /***********************************************************************/
 274 vector<double> VecCalc::power(vector<double> vec, double p){
 275         try {
 276                 vector<double> newVec;
 277                 for(int i = 0; i < vec.size(); i++)
 278                         newVec.push_back(pow(vec.at(i), p));
 279                 return newVec;
 280         }
 281         catch(exception& e) {
 282                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 283                 exit(1);
 284         }
 285         catch(...) {
 286                 cout << "An unknown error has occurred in the VecCalc class function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 287                 exit(1);
 288         }
 289 }
 290 /***********************************************************************/
 291 vector<double> VecCalc::addVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
 292 {       try {
 293                 vector<double> newVec;
 294                 for(int i = 0; i < vec1.size(); i++)
 295                         newVec.push_back(vec1.at(i)+vec2.at(i));
 296                 return newVec;
 297         }
 298         catch(exception& e) {
 299                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 300                 exit(1);
 301         }
 302         catch(...) {
 303                 cout << "An unknown error has occurred in the VecCalc class function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 304                 exit(1);
 305         }
 306 }
 307 /***********************************************************************/
 308 vector<double> VecCalc::multVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
 309 {       try {
 310                 vector<double> newVec;
 311                 for(int i = 0; i < vec1.size(); i++)
 312                         newVec.push_back(vec1.at(i)*vec2.at(i));
 313                 return newVec;
 314         }
 315         catch(exception& e) {
 316                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 317                 exit(1);
 318         }
 319         catch(...) {
 320                 cout << "An unknown error has occurred in the VecCalc class function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 321                 exit(1);
 322         }
 323 }
 324 /***********************************************************************/
 325 vector<double> VecCalc::remDup(vector<double> vec){
 326         try {
 327                 vector<double> newVec;
 328                 newVec.push_back(vec.at(0));
 329                 for(int i = 1; i < vec.size(); i++)
 330                         if(vec.at(i) != vec.at(i-1))
 331                                 newVec.push_back(vec.at(i));
 332                 return newVec;
 333         }
 334         catch(exception& e) {
 335                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 336                 exit(1);
 337         }
 338         catch(...) {
 339                 cout << "An unknown error has occurred in the VecCalc class function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 340                 exit(1);
 341         }
 342 }
 343 /***********************************************************************/
 344 vector<double> VecCalc::genCVec(vector<double> vec1){
 345         try {
 346                 vector<double> vec2;
 347                 double sum = 0;
 348                 for(int i = 0; i < vec1.size(); i++)
 349                 {
 350                         sum += vec1.at(i);
 351                         vec2.push_back(sum);
 352                 }
 353                 return vec2;
 354         }
 355         catch(exception& e) {
 356                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 357                 exit(1);
 358         }
 359         catch(...) {
 360                 cout << "An unknown error has occurred in the VecCalc class function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 361                 exit(1);
 362         }
 363 }
 364 /***********************************************************************/
 365 vector<double> VecCalc::genRelVec(vector<double> vec){
 366         try {
 367                 vector<double> newVec;
 368                 double sum = sumElements(vec);
 369                 for(int i = 0; i < vec.size(); i++)
 370                         newVec.push_back(vec.at(i)/sum);
 371                 return newVec;
 372         }
 373         catch(exception& e) {
 374                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 375                 exit(1);
 376         }
 377         catch(...) {
 378                 cout << "An unknown error has occurred in the VecCalc class function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 379                 exit(1);
 380         }
 381 }
 382 /***********************************************************************/
 383 vector<double> VecCalc::genDiffVec(vector<double> vec1, vector<double> vec2){
 384         try {
 385                 vector<double> newVec;
 386                 for(int i = 0; i < vec1.size(); i++)
 387                         newVec.push_back(vec1.at(i)-vec2.at(i));
 388                 return newVec;
 389         }
 390         catch(exception& e) {
 391                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 392                 exit(1);
 393         }
 394         catch(...) {
 395                 cout << "An unknown error has occurred in the VecCalc class function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 396                 exit(1);
 397         }
 398 }
 399 /***********************************************************************/
 400 vector<double> VecCalc::genCSVec(vector<double> vec){
 401         try {
 402                 vector<double> newVec;
 403                 double curInd = vec.at(vec.size()-1);
 404                 double sumSpec = 0;
 405                 double cSumSpec = 0;
 406                 for(int i = vec.size()-1; i >= 0; i--)
 407                 {
 408                         if(vec.at(i) == curInd)
 409                                 sumSpec++;
 410                         else
 411                         {
 412                                 cSumSpec += sumSpec;
 413                                 newVec.push_back(cSumSpec);
 414                                 sumSpec = 1;
 415                                 curInd = vec.at(i);
 416                         }
 417                         if(i == 0)
 418                                 newVec.push_back(cSumSpec + sumSpec);
 419                 }
 420                 return newVec;
 421         }
 422         catch(exception& e) {
 423                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 424                 exit(1);
 425         }
 426         catch(...) {
 427                 cout << "An unknown error has occurred in the VecCalc class function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 428                 exit(1);
 429         }
 430 }
 431 /***********************************************************************/
 432 vector<double> VecCalc::genTotVec(vector<vector<double> > vec){
 433         try {
 434                 vector<double> newVec = vec.at(0);
 435                 for(int i = 1; i < vec.size(); i++)
 436                         newVec = addVecs(newVec, vec.at(i));
 437                 return newVec;
 438         }
 439         catch(exception& e) {
 440                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 441                 exit(1);
 442         }
 443         catch(...) {
 444                 cout << "An unknown error has occurred in the VecCalc class function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 445                 exit(1);
 446         }
 447 }
 448 /***********************************************************************/
 449 vector<double> VecCalc::quicksort(vector<double> vec){
 450         try {
 451                 sort(vec.rbegin(), vec.rend());
 452                 return vec;
 453         }
 454         catch(exception& e) {
 455                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 456                 exit(1);
 457         }
 458         catch(...) {
 459                 cout << "An unknown error has occurred in the VecCalc class function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 460                 exit(1);
 461         }
 462 }
 463 /***********************************************************************/
 464 vector<vector<double> > VecCalc::gen2DVec(vector<double> vec, int div, int type){
 465         try {
 466                 vector<vector<double> > vec2D;
 467                 int gap = vec.size()/div;
 468                 for(int i = 0; i < div; i++)
 469                 {
 470                         vector<double> inVec;
 471                         for(int j = 0; j < gap; j++)
 472                                 if(type == 0)
 473                                         inVec.push_back(vec.at(j + i*gap)); //Rows
 474                                 else
 475                                         inVec.push_back(vec.at(i + j*div)); //Columns
 476                         vec2D.push_back(inVec);
 477                 }
 478                 return vec2D;
 479         }
 480         catch(exception& e) {
 481                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 482                 exit(1);
 483         }
 484         catch(...) {
 485                 cout << "An unknown error has occurred in the VecCalc class function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 486                 exit(1);
 487         }
 488 }
 489 /***********************************************************************/
 490 vector<string> VecCalc::getSData(char name[]){
 491         try {
 492                 vector<string> vec;
 493                 string data;
 494                 ifstream file(name);
 495                 if(file.is_open())
 496                 {
 497                         while(file >> data)
 498                                 vec.push_back(data);
 499                         file.close();
 500                 }
 501                 return vec;
 502         }
 503         catch(exception& e) {
 504                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 505                 exit(1);
 506         }
 507         catch(...) {
 508                 cout << "An unknown error has occurred in the VecCalc class function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 509                 exit(1);
 510         }
 511 }
 512 /***********************************************************************/
 513 double BDiversity::getWhitt(vector<double> vec1, vector<double> vec2){
 514         try {
 515                 VecCalc vecCalc;
 516                 double numSpec = vecCalc.numNZ(vecCalc.addVecs(vec1,vec2));
 517                 return 2*numSpec/(vecCalc.numNZ(vec1)+vecCalc.numNZ(vec2))-1;
 518         }
 519         catch(exception& e) {
 520                 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getWhitt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 521                 exit(1);
 522         }
 523         catch(...) {
 524                 cout << "An unknown error has occurred in the BDiversity class function getWhitt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 525                 exit(1);
 526         }
 527 }
 528 /***********************************************************************/
 529 double BDiversity::getMS(vector<double> vec1, vector<double> vec2){
 530         try {
 531                 VecCalc vecCalc;
 532                 double a = vecCalc.numNZ(vecCalc.multVecs(vec1, vec2));
 533                 double b = vecCalc.numPos(vecCalc.addVecs(vecCalc.addVecs(vec1, vecCalc.multiply(vec2, -1)), vecCalc.multiply(vecCalc.multVecs(vec1, vec2), -1)));
 534                 double c = vecCalc.numPos(vecCalc.addVecs(vecCalc.addVecs(vec2, vecCalc.multiply(vec1, -1)), vecCalc.multiply(vecCalc.multVecs(vec2, vec1), -1)));
 535                 return a/(a+b+c);
 536         }
 537         catch(exception& e) {
 538                 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getMS. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 539                 exit(1);
 540         }
 541         catch(...) {
 542                 cout << "An unknown error has occurred in the BDiversity class function getMS. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 543                 exit(1);
 544         }
 545 }
 546 /***********************************************************************/
 547 double BDiversity::getSor(vector<double> vec1, vector<double> vec2){
 548         try {
 549                 double sum = 0;
 550                 double asum = 0;
 551                 double bsum = 0;
 552                 for(int i = 0; i < vec1.size(); i++)
 553                 {
 554                         asum += vec1.at(i);
 555                         bsum += vec2.at(i);
 556                         if(vec1.at(i) >= vec2.at(i))
 557                                 sum += vec2.at(i);
 558                         else
 559                                 sum += vec1.at(i);
 560                 }
 561                 return 2*sum/(asum+bsum);
 562         }
 563         catch(exception& e) {
 564                 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getSor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 565                 exit(1);
 566         }
 567         catch(...) {
 568                 cout << "An unknown error has occurred in the BDiversity class function getSor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 569                 exit(1);
 570         }
 571 }
 572 /***********************************************************************/
 573 double BDiversity::getMor(vector<double> vec1, vector<double> vec2){
 574         try {
 575                 double sum = 0;
 576                 double asum = 0;
 577                 double bsum = 0;
 578                 double sum1 = 0;
 579                 double sum2 = 0;
 580                 for(int i = 0; i < vec1.size(); i++)
 581                 {
 582                         sum += vec1.at(i)*vec2.at(i);
 583                         asum += pow(vec1.at(i),2);
 584                         bsum += pow(vec2.at(i),2);
 585                         sum1 += vec1.at(i);
 586                         sum2 += vec2.at(i);
 587                 }
 588                 return 2*sum/((asum/pow(sum1,2)+bsum/pow(sum2,2))*(sum1*sum2));
 589         }
 590         catch(exception& e) {
 591                 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getMor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 592                 exit(1);
 593         }
 594         catch(...) {
 595                 cout << "An unknown error has occurred in the BDiversity class function getMor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 596                 exit(1);
 597         }
 598 }
 599 /***********************************************************************/
 600 void BDiversity::printD(vector<vector<double> > columns, int type){
 601         try {
 602                 cout << "         ";
 603                 for(int i = 0; i < columns.size(); i++)
 604                         cout << "Column:" << i << "    ";
 605                 cout << "\n";
 606                 for(int i = 0; i < columns.size(); i++)
 607                 {
 608                         cout << "Column " << i << ":";
 609                         for(int j = 0; j < columns.size(); j++)
 610                         {
 611                                 if(j > i)
 612                                 {
 613                                         double B;
 614                                         if(type == 1)
 615                                                 B = getWhitt(columns.at(i), columns.at(j));
 616                                         else if(type == 2)
 617                                                 B = 1-getMS(columns.at(i), columns.at(j));
 618                                         else if(type == 3)
 619                                                 B = 1-getSor(columns.at(i), columns.at(j));
 620                                         else if(type == 4)
 621                                                 B = 1-getMor(columns.at(i), columns.at(j));
 622
 623                                         cout << B << "    ";
 624                                 }
 625                                 else
 626                                         cout << "            ";
 627                         }
 628                         cout << "\n";
 629                 }
 630         }
 631         catch(exception& e) {
 632                 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function printD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 633                 exit(1);
 634         }
 635         catch(...) {
 636                 cout << "An unknown error has occurred in the BDiversity class function printD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 637                 exit(1);
 638         }
 639 }
 640 /***********************************************************************/
 641 void BDiversity::doBD(vector<double> vec, double cols)//vec = The data vector if the data table was read left-to-right, top-to-bottom. cols = the number of columns in the data table.
 642 {       try {
 643                 VecCalc vecCalc;
 644                 vector<vector<double> > columns = vecCalc.gen2DVec(vec,cols,1);
 645
 646                 cout.setf(ios_base::fixed, ios_base::floatfield);
 647                 cout.precision(6);//This formats the data so the tables look pretty.
 648
 649                 //Whittaker's measure Bw (presence/absence data)
 650                 cout << "Whittaker's measure Bw (presence/absence data):\n";
 651                 printD(columns, 1);
 652                 double sum = 0;
 653                 for(int i = 0; i < cols; i++)
 654                         sum += vecCalc.numNZ(columns.at(i));
 655                 double meanRichness = sum/cols;
 656                 vector<double> totVec = vecCalc.genTotVec(columns);
 657                 double totRichness = vecCalc.numNZ(totVec);
 658                 cout << "\nOverall B Diversity = " << totRichness/meanRichness << "\n\n\n";//The overall B Diversity
 659
 660                 //Marczewski-Steinhaus(MS) distance(Jaccard index) (presence/abscence data)
 661                 cout << "Marczewski-Steinhaus(MS) distance(Jaccard index) (presence/abscence data):\n";
 662                 printD(columns, 2);
 663                 sum = 0;
 664                 for(int i = 0; i < cols; i++)
 665                         for(int j = 0; j < cols; j++)
 666                                 if(j > i)
 667                                         sum += vecCalc.numNZ(columns.at(i))+vecCalc.numNZ(columns.at(j)) - 2*vecCalc.numNZ(vecCalc.multVecs(columns.at(i), columns.at(j)));
 668                 cout << "\nOverall B Diversity = " << sum/cols << "\n\n\n";//The overall B Diversity
 669
 670                 //Sorensen quantitative index (abundance data)
 671                 cout << "Sorensen quantitative index (abundance data):\n";
 672                 printD(columns, 3);
 673                 cout << "\n\n";
 674
 675                 //Sorensen quantitative index (abundance data)
 676                 cout << "Morisita-Horn index (abundance data):\n";
 677                 printD(columns,4);
 678         }
 679         catch(exception& e) {
 680                 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function doBD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 681                 exit(1);
 682         }
 683         catch(...) {
 684                 cout << "An unknown error has occurred in the BDiversity class function doBD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 685                 exit(1);
 686         }
 687 }
 688
 689 /***********************************************************************/
 690 void BrokenStick::doBStick(vector<double> vec)//vec = The data vector.
 691 {       try {
 692                 VecCalc vecCalc;
 693                 vec = vecCalc.quicksort(vec);
 694                 double numInd = vecCalc.sumElements(vec);
 695                 double numSpec = vecCalc.numNZ(vec);
 696                 vector<double> cObsVec = vecCalc.genCVec(vec);
 697                 vector<double> cExpVec;
 698                 vec = vecCalc.power(vec, -1);
 699                 double sumExp = 0;
 700                 for(int i = 0; i < numSpec; i++)
 701                 {
 702                         double n = 0;
 703                         for(int j = i; j < numSpec; j++)
 704                                 n += 1/(double)(j+1);
 705
 706                         sumExp += numInd/numSpec*n;
 707                         cExpVec.push_back(sumExp);
 708                 }
 709
 710                 //Statistical analysis
 711                 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
 712                 double DStatistic = maxDiff/numInd;
 713
 714                 cout << "D-Statistic = " << DStatistic << "\n";
 715                 if(vec.size() > 30)
 716                         cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
 717         }
 718         catch(exception& e) {
 719                 cout << "Standard Error: " << e.what() << " has occurred in the BrokenStick class Function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 720                 exit(1);
 721         }
 722         catch(...) {
 723                 cout << "An unknown error has occurred in the BrokenStick class function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 724                 exit(1);
 725         }
 726 }
 727
 728 /***********************************************************************/
 729 double kEq(double k, double spec)
 730 {
 731         return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec));
 732 }
 733 /***********************************************************************/
 734 void GeometricSeries::doGeomTest(vector<double> vec)//vec = the data vector
 735 {       try {
 736                 VecCalc vecCalc;
 737                 vec = vecCalc.quicksort(vec);//Sorts vec
 738                 double numInd = vecCalc.sumElements(vec);//numInd = The total number of individuals in the data set.
 739                 double numSpec = vecCalc.numNZ(vec);//numSpec = the number of nonzero elements in the data set.
 740                 double min = -1.0*vecCalc.findMax(vecCalc.multiply(vec, -1));
 741                 double k = .5;
 742                 double step = .49999;
 743                 while(fabs(min - numInd*kEq(k, numSpec)) > .0001)//This uses a binary search to find the value of k.
 744                 {
 745                         if(numInd*kEq(k, numSpec) > min)
 746                                 k += step;
 747                         else
 748                                 k -= step;
 749                         step /= 2;
 750                 }
 751                 cout << "k = " << k << "\n";
 752                 double cK = 1/(1-pow(1-k, numSpec));
 753
 754                 vector<double> cObsVec = vecCalc.genCVec(vec);
 755                 vector<double> cExpVec;
 756                 double sumExp = 0;
 757                 for(int i = 0; i < vec.size(); i++)
 758                 {
 759                         sumExp += numInd*cK*k*pow(1-k, i);
 760                         cExpVec.push_back(sumExp);
 761                 }
 762                 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
 763                 double DStatistic = maxDiff/numInd;
 764
 765                 //Statistical Analysis
 766                 cout << "D-Statistic = " << DStatistic << "\n";
 767                 if(vec.size() > 30)
 768                         cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
 769         }
 770         catch(exception& e) {
 771                 cout << "Standard Error: " << e.what() << " has occurred in the GeometricSeries class Function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 772                 exit(1);
 773         }
 774         catch(...) {
 775                 cout << "An unknown error has occurred in the GeometricSeries class function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 776                 exit(1);
 777         }
 778 }
 779
 780 /***********************************************************************/
 781 void Jackknifing::doJK(vector<double> vec, double cols)//vec = the data vector if the data table was read left-to-right, top-to-bottom. cols = The number of columns in the data table.
 782 {       try {
 783                 VecCalc vecCalc;
 784                 SSBPDiversityIndices ssbp;
 785                 double rows = vec.size()/cols;
 786                 vector<vector<double> > species = vecCalc.gen2DVec(vec, rows, 0);//converts vec into a 2 dimensional vector
 787                 vector<double> sumRows;
 788                 vector<double> pVals;
 789
 790                 for(int i = 0; i < rows; i++)
 791                         sumRows.push_back(vecCalc.sumElements(species.at(i)));
 792                 double st = 1/ssbp.getSimp(sumRows);//The observed estimate using the Simpson Index. Can be changed to use other indexes of diversity.
 793                 for(int i = 0; i < cols; i++)
 794                 {
 795                         vector<double> newVec;
 796                         for(int j = 0; j < rows; j++)
 797                                 newVec.push_back(vecCalc.sumElements(species.at(j))-species.at(j).at(i));
 798                         pVals.push_back(cols*st-((cols-1)/ssbp.getSimp(newVec)));
 799                 }
 800
 801                 double mean = vecCalc.mean(pVals);
 802                 double stErr = vecCalc.stError(pVals);//stErr = standard error
 803                 TDTable table;
 804                 double confidence = .95;
 805                 double confLimit;
 806                 cout << "Confidence Level (.8, .9, .95, .98, .99, .998, .999): ";
 807                 cin >> confidence;
 808                 double confLevels[] = {.80,.90,.95,.98,.99,.998,.999};
 809                 for(int i = 0; i < 7; i++)
 810                         if(confidence == confLevels[i])
 811                         {
 812                                 confLimit = table.getConfLimit(cols-2,i);
 813                                 i = 7;
 814                         }
 815
 816                 //Statistical Analysis
 817                 cout << "Lower limit = " << mean - confLimit*stErr << "\n";
 818                 cout << "Upper limit = " << mean + confLimit*stErr << "\n";
 819                 cout << "Observed estimate = " << st << "\n\n";
 820                 cout << "Jackknifed estimate = " << mean << "\n";
 821         }
 822         catch(exception& e) {
 823                 cout << "Standard Error: " << e.what() << " has occurred in the Jackknifing class Function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 824                 exit(1);
 825         }
 826         catch(...) {
 827                 cout << "An unknown error has occurred in the Jackknifing class function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 828                 exit(1);
 829         }
 830 }
 831 /***********************************************************************/
 832 void KS2SampleTest::doKSTest(vector<double> abun1, vector<double> abun2)//abun1 = 1st vector of abundances, abun2 = 2nd vector of abundances
 833 {       try {
 834                 VecCalc vecCalc;
 835                 double numNZ1 = vecCalc.numNZ(abun1);//Number of nonzero elements in abun1
 836                 double numNZ2 = vecCalc.numNZ(abun2);//Number of nonzero elements in abun2
 837                 abun1 = vecCalc.quicksort(abun1);//Sorts abun1
 838                 abun2 = vecCalc.quicksort(abun2);//Sorts abun2
 839                 abun1 = vecCalc.genRelVec(abun1);//Generates the relative vector for abun1
 840                 abun2 = vecCalc.genRelVec(abun2);//Generates the relative vector for abun2
 841
 842                 abun1 = vecCalc.genCVec(abun1);//Generates the cumulative vector for abun1
 843                 abun2 = vecCalc.genCVec(abun2);//Generates the cumulative vector for abun2
 844
 845                 double maxDiff = vecCalc.findMaxDiff(abun1, abun2);
 846                 double DStatistic = maxDiff*numNZ1*numNZ2;
 847
 848                 cout << "Null Hypothesis = There is no difference.\n";
 849                 cout << "D-Statistic = " << DStatistic << "\n";
 850
 851                 double a = pow((numNZ1 + numNZ2)/(numNZ1*numNZ2),.5);
 852                 double pVal = exp(-2*pow(maxDiff/a,2));
 853
 854                 if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25.
 855                         cout << "P-Value = " << pVal << "\n\n";
 856                 else
 857                 {
 858                         cout << "D-Statistic must be higher than the critical value to reject the null hypothesis.\n" ;
 859                         cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n";
 860                         cout << "95% Confidence Critical Value = " << 1.36*a*numNZ1*numNZ2 << "\n";
 861                 }
 862         }
 863         catch(exception& e) {
 864                 cout << "Standard Error: " << e.what() << " has occurred in the KS2SampleTest class Function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 865                 exit(1);
 866         }
 867         catch(...) {
 868                 cout << "An unknown error has occurred in the KS2SampleTest class function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 869                 exit(1);
 870         }
 871 }
 872
 873 /***********************************************************************/
 874 double logS(double num)
 875 {
 876         return -(1-num)*log(1-num)/num;
 877 }
 878
 879 /***********************************************************************/
 880 void LogSD::doLogSD(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
 881 {       try {
 882                 VecCalc vecCalc;
 883                 double numSpec = vecCalc.sumElements(specVec);//numSpec = The total number of species
 884                 cout << "number of species = " << numSpec << "\n";
 885                 double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec, specVec));
 886                 double snRatio = numSpec/numInd;
 887                 double x = .5;
 888                 double step = .4999999999;
 889                 while(fabs(snRatio - logS(x)) > .00001) //This uses a binary search to find the value of x.
 890                 {
 891                         if(logS(x) > snRatio)
 892                                 x += step;
 893                         else
 894                                 x -= step;
 895                         step /= 2;
 896                 }
 897                 double alpha = numInd*(1-x)/x;
 898
 899                 int ind;
 900                 cout << "Number of individuals:"; //Ask the user for the number of individuals.
 901                 cin >> ind;
 902                 double spec = alpha*pow(x, ind)/ind;
 903                 cout << "Number of species expected = " << spec << "\n" << "X value = " << x << "\n" << "Alpha value= " << alpha << "\n";//Outputs the number of species expected with the given number of individuals.
 904
 905                 vector<double> obsSpec;
 906                 vector<double> cObsSpec;
 907                 vector<double> expSpec;
 908                 vector<double> cExpSpec;
 909                 vector<double> cDiff;
 910
 911                 // Generates the cumulative observed species vector.
 912                 int oct = 1;
 913                 double octSumObs = 0;
 914                 for(int y = 0; y < specVec.size(); y++)
 915                 {
 916                         if(indVec.at(y) - .5 < pow(2.0, oct))
 917                                 octSumObs += specVec.at(y);
 918                         else
 919                         {
 920                                 obsSpec.push_back(octSumObs);
 921                                 octSumObs = specVec.at(y);
 922                                 oct++;
 923                         }
 924                         if(y == specVec.size()-1)
 925                                 obsSpec.push_back(octSumObs);
 926                 }
 927                 cObsSpec = vecCalc.genCVec(obsSpec);
 928                 cObsSpec = vecCalc.add(cObsSpec,-.5);
 929
 930                 // Generates the cumulative expected species vector.
 931                 oct = 1;
 932                 double octSumExp = 0;
 933                 for(int g = 1; g <= indVec.at(indVec.size()-1); g++)
 934                 {
 935                         if(g - .5 < pow(2.0, oct))
 936                                 octSumExp += alpha*pow(x,g)/(g);
 937                         else
 938                         {
 939                                 expSpec.push_back(octSumExp);
 940                                 octSumExp = alpha*pow(x,g)/(g);
 941                                 oct++;
 942                         }
 943                         if(g == indVec.at(indVec.size()-1))
 944                                 expSpec.push_back(octSumExp);
 945                 }
 946                 cExpSpec = vecCalc.genCVec(expSpec);
 947
 948                 // Statistical Analysis
 949                 double dTStat = vecCalc.findDStat(cObsSpec, cExpSpec, numSpec);
 950                 cout << "D Test Statistic = " << dTStat << "\n";
 951                 cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
 952                 cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
 953         }
 954         catch(exception& e) {
 955                 cout << "Standard Error: " << e.what() << " has occurred in the LogSD class Function doLogSD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 956                 exit(1);
 957         }
 958         catch(...) {
 959                 cout << "An unknown error has occurred in the LogSD class function doLogSD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 960                 exit(1);
 961         }
 962 }
 963 /***********************************************************************/
 964 void QStatistic::doQStat(vector<double> vec)//vec = The data vector.
 965 {       try {
 966                 VecCalc vecCalc;
 967                 vector<double> cVec = vecCalc.genCSVec(vec);
 968                 vector<int> q = vecCalc.findQuartiles(cVec);//q.at(0) is the index of the first quartile; q.at(1) is the index of the second quartile.
 969                 vector<double> nDupVec = vecCalc.remDup(vec);//nDupVec only contains one of every unique element in cVec.
 970                 double Q;
 971                 if(q.at(0) != 0)//The case if neither quartile is 0 or 1
 972                         Q = (.5*(cVec.at(q.at(0))-cVec.at(q.at(0)-1)) + (cVec.at(q.at(1)-1)-cVec.at(q.at(0))) + .5*(cVec.at(q.at(1))-cVec.at(q.at(1)-1)))/log(nDupVec.at(nDupVec.size()-1-q.at(1))/nDupVec.at(nDupVec.size()-1-q.at(0)));
 973                 else if(q.at(0) == 0 && (q.at(1) == 0 || q.at(1) == 1))//The case if the quartiles are both at index 0 or one is at 0 and the other at 1.
 974                         Q = (.5*cVec.at(0) + .5*(cVec.at(1)-cVec.at(0)))/log(nDupVec.at(nDupVec.size()-2)/nDupVec.at(nDupVec.size()-1));
 975                 else if(q.at(0) == 0 && q.at(1) > 1) //The case if the lower quartile is at index 0 and upper quartile index is above index 1.
 976                         Q = (.5*cVec.at(0) + (cVec.at(q.at(1)-1)-cVec.at(q.at(0))) + .5*(cVec.at(q.at(1))-cVec.at(q.at(1)-1)))/log(nDupVec.at(nDupVec.size()-1-q.at(1))/nDupVec.at(nDupVec.size()-1));
 977
 978                 cout << "Q = " << Q << "\n";
 979         }
 980         catch(exception& e) {
 981                 cout << "Standard Error: " << e.what() << " has occurred in the QStatistic class Function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 982                 exit(1);
 983         }
 984         catch(...) {
 985                 cout << "An unknown error has occurred in the QStatistic class function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 986                 exit(1);
 987         }
 988 }
 989 /***********************************************************************/
 990 double SSBPDiversityIndices::getShan(vector<double> vec)//vec = The data vector.
 991 {       try {
 992                 VecCalc vecCalc;
 993                 double nz = vecCalc.numNZ(vec);
 994                 double nSum = vecCalc.sumElements(vec);
 995                 double H = 0;
 996                 for(int i = 0; i < nz; i++)
 997                         H += vec.at(i)/nSum*log(vec.at(i)/nSum);
 998                 H *= -1;
 999                 return H;
1000         }
1001         catch(exception& e) {
1002                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1003                 exit(1);
1004         }
1005         catch(...) {
1006                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1007                 exit(1);
1008         }
1009 }
1010 /***********************************************************************/
1011 double SSBPDiversityIndices::getSimp(vector<double> vec)//vec = The data vector.
1012 {       try {
1013                 VecCalc vecCalc;
1014                 double nSum = vecCalc.sumElements(vec);
1015                 double D = 0;
1016                 for(int j = 0; j < vec.size(); j++)
1017                         D += vec.at(j)*(vec.at(j)-1)/(nSum*(nSum-1));
1018                 return D;
1019         }
1020         catch(exception& e) {
1021                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1022                 exit(1);
1023         }
1024         catch(...) {
1025                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1026                 exit(1);
1027         }
1028 }
1029 /***********************************************************************/
1030 double SSBPDiversityIndices::getBP(vector<double> vec)//vec = The data vector.
1031 {       try {
1032                 VecCalc vecCalc;
1033                 double nSum = vecCalc.sumElements(vec);
1034                 return vecCalc.findMax(vec)/nSum;
1035         }
1036         catch(exception& e) {
1037                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1038                 exit(1);
1039         }
1040         catch(...) {
1041                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1042                 exit(1);
1043         }
1044 }
1045 /***********************************************************************/
1046 void SSBPDiversityIndices::doSSBP(vector<double> vec)//vec = The data vector.
1047 {       try {
1048                 VecCalc vecCalc;
1049                 double nz = vecCalc.numNZ(vec);
1050
1051                 //Shannon index
1052                 double H = getShan(vec);
1053                 cout << "H = " << H << "\n";
1054                 cout << "Eveness = " << H/log(nz) << "\n\n";
1055
1056                 //Simpson index
1057                 double D = getSimp(vec);
1058                 cout << "D diversity = " << 1/D << "\n";
1059                 cout << "Eveness = " << 1/D/nz << "\n\n";
1060
1061                 //Berger-Parker index
1062                 double BP = getBP(vec);
1063                 cout << "BP index = " << BP << "\n\n";
1064         }
1065         catch(exception& e) {
1066                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1067                 exit(1);
1068         }
1069         catch(...) {
1070                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1071                 exit(1);
1072         }
1073 }
1074 /***********************************************************************/
1075 double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom
1076 {    try {
1077                 //Confidence Level        .80    .90     .95     .98     .99     .998     .999
1078                 double values[33][7] = {{3.078, 6.314,  12.706, 31.821, 63.656, 318.289, 636.578},
1079                                                     {1.886,     2.920,  4.303,  6.965,  9.925,  22.328, 31.600},
1080                                                         {1.638, 2.353,  3.182,  4.541,  5.841,  10.214, 12.924},
1081                                                         {1.533, 2.132,  2.776,  3.747,  4.604,  7.173,  8.610},
1082                                                         {1.476, 2.015,  2.571,  3.365,  4.032,  5.894,  6.869},
1083                                                         {1.440, 1.943,  2.447,  3.143,  3.707,  5.208,  5.959},
1084                                                         {1.415, 1.895,  2.365,  2.998,  3.499,  4.785,  5.408},
1085                                                         {1.397, 1.860,  2.306,  2.896,  3.355,  4.501,  5.041},
1086                                                         {1.383, 1.833,  2.262,  2.821,  3.250,  4.297,  4.781},
1087                                                         {1.372, 1.812,  2.228,  2.764,  3.169,  4.144,  4.587},
1088                                                         {1.363, 1.796,  2.201,  2.718,  3.106,  4.025,  4.437},
1089                                                         {1.356, 1.782,  2.179,  2.681,  3.055,  3.930,  4.318},
1090                                                         {1.350, 1.771,  2.160,  2.650,  3.012,  3.852,  4.221},
1091                                                         {1.345, 1.761,  2.145,  2.624,  2.977,  3.787,  4.140},
1092                                                         {1.341, 1.753,  2.131,  2.602,  2.947,  3.733,  4.073},
1093                                                         {1.337, 1.746,  2.120,  2.583,  2.921,  3.686,  4.015},
1094                                                         {1.333, 1.740,  2.110,  2.567,  2.898,  3.646,  3.965},
1095                                                         {1.330, 1.734,  2.101,  2.552,  2.878,  3.610,  3.922},
1096                                                         {1.328, 1.729,  2.093,  2.539,  2.861,  3.579,  3.883},
1097                                                         {1.325, 1.725,  2.086,  2.528,  2.845,  3.552,  3.850},
1098                                                         {1.323, 1.721,  2.080,  2.518,  2.831,  3.527,  3.819},
1099                                                         {1.321, 1.717,  2.074,  2.508,  2.819,  3.505,  3.792},
1100                                                         {1.319, 1.714,  2.069,  2.500,  2.807,  3.485,  3.768},
1101                                                         {1.318, 1.711,  2.064,  2.492,  2.797,  3.467,  3.745},
1102                                                         {1.316, 1.708,  2.060,  2.485,  2.787,  3.450,  3.725},
1103                                                         {1.315, 1.706,  2.056,  2.479,  2.779,  3.435,  3.707},
1104                                                         {1.314, 1.703,  2.052,  2.473,  2.771,  3.421,  3.689},
1105                                                         {1.313, 1.701,  2.048,  2.467,  2.763,  3.408,  3.674},
1106                                                         {1.311, 1.699,  2.045,  2.462,  2.756,  3.396,  3.660},
1107                                                         {1.310, 1.697,  2.042,  2.457,  2.750,  3.385,  3.646},
1108                                                         {1.296, 1.671,  2.000,  2.390,  2.660,  3.232,  3.460},
1109                                                         {1.289, 1.658,  1.980,  2.358,  2.617,  3.160,  3.373},
1110                                                         {1.282, 1.645,  1.960,  2.326,  2.576,  3.091,  3.291}};
1111                 return values[row][col];
1112         }
1113         catch(exception& e) {
1114                 cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1115                 exit(1);
1116         }
1117         catch(...) {
1118                 cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1119                 exit(1);
1120         }
1121 }
1122 /***********************************************************************
1123 void TrunLN::doTrunLN(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
1124 {
1125         VecCalc vecCalc;
1126         double numSpec = vecCalc.sumElements(specVec); //numSpec = The total number of species
1127         double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec,specVec)); //numInd = The total number of individuals
1128
1129         double obsMean = 0;
1130         for(int i = 0; i < indVec.size(); i++)
1131                 obsMean += log10(indVec.at(i));
1132         obsMean /= numSpec; //obsMean = observed mean of the individuals vector
1133         cout << "obsMean = " << obsMean << "\n";
1134         double variance = 0;
1135         for(int t = 0; t < indVec.size(); t++)
1136                 variance += pow(log10(indVec.at(t))-obsMean,2)/numSpec;
1137
1138          double rO = 0;
1139          for(int k = 0; k < indVec.size(); k++)
1140                 rO += log10(indVec.at(k));
1141          rO /= indVec.size();
1142          double veilLine = .5;//The desired veil line.
1143          double auxFunc = -(obsMean-rO)/(obsMean-log10(veilLine));
1144          double uX = obsMean-auxFunc*(obsMean-log10(veilLine));
1145          double vX = variance + auxFunc*pow(obsMean-log10(veilLine),2);
1146          double z = (log10(veilLine)-uX)/pow(vX, .5);
1147          double p = .5*(erf(z)+1);
1148          double specRichness = numSpec/(1-p);
1149
1150          double numUnseen = .5*(erf((log10(.5)-uX)/pow(vX,.5))+1)*specRichness;
1151
1152
1153          vector<double> cExp;
1154          for(int i = 1; i < 8; i++)
1155          {
1156                 double a = pow(10, i-1)+.5;
1157                 double b = log10(a);
1158                 double c = (b - uX)/pow(vX,.5);
1159                 double d = .5*(erf(c)+1);
1160                 double numS = d*specRichness;
1161                 double toPush = numS - numUnseen;
1162                 cExp.push_back(toPush);
1163         }
1164         vector<double> cObs;
1165         double sumOct = 0;
1166         for(int i = 0; i < 8; i++)
1167         {
1168                 sumOct = 0;
1169                 for(int r = 0; r < indVec.size(); r++)
1170                 {
1171                         if(indVec.at(r) < pow(10, i-1)+.5)
1172                                 sumOct += specVec.at(r);
1173                         else
1174                         {
1175                                 cObs.push_back(sumOct);
1176                                 sumOct = specVec.at(r);
1177                                 r = indVec.size();
1178                         }
1179                         if(r == indVec.size()-1)
1180                                 cObs.push_back(sumOct);
1181                 }
1182         }
1183
1184         //Statistical Analysis
1185         double d = vecCalc.findDStat(cExp, cObs, numSpec);
1186         cout << "DStat = " << d << "\n";
1187         cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
1188         cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
1189 }
1190 /***********************************************************************/