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
 750 void QStatistic::doQStat(vector<double> vec)//vec = The data vector.
 751 {       try {
 752                 VecCalc vecCalc;
 753                 vector<double> cVec = vecCalc.genCSVec(vec);
 754                 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.
 755                 vector<double> nDupVec = vecCalc.remDup(vec);//nDupVec only contains one of every unique element in cVec.
 756                 double Q;
 757                 if(q.at(0) != 0)//The case if neither quartile is 0 or 1
 758                         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)));
 759                 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.
 760                         Q = (.5*cVec.at(0) + .5*(cVec.at(1)-cVec.at(0)))/log(nDupVec.at(nDupVec.size()-2)/nDupVec.at(nDupVec.size()-1));
 761                 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.
 762                         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));
 763
 764                 cout << "Q = " << Q << "\n";
 765         }
 766         catch(exception& e) {
 767                 cout << "Standard Error: " << e.what() << " has occurred in the QStatistic class Function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 768                 exit(1);
 769         }
 770         catch(...) {
 771                 cout << "An unknown error has occurred in the QStatistic class function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 772                 exit(1);
 773         }
 774 }
 775 /***********************************************************************/
 776 double SSBPDiversityIndices::getShan(vector<double> vec)//vec = The data vector.
 777 {       try {
 778                 VecCalc vecCalc;
 779                 double nz = vecCalc.numNZ(vec);
 780                 double nSum = vecCalc.sumElements(vec);
 781                 double H = 0;
 782                 for(int i = 0; i < nz; i++)
 783                         H += vec.at(i)/nSum*log(vec.at(i)/nSum);
 784                 H *= -1;
 785                 return H;
 786         }
 787         catch(exception& e) {
 788                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 789                 exit(1);
 790         }
 791         catch(...) {
 792                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 793                 exit(1);
 794         }
 795 }
 796 /***********************************************************************/
 797 double SSBPDiversityIndices::getSimp(vector<double> vec)//vec = The data vector.
 798 {       try {
 799                 VecCalc vecCalc;
 800                 double nSum = vecCalc.sumElements(vec);
 801                 double D = 0;
 802                 for(int j = 0; j < vec.size(); j++)
 803                         D += vec.at(j)*(vec.at(j)-1)/(nSum*(nSum-1));
 804                 return D;
 805         }
 806         catch(exception& e) {
 807                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 808                 exit(1);
 809         }
 810         catch(...) {
 811                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 812                 exit(1);
 813         }
 814 }
 815 /***********************************************************************/
 816 double SSBPDiversityIndices::getBP(vector<double> vec)//vec = The data vector.
 817 {       try {
 818                 VecCalc vecCalc;
 819                 double nSum = vecCalc.sumElements(vec);
 820                 return vecCalc.findMax(vec)/nSum;
 821         }
 822         catch(exception& e) {
 823                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 824                 exit(1);
 825         }
 826         catch(...) {
 827                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 828                 exit(1);
 829         }
 830 }
 831 /***********************************************************************/
 832 void SSBPDiversityIndices::doSSBP(vector<double> vec)//vec = The data vector.
 833 {       try {
 834                 VecCalc vecCalc;
 835                 double nz = vecCalc.numNZ(vec);
 836
 837                 //Shannon index
 838                 double H = getShan(vec);
 839                 cout << "H = " << H << "\n";
 840                 cout << "Eveness = " << H/log(nz) << "\n\n";
 841
 842                 //Simpson index
 843                 double D = getSimp(vec);
 844                 cout << "D diversity = " << 1/D << "\n";
 845                 cout << "Eveness = " << 1/D/nz << "\n\n";
 846
 847                 //Berger-Parker index
 848                 double BP = getBP(vec);
 849                 cout << "BP index = " << BP << "\n\n";
 850         }
 851         catch(exception& e) {
 852                 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 853                 exit(1);
 854         }
 855         catch(...) {
 856                 cout << "An unknown error has occurred in the SSBPDiversityIndices class function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 857                 exit(1);
 858         }
 859 }
 860 /***********************************************************************/
 861 double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom
 862 {    try {
 863                 //Found on http://www.vgtu.lt/leidiniai/elektroniniai/Probability.pdf/Table%203.pdf
 864
 865                 //Confidence Level        .90    .95     .975     .99    .995     .999    .9995
 866                 double values[33][7] = {{3.078, 6.314,  12.706, 31.821, 63.656, 318.289, 636.578},
 867                                                         {1.886, 2.920,  4.303,  6.965,  9.925,  22.328, 31.600},
 868                                                             {1.638,     2.353,  3.182,  4.541,  5.841,  10.214, 12.924},
 869                                                             {1.533,     2.132,  2.776,  3.747,  4.604,  7.173,  8.610},
 870                                                                 {1.476, 2.015,  2.571,  3.365,  4.032,  5.894,  6.869},
 871                                                                 {1.440, 1.943,  2.447,  3.143,  3.707,  5.208,  5.959},
 872                                                                 {1.415, 1.895,  2.365,  2.998,  3.499,  4.785,  5.408},
 873                                                                 {1.397, 1.860,  2.306,  2.896,  3.355,  4.501,  5.041},
 874                                                                 {1.383, 1.833,  2.262,  2.821,  3.250,  4.297,  4.781},
 875                                                                 {1.372, 1.812,  2.228,  2.764,  3.169,  4.144,  4.587},
 876                                                                 {1.363, 1.796,  2.201,  2.718,  3.106,  4.025,  4.437},
 877                                                                 {1.356, 1.782,  2.179,  2.681,  3.055,  3.930,  4.318},
 878                                                                 {1.350, 1.771,  2.160,  2.650,  3.012,  3.852,  4.221},
 879                                                                 {1.345, 1.761,  2.145,  2.624,  2.977,  3.787,  4.140},
 880                                                                 {1.341, 1.753,  2.131,  2.602,  2.947,  3.733,  4.073},
 881                                                                 {1.337, 1.746,  2.120,  2.583,  2.921,  3.686,  4.015},
 882                                                                 {1.333, 1.740,  2.110,  2.567,  2.898,  3.646,  3.965},
 883                                                                 {1.330, 1.734,  2.101,  2.552,  2.878,  3.610,  3.922},
 884                                                                 {1.328, 1.729,  2.093,  2.539,  2.861,  3.579,  3.883},
 885                                                                 {1.325, 1.725,  2.086,  2.528,  2.845,  3.552,  3.850},
 886                                                                 {1.323, 1.721,  2.080,  2.518,  2.831,  3.527,  3.819},
 887                                                                 {1.321, 1.717,  2.074,  2.508,  2.819,  3.505,  3.792},
 888                                                                 {1.319, 1.714,  2.069,  2.500,  2.807,  3.485,  3.768},
 889                                                                 {1.318, 1.711,  2.064,  2.492,  2.797,  3.467,  3.745},
 890                                                                 {1.316, 1.708,  2.060,  2.485,  2.787,  3.450,  3.725},
 891                                                                 {1.315, 1.706,  2.056,  2.479,  2.779,  3.435,  3.707},
 892                                                                 {1.314, 1.703,  2.052,  2.473,  2.771,  3.421,  3.689},
 893                                                                 {1.313, 1.701,  2.048,  2.467,  2.763,  3.408,  3.674},
 894                                                                 {1.311, 1.699,  2.045,  2.462,  2.756,  3.396,  3.660},
 895                                                                 {1.310, 1.697,  2.042,  2.457,  2.750,  3.385,  3.646},
 896                                                                 {1.296, 1.671,  2.000,  2.390,  2.660,  3.232,  3.460},
 897                                                                 {1.289, 1.658,  1.980,  2.358,  2.617,  3.160,  3.373},
 898                                                                 {1.282, 1.645,  1.960,  2.326,  2.576,  3.091,  3.291}};
 899                 return values[row][col];
 900         }
 901         catch(exception& e) {
 902                 cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 903                 exit(1);
 904         }
 905         catch(...) {
 906                 cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 907                 exit(1);
 908         }
 909 }
 910
 911 /***********************************************************************/
 912 double KOSTable::getConfLimit(int index) //Table of Critical values for N = 4-20 for One-Sample Kolmogorov-Smirnov Test
 913 {    try {
 914                 //Confidence Level = .05
 915                 //Sample size = 4-20.
 916                 //Found on http://www.utdallas.edu/~herve/MolinAbdi1998-LillieforsTechReport.pdf
 917
 918                 double values[21] = {.9011,
 919                                                          .6372,
 920                                                          .5203,
 921                                                          .4506,
 922                                                          0.3754,
 923                                                  0.3427,
 924                                                          0.3245,
 925                                                  0.3041,
 926                                                  0.2875,
 927                                                  0.2744,
 928                                                  0.2616,
 929                                                  0.2506,
 930                                                  0.2426,
 931                                                  0.2337,
 932                                                  0.2257,
 933                                                  0.2196,
 934                                                  0.2128,
 935                                                  0.2071,
 936                                                  0.2018,
 937                                                  0.1965,
 938                                                  0.1920,
 939                                                          };
 940                 return values[index];
 941         }
 942         catch(exception& e) {
 943                 cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 944                 exit(1);
 945         }
 946         catch(...) {
 947                 cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 948                 exit(1);
 949         }
 950 }
 951
 952 /***********************************************************************
 953 void TrunLN::doTrunLN(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
 954 {
 955         VecCalc vecCalc;
 956         double numSpec = vecCalc.sumElements(specVec); //numSpec = The total number of species
 957         double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec,specVec)); //numInd = The total number of individuals
 958
 959         double obsMean = 0;
 960         for(int i = 0; i < indVec.size(); i++)
 961                 obsMean += log10(indVec.at(i));
 962         obsMean /= numSpec; //obsMean = observed mean of the individuals vector
 963         cout << "obsMean = " << obsMean << "\n";
 964         double variance = 0;
 965         for(int t = 0; t < indVec.size(); t++)
 966                 variance += pow(log10(indVec.at(t))-obsMean,2)/numSpec;
 967
 968          double rO = 0;
 969          for(int k = 0; k < indVec.size(); k++)
 970                 rO += log10(indVec.at(k));
 971          rO /= indVec.size();
 972          double veilLine = .5;//The desired veil line.
 973          double auxFunc = -(obsMean-rO)/(obsMean-log10(veilLine));
 974          double uX = obsMean-auxFunc*(obsMean-log10(veilLine));
 975          double vX = variance + auxFunc*pow(obsMean-log10(veilLine),2);
 976          double z = (log10(veilLine)-uX)/pow(vX, .5);
 977          double p = .5*(erf(z)+1);
 978          double specRichness = numSpec/(1-p);
 979
 980          double numUnseen = .5*(erf((log10(.5)-uX)/pow(vX,.5))+1)*specRichness;
 981
 982
 983          vector<double> cExp;
 984          for(int i = 1; i < 8; i++)
 985          {
 986                 double a = pow(10, i-1)+.5;
 987                 double b = log10(a);
 988                 double c = (b - uX)/pow(vX,.5);
 989                 double d = .5*(erf(c)+1);
 990                 double numS = d*specRichness;
 991                 double toPush = numS - numUnseen;
 992                 cExp.push_back(toPush);
 993         }
 994         vector<double> cObs;
 995         double sumOct = 0;
 996         for(int i = 0; i < 8; i++)
 997         {
 998                 sumOct = 0;
 999                 for(int r = 0; r < indVec.size(); r++)
1000                 {
1001                         if(indVec.at(r) < pow(10, i-1)+.5)
1002                                 sumOct += specVec.at(r);
1003                         else
1004                         {
1005                                 cObs.push_back(sumOct);
1006                                 sumOct = specVec.at(r);
1007                                 r = indVec.size();
1008                         }
1009                         if(r == indVec.size()-1)
1010                                 cObs.push_back(sumOct);
1011                 }
1012         }
1013
1014         //Statistical Analysis
1015         double d = vecCalc.findDStat(cExp, cObs, numSpec);
1016         cout << "DStat = " << d << "\n";
1017         cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
1018         cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
1019 }
1020 /***********************************************************************/