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 int VecCalc::sumElements(vector<int> vec){
  87         try {
  88                 return sumElements(vec,0);
  89         }
  90         catch(exception& e) {
  91                 errorOut(e, "VecCalc", "sumElements");
  92                 exit(1);
  93         }
  94 }
  95 /***********************************************************************/
  96 int VecCalc::sumElements(vector<int> vec, int index){
  97         try {
  98                 int sum = 0;
  99                 for(int i = index; i < vec.size(); i++)
 100                         sum += vec.at(i);
 101                 return sum;
 102         }
 103         catch(exception& e) {
 104                 errorOut(e, "VecCalc", "sumElements");
 105                 exit(1);
 106         }
 107 }
 108
 109 /***********************************************************************/
 110 double VecCalc::sumElements(vector<double> vec){
 111         try {
 112                 double sum = 0;
 113                 for(int i = 0; i < vec.size(); i++)
 114                         sum += vec.at(i);
 115                 return sum;
 116         }
 117         catch(exception& e) {
 118                 errorOut(e, "VecCalc", "sumElements");
 119                 exit(1);
 120         }
 121 }
 122 /***********************************************************************/
 123 double VecCalc::sumElements(vector<double> vec, int index){
 124         try {
 125                 double sum = 0;
 126                 for(int i = index; i < vec.size(); i++)
 127                         sum += vec.at(i);
 128                 return sum;
 129         }
 130         catch(exception& e) {
 131                 errorOut(e, "VecCalc", "sumElements");
 132                 exit(1);
 133         }
 134 }
 135 /***********************************************************************
 136 double VecCalc::findMax(vector<double> vec){
 137         try {
 138                 double max = -1000000.0;
 139                 for(int i = 0; i < vec.size(); i++)
 140                         if(vec.at(i) > max)
 141                                 max = vec.at(i);
 142                 return max;
 143         }
 144         catch(exception& e) {
 145                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 146                 exit(1);
 147         }
 148         catch(...) {
 149                 cout << "An unknown error has occurred in the VecCalc class function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 150                 exit(1);
 151         }
 152 }
 153 /***********************************************************************/
 154 int VecCalc::numNZ(vector<int> vec){
 155         try {
 156                 int numNZ = 0;
 157                 for(int i = 0; i < vec.size(); i++)
 158                         if(vec.at(i) != 0)
 159                                 numNZ++;
 160                 return numNZ;
 161         }
 162         catch(exception& e) {
 163                 errorOut(e, "VecCalc", "numNZ");
 164                 exit(1);
 165         }
 166 }
 167 /***********************************************************************/
 168 double VecCalc::numNZ(vector<double> vec){
 169         try {
 170                 double numNZ = 0;
 171                 for(int i = 0; i < vec.size(); i++)
 172                         if(vec.at(i) != 0)
 173                                 numNZ++;
 174                 return numNZ;
 175         }
 176         catch(exception& e) {
 177                 errorOut(e, "VecCalc", "numNZ");
 178                 exit(1);
 179         }
 180 }
 181 /***********************************************************************
 182 double VecCalc::numPos(vector<double> vec){
 183         try {
 184                 double numPos = 0;
 185                 for(int i = 0 ; i < vec.size(); i++)
 186                         if(vec.at(i) > 0)
 187                                 numPos++;
 188                 return numPos;
 189         }
 190         catch(exception& e) {
 191                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 192                 exit(1);
 193         }
 194         catch(...) {
 195                 cout << "An unknown error has occurred in the VecCalc class function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 196                 exit(1);
 197         }
 198 }
 199 /***********************************************************************
 200 double VecCalc::findMaxDiff(vector<double> vec1, vector<double> vec2){
 201         try {
 202                 double maxDiff = -10000000000.0;
 203                 for(int i = 0; i < vec1.size(); i++)
 204                         if(fabs(vec1.at(i)-vec2.at(i)) > maxDiff)
 205                                 maxDiff = fabs(vec1.at(i)-vec2.at(i));
 206                 return maxDiff;
 207         }
 208         catch(exception& e) {
 209                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 210                 exit(1);
 211         }
 212         catch(...) {
 213                 cout << "An unknown error has occurred in the VecCalc class function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 214                 exit(1);
 215         }
 216 }
 217 /***********************************************************************
 218 double VecCalc::findDStat(vector<double> vec1, vector<double> vec2, double num){
 219         try {
 220                 double mDiff = findMaxDiff(vec1, vec2);
 221                 return (mDiff+.5)/num;
 222         }
 223         catch(exception& e) {
 224                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 225                 exit(1);
 226         }
 227         catch(...) {
 228                 cout << "An unknown error has occurred in the VecCalc class function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 229                 exit(1);
 230         }
 231 }
 232 /***********************************************************************
 233 vector<int> VecCalc::findQuartiles(vector<double> vec){
 234         try {
 235                 vector<int> quartiles;
 236                 double max = vec.at(vec.size()-1);
 237                 double r1 = max/4;
 238                 double r2 = max*3/4;
 239                 bool r1found = false;
 240                 bool r2found = false;
 241                 for(int i = 0; i < vec.size(); i++)
 242                 {
 243                         if(vec.at(i) > r1 && !r1found)
 244                         {
 245                                 quartiles.push_back(i);
 246                                 r1found = true;
 247                         }
 248                         if(vec.at(i) > r2 && !r2found)
 249                         {
 250                                 quartiles.push_back(i);
 251                                 i = vec.size();
 252                         }
 253                 }
 254                 return quartiles;
 255         }
 256         catch(exception& e) {
 257                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 258                 exit(1);
 259         }
 260         catch(...) {
 261                 cout << "An unknown error has occurred in the VecCalc class function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 262                 exit(1);
 263         }
 264 }
 265 /***********************************************************************
 266 vector<double> VecCalc::add(vector<double> vec, double x){
 267         try {
 268                 vector<double> newVec;
 269                 for(int i = 0; i < vec.size(); i++)
 270                         newVec.push_back(vec.at(i)+x);
 271                 return newVec;
 272         }
 273         catch(exception& e) {
 274                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 275                 exit(1);
 276         }
 277         catch(...) {
 278                 cout << "An unknown error has occurred in the VecCalc class function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 279                 exit(1);
 280         }
 281 }
 282 /***********************************************************************
 283 vector<double> VecCalc::multiply(vector<double> vec, double x){
 284         try {
 285                 vector<double> newVec;
 286                 for(int i = 0; i < vec.size(); i++)
 287                         newVec.push_back(vec.at(i)*x);
 288                 return newVec;
 289         }
 290         catch(exception& e) {
 291                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 292                 exit(1);
 293         }
 294         catch(...) {
 295                 cout << "An unknown error has occurred in the VecCalc class function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 296                 exit(1);
 297         }
 298 }
 299 /***********************************************************************
 300 vector<double> VecCalc::power(vector<double> vec, double p){
 301         try {
 302                 vector<double> newVec;
 303                 for(int i = 0; i < vec.size(); i++)
 304                         newVec.push_back(pow(vec.at(i), p));
 305                 return newVec;
 306         }
 307         catch(exception& e) {
 308                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 309                 exit(1);
 310         }
 311         catch(...) {
 312                 cout << "An unknown error has occurred in the VecCalc class function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 313                 exit(1);
 314         }
 315 }
 316 /***********************************************************************
 317 vector<double> VecCalc::addVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
 318 {       try {
 319                 vector<double> newVec;
 320                 for(int i = 0; i < vec1.size(); i++)
 321                         newVec.push_back(vec1.at(i)+vec2.at(i));
 322                 return newVec;
 323         }
 324         catch(exception& e) {
 325                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 326                 exit(1);
 327         }
 328         catch(...) {
 329                 cout << "An unknown error has occurred in the VecCalc class function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 330                 exit(1);
 331         }
 332 }
 333 /***********************************************************************
 334 vector<double> VecCalc::multVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
 335 {       try {
 336                 vector<double> newVec;
 337                 for(int i = 0; i < vec1.size(); i++)
 338                         newVec.push_back(vec1.at(i)*vec2.at(i));
 339                 return newVec;
 340         }
 341         catch(exception& e) {
 342                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 343                 exit(1);
 344         }
 345         catch(...) {
 346                 cout << "An unknown error has occurred in the VecCalc class function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 347                 exit(1);
 348         }
 349 }
 350 /***********************************************************************
 351 vector<double> VecCalc::remDup(vector<double> vec){
 352         try {
 353                 vector<double> newVec;
 354                 newVec.push_back(vec.at(0));
 355                 for(int i = 1; i < vec.size(); i++)
 356                         if(vec.at(i) != vec.at(i-1))
 357                                 newVec.push_back(vec.at(i));
 358                 return newVec;
 359         }
 360         catch(exception& e) {
 361                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 362                 exit(1);
 363         }
 364         catch(...) {
 365                 cout << "An unknown error has occurred in the VecCalc class function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 366                 exit(1);
 367         }
 368 }
 369 /***********************************************************************
 370 vector<double> VecCalc::genCVec(vector<double> vec1){
 371         try {
 372                 vector<double> vec2;
 373                 double sum = 0;
 374                 for(int i = 0; i < vec1.size(); i++)
 375                 {
 376                         sum += vec1.at(i);
 377                         vec2.push_back(sum);
 378                 }
 379                 return vec2;
 380         }
 381         catch(exception& e) {
 382                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 383                 exit(1);
 384         }
 385         catch(...) {
 386                 cout << "An unknown error has occurred in the VecCalc class function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 387                 exit(1);
 388         }
 389 }
 390 /***********************************************************************
 391 vector<double> VecCalc::genRelVec(vector<double> vec){
 392         try {
 393                 vector<double> newVec;
 394                 double sum = sumElements(vec);
 395                 for(int i = 0; i < vec.size(); i++)
 396                         newVec.push_back(vec.at(i)/sum);
 397                 return newVec;
 398         }
 399         catch(exception& e) {
 400                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 401                 exit(1);
 402         }
 403         catch(...) {
 404                 cout << "An unknown error has occurred in the VecCalc class function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 405                 exit(1);
 406         }
 407 }
 408 /***********************************************************************
 409 vector<double> VecCalc::genDiffVec(vector<double> vec1, vector<double> vec2){
 410         try {
 411                 vector<double> newVec;
 412                 for(int i = 0; i < vec1.size(); i++)
 413                         newVec.push_back(vec1.at(i)-vec2.at(i));
 414                 return newVec;
 415         }
 416         catch(exception& e) {
 417                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 418                 exit(1);
 419         }
 420         catch(...) {
 421                 cout << "An unknown error has occurred in the VecCalc class function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 422                 exit(1);
 423         }
 424 }
 425 /***********************************************************************
 426 vector<double> VecCalc::genCSVec(vector<double> vec){
 427         try {
 428                 vector<double> newVec;
 429                 double curInd = vec.at(vec.size()-1);
 430                 double sumSpec = 0;
 431                 double cSumSpec = 0;
 432                 for(int i = vec.size()-1; i >= 0; i--)
 433                 {
 434                         if(vec.at(i) == curInd)
 435                                 sumSpec++;
 436                         else
 437                         {
 438                                 cSumSpec += sumSpec;
 439                                 newVec.push_back(cSumSpec);
 440                                 sumSpec = 1;
 441                                 curInd = vec.at(i);
 442                         }
 443                         if(i == 0)
 444                                 newVec.push_back(cSumSpec + sumSpec);
 445                 }
 446                 return newVec;
 447         }
 448         catch(exception& e) {
 449                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 450                 exit(1);
 451         }
 452         catch(...) {
 453                 cout << "An unknown error has occurred in the VecCalc class function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 454                 exit(1);
 455         }
 456 }
 457 /***********************************************************************
 458 vector<double> VecCalc::genTotVec(vector<vector<double> > vec){
 459         try {
 460                 vector<double> newVec = vec.at(0);
 461                 for(int i = 1; i < vec.size(); i++)
 462                         newVec = addVecs(newVec, vec.at(i));
 463                 return newVec;
 464         }
 465         catch(exception& e) {
 466                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 467                 exit(1);
 468         }
 469         catch(...) {
 470                 cout << "An unknown error has occurred in the VecCalc class function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 471                 exit(1);
 472         }
 473 }
 474 /***********************************************************************
 475 vector<double> VecCalc::quicksort(vector<double> vec){
 476         try {
 477                 sort(vec.rbegin(), vec.rend());
 478                 return vec;
 479         }
 480         catch(exception& e) {
 481                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function quicksort. 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 quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 486                 exit(1);
 487         }
 488 }
 489 /***********************************************************************
 490 vector<vector<double> > VecCalc::gen2DVec(vector<double> vec, int div, int type){
 491         try {
 492                 vector<vector<double> > vec2D;
 493                 int gap = vec.size()/div;
 494                 for(int i = 0; i < div; i++)
 495                 {
 496                         vector<double> inVec;
 497                         for(int j = 0; j < gap; j++)
 498                                 if(type == 0)
 499                                         inVec.push_back(vec.at(j + i*gap)); //Rows
 500                                 else
 501                                         inVec.push_back(vec.at(i + j*div)); //Columns
 502                         vec2D.push_back(inVec);
 503                 }
 504                 return vec2D;
 505         }
 506         catch(exception& e) {
 507                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 508                 exit(1);
 509         }
 510         catch(...) {
 511                 cout << "An unknown error has occurred in the VecCalc class function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 512                 exit(1);
 513         }
 514 }
 515 /***********************************************************************
 516 vector<string> VecCalc::getSData(char name[]){
 517         try {
 518                 vector<string> vec;
 519                 string data;
 520                 ifstream file(name);
 521                 if(file.is_open())
 522                 {
 523                         while(file >> data)
 524                                 vec.push_back(data);
 525                         file.close();
 526                 }
 527                 return vec;
 528         }
 529         catch(exception& e) {
 530                 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 531                 exit(1);
 532         }
 533         catch(...) {
 534                 cout << "An unknown error has occurred in the VecCalc class function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 535                 exit(1);
 536         }
 537 }
 538 /***********************************************************************/
 539
 540 /***********************************************************************
 541 void BrokenStick::doBStick(vector<double> vec)//vec = The data vector.
 542 {       try {
 543                 VecCalc vecCalc;
 544                 vec = vecCalc.quicksort(vec);
 545                 double numInd = vecCalc.sumElements(vec);
 546                 double numSpec = vecCalc.numNZ(vec);
 547                 vector<double> cObsVec = vecCalc.genCVec(vec);
 548                 vector<double> cExpVec;
 549                 vec = vecCalc.power(vec, -1);
 550                 double sumExp = 0;
 551                 for(int i = 0; i < numSpec; i++)
 552                 {
 553                         double n = 0;
 554                         for(int j = i; j < numSpec; j++)
 555                                 n += 1/(double)(j+1);
 556
 557                         sumExp += numInd/numSpec*n;
 558                         cExpVec.push_back(sumExp);
 559                 }
 560
 561                 //Statistical analysis
 562                 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
 563                 double DStatistic = maxDiff/numInd;
 564
 565                 cout << "D-Statistic = " << DStatistic << "\n";
 566                 if(vec.size() > 30)
 567                         cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
 568         }
 569         catch(exception& e) {
 570                 cout << "Standard Error: " << e.what() << " has occurred in the BrokenStick class Function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 571                 exit(1);
 572         }
 573         catch(...) {
 574                 cout << "An unknown error has occurred in the BrokenStick class function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 575                 exit(1);
 576         }
 577 }
 578
 579 /***********************************************************************
 580 double kEq(double k, double spec)
 581 {
 582         return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec));
 583 }
 584 /***********************************************************************/
 585 /*void GeometricSeries::doGeomTest(vector<double> vec)//vec = the data vector
 586 {       try {
 587                 VecCalc vecCalc;
 588                 vec = vecCalc.quicksort(vec);//Sorts vec
 589                 double numInd = vecCalc.sumElements(vec);//numInd = The total number of individuals in the data set.
 590                 double numSpec = vecCalc.numNZ(vec);//numSpec = the number of nonzero elements in the data set.
 591                 double min = -1.0*vecCalc.findMax(vecCalc.multiply(vec, -1));
 592                 double k = .5;
 593                 double step = .49999;
 594                 while(fabs(min - numInd*kEq(k, numSpec)) > .0001)//This uses a binary search to find the value of k.
 595                 {
 596                         if(numInd*kEq(k, numSpec) > min)
 597                                 k += step;
 598                         else
 599                                 k -= step;
 600                         step /= 2;
 601                 }
 602                 cout << "k = " << k << "\n";
 603                 double cK = 1/(1-pow(1-k, numSpec));
 604
 605                 vector<double> cObsVec = vecCalc.genCVec(vec);
 606                 vector<double> cExpVec;
 607                 double sumExp = 0;
 608                 for(int i = 0; i < vec.size(); i++)
 609                 {
 610                         sumExp += numInd*cK*k*pow(1-k, i);
 611                         cExpVec.push_back(sumExp);
 612                 }
 613                 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
 614                 double DStatistic = maxDiff/numInd;
 615
 616                 //Statistical Analysis
 617                 cout << "D-Statistic = " << DStatistic << "\n";
 618                 if(vec.size() > 30)
 619                         cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
 620         }
 621         catch(exception& e) {
 622                 cout << "Standard Error: " << e.what() << " has occurred in the GeometricSeries class Function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 623                 exit(1);
 624         }
 625         catch(...) {
 626                 cout << "An unknown error has occurred in the GeometricSeries class function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 627                 exit(1);
 628         }
 629 }*/
 630
 631 /***********************************************************************
 632 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.
 633 {       try {
 634                 VecCalc vecCalc;
 635                 SSBPDiversityIndices ssbp;
 636                 double rows = vec.size()/cols;
 637                 vector<vector<double> > species = vecCalc.gen2DVec(vec, rows, 0);//converts vec into a 2 dimensional vector
 638                 vector<double> sumRows;
 639                 vector<double> pVals;
 640
 641                 for(int i = 0; i < rows; i++)
 642                         sumRows.push_back(vecCalc.sumElements(species.at(i)));
 643                 double st = 1/ssbp.getSimp(sumRows);//The observed estimate using the Simpson Index. Can be changed to use other indexes of diversity.
 644                 for(int i = 0; i < cols; i++)
 645                 {
 646                         vector<double> newVec;
 647                         for(int j = 0; j < rows; j++)
 648                                 newVec.push_back(vecCalc.sumElements(species.at(j))-species.at(j).at(i));
 649                         pVals.push_back(cols*st-((cols-1)/ssbp.getSimp(newVec)));
 650                 }
 651
 652                 double mean = vecCalc.mean(pVals);
 653                 double stErr = vecCalc.stError(pVals);//stErr = standard error
 654                 TDTable table;
 655                 double confidence = .95;
 656                 double confLimit;
 657                 cout << "Confidence Level (.8, .9, .95, .98, .99, .998, .999): ";
 658                 cin >> confidence;
 659                 double confLevels[] = {.80,.90,.95,.98,.99,.998,.999};
 660                 for(int i = 0; i < 7; i++)
 661                         if(confidence == confLevels[i])
 662                         {
 663                                 confLimit = table.getConfLimit(cols-2,i);
 664                                 i = 7;
 665                         }
 666
 667                 //Statistical Analysis
 668                 cout << "Lower limit = " << mean - confLimit*stErr << "\n";
 669                 cout << "Upper limit = " << mean + confLimit*stErr << "\n";
 670                 cout << "Observed estimate = " << st << "\n\n";
 671                 cout << "Jackknifed estimate = " << mean << "\n";
 672         }
 673         catch(exception& e) {
 674                 cout << "Standard Error: " << e.what() << " has occurred in the Jackknifing class Function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 675                 exit(1);
 676         }
 677         catch(...) {
 678                 cout << "An unknown error has occurred in the Jackknifing class function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 679                 exit(1);
 680         }
 681 }
 682 /***********************************************************************
 683 void KS2SampleTest::doKSTest(vector<double> abun1, vector<double> abun2)//abun1 = 1st vector of abundances, abun2 = 2nd vector of abundances
 684 {       try {
 685                 VecCalc vecCalc;
 686                 double numNZ1 = vecCalc.numNZ(abun1);//Number of nonzero elements in abun1
 687                 double numNZ2 = vecCalc.numNZ(abun2);//Number of nonzero elements in abun2
 688                 abun1 = vecCalc.quicksort(abun1);//Sorts abun1
 689                 abun2 = vecCalc.quicksort(abun2);//Sorts abun2
 690                 abun1 = vecCalc.genRelVec(abun1);//Generates the relative vector for abun1
 691                 abun2 = vecCalc.genRelVec(abun2);//Generates the relative vector for abun2
 692
 693                 abun1 = vecCalc.genCVec(abun1);//Generates the cumulative vector for abun1
 694                 abun2 = vecCalc.genCVec(abun2);//Generates the cumulative vector for abun2
 695
 696                 double maxDiff = vecCalc.findMaxDiff(abun1, abun2);
 697                 double DStatistic = maxDiff*numNZ1*numNZ2;
 698
 699                 cout << "Null Hypothesis = There is no difference.\n";
 700                 cout << "D-Statistic = " << DStatistic << "\n";
 701
 702                 double a = pow((numNZ1 + numNZ2)/(numNZ1*numNZ2),.5);
 703                 double pVal = exp(-2*pow(maxDiff/a,2));
 704
 705                 if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25.
 706                         cout << "P-Value = " << pVal << "\n\n";
 707                 else
 708                 {
 709                         cout << "D-Statistic must be higher than the critical value to reject the null hypothesis.\n" ;
 710                         cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n";
 711                         cout << "95% Confidence Critical Value = " << 1.36*a*numNZ1*numNZ2 << "\n";
 712                 }
 713         }
 714         catch(exception& e) {
 715                 cout << "Standard Error: " << e.what() << " has occurred in the KS2SampleTest class Function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 716                 exit(1);
 717         }
 718         catch(...) {
 719                 cout << "An unknown error has occurred in the KS2SampleTest class function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 720                 exit(1);
 721         }
 722 }
 723
 724 /***********************************************************************
 725
 726 void QStatistic::doQStat(vector<double> vec)//vec = The data vector.
 727 {       try {
 728                 VecCalc vecCalc;
 729                 vector<double> cVec = vecCalc.genCSVec(vec);
 730                 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.
 731                 vector<double> nDupVec = vecCalc.remDup(vec);//nDupVec only contains one of every unique element in cVec.
 732                 double Q;
 733                 if(q.at(0) != 0)//The case if neither quartile is 0 or 1
 734                         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)));
 735                 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.
 736                         Q = (.5*cVec.at(0) + .5*(cVec.at(1)-cVec.at(0)))/log(nDupVec.at(nDupVec.size()-2)/nDupVec.at(nDupVec.size()-1));
 737                 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.
 738                         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));
 739
 740                 cout << "Q = " << Q << "\n";
 741         }
 742         catch(exception& e) {
 743                 cout << "Standard Error: " << e.what() << " has occurred in the QStatistic class Function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 744                 exit(1);
 745         }
 746         catch(...) {
 747                 cout << "An unknown error has occurred in the QStatistic class function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 748                 exit(1);
 749         }
 750 }
 751 /***********************************************************************
 752 double SSBPDiversityIndices::getShan(vector<double> vec)//vec = The data vector.
 753 {       try {
 754                 VecCalc vecCalc;
 755                 double nz = vecCalc.numNZ(vec);
 756                 double nSum = vecCalc.sumElements(vec);
 757                 double H = 0;
 758                 for(int i = 0; i < nz; i++)
 759                         H += vec.at(i)/nSum*log(vec.at(i)/nSum);
 760                 H *= -1;
 761                 return H;
 762         }
 763         catch(exception& e) {
 764                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 765                 exit(1);
 766         }
 767         catch(...) {
 768                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 769                 exit(1);
 770         }
 771 }
 772 /***********************************************************************
 773 double SSBPDiversityIndices::getSimp(vector<double> vec)//vec = The data vector.
 774 {       try {
 775                 VecCalc vecCalc;
 776                 double nSum = vecCalc.sumElements(vec);
 777                 double D = 0;
 778                 for(int j = 0; j < vec.size(); j++)
 779                         D += vec.at(j)*(vec.at(j)-1)/(nSum*(nSum-1));
 780                 return D;
 781         }
 782         catch(exception& e) {
 783                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 784                 exit(1);
 785         }
 786         catch(...) {
 787                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 788                 exit(1);
 789         }
 790 }
 791 /***********************************************************************
 792 double SSBPDiversityIndices::getBP(vector<double> vec)//vec = The data vector.
 793 {       try {
 794                 VecCalc vecCalc;
 795                 double nSum = vecCalc.sumElements(vec);
 796                 return vecCalc.findMax(vec)/nSum;
 797         }
 798         catch(exception& e) {
 799                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 800                 exit(1);
 801         }
 802         catch(...) {
 803                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 804                 exit(1);
 805         }
 806 }
 807 /***********************************************************************
 808 void SSBPDiversityIndices::doSSBP(vector<double> vec)//vec = The data vector.
 809 {       try {
 810                 VecCalc vecCalc;
 811                 double nz = vecCalc.numNZ(vec);
 812
 813                 //Shannon index
 814                 double H = getShan(vec);
 815                 cout << "H = " << H << "\n";
 816                 cout << "Eveness = " << H/log(nz) << "\n\n";
 817
 818                 //Simpson index
 819                 double D = getSimp(vec);
 820                 cout << "D diversity = " << 1/D << "\n";
 821                 cout << "Eveness = " << 1/D/nz << "\n\n";
 822
 823                 //Berger-Parker index
 824                 double BP = getBP(vec);
 825                 cout << "BP index = " << BP << "\n\n";
 826         }
 827         catch(exception& e) {
 828                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 829                 exit(1);
 830         }
 831         catch(...) {
 832                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 833                 exit(1);
 834         }
 835 }
 836 /***********************************************************************/
 837 double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom
 838 {    try {
 839                 //Found on http://www.vgtu.lt/leidiniai/elektroniniai/Probability.pdf/Table%203.pdf
 840
 841                 //Confidence Level        .90    .95     .975     .99    .995     .999    .9995
 842                 double values[30][7] = {{3.078, 6.314,  12.706, 31.821, 63.656, 318.289, 636.578},
 843                                                         {1.886, 2.920,  4.303,  6.965,  9.925,  22.328, 31.600},
 844                                                             {1.638,     2.353,  3.182,  4.541,  5.841,  10.214, 12.924},
 845                                                             {1.533,     2.132,  2.776,  3.747,  4.604,  7.173,  8.610},
 846                                                                 {1.476, 2.015,  2.571,  3.365,  4.032,  5.894,  6.869},
 847                                                                 {1.440, 1.943,  2.447,  3.143,  3.707,  5.208,  5.959},
 848                                                                 {1.415, 1.895,  2.365,  2.998,  3.499,  4.785,  5.408},
 849                                                                 {1.397, 1.860,  2.306,  2.896,  3.355,  4.501,  5.041},
 850                                                                 {1.383, 1.833,  2.262,  2.821,  3.250,  4.297,  4.781},
 851                                                                 {1.372, 1.812,  2.228,  2.764,  3.169,  4.144,  4.587},
 852                                                                 {1.363, 1.796,  2.201,  2.718,  3.106,  4.025,  4.437},
 853                                                                 {1.356, 1.782,  2.179,  2.681,  3.055,  3.930,  4.318},
 854                                                                 {1.350, 1.771,  2.160,  2.650,  3.012,  3.852,  4.221},
 855                                                                 {1.345, 1.761,  2.145,  2.624,  2.977,  3.787,  4.140},
 856                                                                 {1.341, 1.753,  2.131,  2.602,  2.947,  3.733,  4.073},
 857                                                                 {1.337, 1.746,  2.120,  2.583,  2.921,  3.686,  4.015},
 858                                                                 {1.333, 1.740,  2.110,  2.567,  2.898,  3.646,  3.965},
 859                                                                 {1.330, 1.734,  2.101,  2.552,  2.878,  3.610,  3.922},
 860                                                                 {1.328, 1.729,  2.093,  2.539,  2.861,  3.579,  3.883},
 861                                                                 {1.325, 1.725,  2.086,  2.528,  2.845,  3.552,  3.850},
 862                                                                 {1.323, 1.721,  2.080,  2.518,  2.831,  3.527,  3.819},
 863                                                                 {1.321, 1.717,  2.074,  2.508,  2.819,  3.505,  3.792},
 864                                                                 {1.319, 1.714,  2.069,  2.500,  2.807,  3.485,  3.768},
 865                                                                 {1.318, 1.711,  2.064,  2.492,  2.797,  3.467,  3.745},
 866                                                                 {1.316, 1.708,  2.060,  2.485,  2.787,  3.450,  3.725},
 867                                                                 {1.315, 1.706,  2.056,  2.479,  2.779,  3.435,  3.707},
 868                                                                 {1.314, 1.703,  2.052,  2.473,  2.771,  3.421,  3.689},
 869                                                                 {1.313, 1.701,  2.048,  2.467,  2.763,  3.408,  3.674},
 870                                                                 {1.311, 1.699,  2.045,  2.462,  2.756,  3.396,  3.660},
 871                                                                 {1.310, 1.697,  2.042,  2.457,  2.750,  3.385,  3.646}};
 872
 873                 return values[row][col];
 874         }
 875         catch(exception& e) {
 876                 errorOut(e, "TDTable", "getConfLimit");
 877                 exit(1);
 878         }
 879 }
 880
 881 /***********************************************************************
 882 double KOSTable::getConfLimit(int index) //Table of Critical values for N = 4-20 for One-Sample Kolmogorov-Smirnov Test
 883 {    try {
 884                 //Confidence Level = .05
 885                 //Sample size = 4-20.
 886                 //Found on http://www.utdallas.edu/~herve/MolinAbdi1998-LillieforsTechReport.pdf
 887
 888                 double values[21] = {.9011,
 889                                                          .6372,
 890                                                          .5203,
 891                                                          .4506,
 892                                                          0.3754,
 893                                                  0.3427,
 894                                                          0.3245,
 895                                                  0.3041,
 896                                                  0.2875,
 897                                                  0.2744,
 898                                                  0.2616,
 899                                                  0.2506,
 900                                                  0.2426,
 901                                                  0.2337,
 902                                                  0.2257,
 903                                                  0.2196,
 904                                                  0.2128,
 905                                                  0.2071,
 906                                                  0.2018,
 907                                                  0.1965,
 908                                                  0.1920,
 909                                                          };
 910                 return values[index];
 911         }
 912         catch(exception& e) {
 913                 cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 914                 exit(1);
 915         }
 916         catch(...) {
 917                 cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 918                 exit(1);
 919         }
 920 }
 921
 922 /***********************************************************************
 923 void TrunLN::doTrunLN(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
 924 {
 925         VecCalc vecCalc;
 926         double numSpec = vecCalc.sumElements(specVec); //numSpec = The total number of species
 927         double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec,specVec)); //numInd = The total number of individuals
 928
 929         double obsMean = 0;
 930         for(int i = 0; i < indVec.size(); i++)
 931                 obsMean += log10(indVec.at(i));
 932         obsMean /= numSpec; //obsMean = observed mean of the individuals vector
 933         cout << "obsMean = " << obsMean << "\n";
 934         double variance = 0;
 935         for(int t = 0; t < indVec.size(); t++)
 936                 variance += pow(log10(indVec.at(t))-obsMean,2)/numSpec;
 937
 938          double rO = 0;
 939          for(int k = 0; k < indVec.size(); k++)
 940                 rO += log10(indVec.at(k));
 941          rO /= indVec.size();
 942          double veilLine = .5;//The desired veil line.
 943          double auxFunc = -(obsMean-rO)/(obsMean-log10(veilLine));
 944          double uX = obsMean-auxFunc*(obsMean-log10(veilLine));
 945          double vX = variance + auxFunc*pow(obsMean-log10(veilLine),2);
 946          double z = (log10(veilLine)-uX)/pow(vX, .5);
 947          double p = .5*(erf(z)+1);
 948          double specRichness = numSpec/(1-p);
 949
 950          double numUnseen = .5*(erf((log10(.5)-uX)/pow(vX,.5))+1)*specRichness;
 951
 952
 953          vector<double> cExp;
 954          for(int i = 1; i < 8; i++)
 955          {
 956                 double a = pow(10, i-1)+.5;
 957                 double b = log10(a);
 958                 double c = (b - uX)/pow(vX,.5);
 959                 double d = .5*(erf(c)+1);
 960                 double numS = d*specRichness;
 961                 double toPush = numS - numUnseen;
 962                 cExp.push_back(toPush);
 963         }
 964         vector<double> cObs;
 965         double sumOct = 0;
 966         for(int i = 0; i < 8; i++)
 967         {
 968                 sumOct = 0;
 969                 for(int r = 0; r < indVec.size(); r++)
 970                 {
 971                         if(indVec.at(r) < pow(10, i-1)+.5)
 972                                 sumOct += specVec.at(r);
 973                         else
 974                         {
 975                                 cObs.push_back(sumOct);
 976                                 sumOct = specVec.at(r);
 977                                 r = indVec.size();
 978                         }
 979                         if(r == indVec.size()-1)
 980                                 cObs.push_back(sumOct);
 981                 }
 982         }
 983
 984         //Statistical Analysis
 985         double d = vecCalc.findDStat(cExp, cObs, numSpec);
 986         cout << "DStat = " << d << "\n";
 987         cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
 988         cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
 989 }
 990 /***********************************************************************/