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