5 * Created by Sarah Westcott on 11/18/08.
6 * Copyright 2008 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "calculator.h"
12 /***********************************************************************/
13 void VecCalc::printElements(vector<double> vec){
15 for(int i = 0; i < vec.size(); i++)
16 cout << vec.at(i) << " ";
20 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
24 cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
29 /***********************************************************************/
31 void VecCalc::printElements(vector<string> vec){
33 for(int i = 0; i < vec.size(); i++)
34 cout << vec.at(i) << " ";
38 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
42 cout << "An unknown error has occurred in the VecCalc class function printElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
46 /***********************************************************************/
47 int VecCalc::findString(vector<string> vec, string str){
49 for(int i = 0; i < vec.size(); i++)
50 if(str.compare(vec.at(i)) == 0)
55 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
59 cout << "An unknown error has occurred in the VecCalc class function findString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
63 /***********************************************************************/
64 double VecCalc::mean(vector<double> vec){
65 return sumElements(vec)/vec.size();
67 /***********************************************************************/
68 double VecCalc::stError(vector<double> 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);
77 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
81 cout << "An unknown error has occurred in the VecCalc class function stError. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
85 /***********************************************************************/
86 double VecCalc::sumElements(vector<double> vec){
89 for(int i = 0; i < vec.size(); i++)
94 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
98 cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
102 /***********************************************************************/
103 double VecCalc::sumElements(vector<double> vec, int index){
106 for(int i = index; i < vec.size(); i++)
110 catch(exception& e) {
111 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
115 cout << "An unknown error has occurred in the VecCalc class function sumElements. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
119 /***********************************************************************/
120 double VecCalc::findMax(vector<double> vec){
122 double max = -1000000.0;
123 for(int i = 0; i < vec.size(); i++)
128 catch(exception& e) {
129 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
133 cout << "An unknown error has occurred in the VecCalc class function findMax. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
137 /***********************************************************************/
138 double VecCalc::numNZ(vector<double> vec){
141 for(int i = 0; i < vec.size(); i++)
146 catch(exception& e) {
147 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
151 cout << "An unknown error has occurred in the VecCalc class function numNZ. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
155 /***********************************************************************/
156 double VecCalc::numPos(vector<double> vec){
159 for(int i = 0 ; i < vec.size(); i++)
164 catch(exception& e) {
165 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
169 cout << "An unknown error has occurred in the VecCalc class function numPos. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
173 /***********************************************************************/
174 double VecCalc::findMaxDiff(vector<double> vec1, vector<double> vec2){
176 double maxDiff = -10000000000.0;
177 for(int i = 0; i < vec1.size(); i++)
178 if(fabs(vec1.at(i)-vec2.at(i)) > maxDiff)
179 maxDiff = fabs(vec1.at(i)-vec2.at(i));
182 catch(exception& e) {
183 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
187 cout << "An unknown error has occurred in the VecCalc class function findMaxDiff. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
191 /***********************************************************************/
192 double VecCalc::findDStat(vector<double> vec1, vector<double> vec2, double num){
194 double mDiff = findMaxDiff(vec1, vec2);
195 return (mDiff+.5)/num;
197 catch(exception& e) {
198 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
202 cout << "An unknown error has occurred in the VecCalc class function findDStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
206 /***********************************************************************/
207 vector<int> VecCalc::findQuartiles(vector<double> vec){
209 vector<int> quartiles;
210 double max = vec.at(vec.size()-1);
213 bool r1found = false;
214 bool r2found = false;
215 for(int i = 0; i < vec.size(); i++)
217 if(vec.at(i) > r1 && !r1found)
219 quartiles.push_back(i);
222 if(vec.at(i) > r2 && !r2found)
224 quartiles.push_back(i);
230 catch(exception& e) {
231 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
235 cout << "An unknown error has occurred in the VecCalc class function findQuartiles. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
239 /***********************************************************************/
240 vector<double> VecCalc::add(vector<double> vec, double x){
242 vector<double> newVec;
243 for(int i = 0; i < vec.size(); i++)
244 newVec.push_back(vec.at(i)+x);
247 catch(exception& e) {
248 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
252 cout << "An unknown error has occurred in the VecCalc class function add. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
256 /***********************************************************************/
257 vector<double> VecCalc::multiply(vector<double> vec, double x){
259 vector<double> newVec;
260 for(int i = 0; i < vec.size(); i++)
261 newVec.push_back(vec.at(i)*x);
264 catch(exception& e) {
265 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
269 cout << "An unknown error has occurred in the VecCalc class function multiply. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
273 /***********************************************************************/
274 vector<double> VecCalc::power(vector<double> vec, double p){
276 vector<double> newVec;
277 for(int i = 0; i < vec.size(); i++)
278 newVec.push_back(pow(vec.at(i), p));
281 catch(exception& e) {
282 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
286 cout << "An unknown error has occurred in the VecCalc class function power. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
290 /***********************************************************************/
291 vector<double> VecCalc::addVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
293 vector<double> newVec;
294 for(int i = 0; i < vec1.size(); i++)
295 newVec.push_back(vec1.at(i)+vec2.at(i));
298 catch(exception& e) {
299 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
303 cout << "An unknown error has occurred in the VecCalc class function addVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
307 /***********************************************************************/
308 vector<double> VecCalc::multVecs(vector<double> vec1, vector<double> vec2) //Vectors must be the same size.
310 vector<double> newVec;
311 for(int i = 0; i < vec1.size(); i++)
312 newVec.push_back(vec1.at(i)*vec2.at(i));
315 catch(exception& e) {
316 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
320 cout << "An unknown error has occurred in the VecCalc class function multVecs. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
324 /***********************************************************************/
325 vector<double> VecCalc::remDup(vector<double> vec){
327 vector<double> newVec;
328 newVec.push_back(vec.at(0));
329 for(int i = 1; i < vec.size(); i++)
330 if(vec.at(i) != vec.at(i-1))
331 newVec.push_back(vec.at(i));
334 catch(exception& e) {
335 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
339 cout << "An unknown error has occurred in the VecCalc class function remDup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
343 /***********************************************************************/
344 vector<double> VecCalc::genCVec(vector<double> vec1){
348 for(int i = 0; i < vec1.size(); i++)
355 catch(exception& e) {
356 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
360 cout << "An unknown error has occurred in the VecCalc class function genCVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
364 /***********************************************************************/
365 vector<double> VecCalc::genRelVec(vector<double> vec){
367 vector<double> newVec;
368 double sum = sumElements(vec);
369 for(int i = 0; i < vec.size(); i++)
370 newVec.push_back(vec.at(i)/sum);
373 catch(exception& e) {
374 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
378 cout << "An unknown error has occurred in the VecCalc class function genRelVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
382 /***********************************************************************/
383 vector<double> VecCalc::genDiffVec(vector<double> vec1, vector<double> vec2){
385 vector<double> newVec;
386 for(int i = 0; i < vec1.size(); i++)
387 newVec.push_back(vec1.at(i)-vec2.at(i));
390 catch(exception& e) {
391 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
395 cout << "An unknown error has occurred in the VecCalc class function genDiffVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
399 /***********************************************************************/
400 vector<double> VecCalc::genCSVec(vector<double> vec){
402 vector<double> newVec;
403 double curInd = vec.at(vec.size()-1);
406 for(int i = vec.size()-1; i >= 0; i--)
408 if(vec.at(i) == curInd)
413 newVec.push_back(cSumSpec);
418 newVec.push_back(cSumSpec + sumSpec);
422 catch(exception& e) {
423 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
427 cout << "An unknown error has occurred in the VecCalc class function genCSVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
431 /***********************************************************************/
432 vector<double> VecCalc::genTotVec(vector<vector<double> > vec){
434 vector<double> newVec = vec.at(0);
435 for(int i = 1; i < vec.size(); i++)
436 newVec = addVecs(newVec, vec.at(i));
439 catch(exception& e) {
440 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
444 cout << "An unknown error has occurred in the VecCalc class function genTotVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
448 /***********************************************************************/
449 vector<double> VecCalc::quicksort(vector<double> vec){
451 sort(vec.rbegin(), vec.rend());
454 catch(exception& e) {
455 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
459 cout << "An unknown error has occurred in the VecCalc class function quicksort. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
463 /***********************************************************************/
464 vector<vector<double> > VecCalc::gen2DVec(vector<double> vec, int div, int type){
466 vector<vector<double> > vec2D;
467 int gap = vec.size()/div;
468 for(int i = 0; i < div; i++)
470 vector<double> inVec;
471 for(int j = 0; j < gap; j++)
473 inVec.push_back(vec.at(j + i*gap)); //Rows
475 inVec.push_back(vec.at(i + j*div)); //Columns
476 vec2D.push_back(inVec);
480 catch(exception& e) {
481 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
485 cout << "An unknown error has occurred in the VecCalc class function gen2DVec. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
489 /***********************************************************************/
490 vector<string> VecCalc::getSData(char name[]){
503 catch(exception& e) {
504 cout << "Standard Error: " << e.what() << " has occurred in the VecCalc class Function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
508 cout << "An unknown error has occurred in the VecCalc class function getSData. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
512 /***********************************************************************/
513 double BDiversity::getWhitt(vector<double> vec1, vector<double> vec2){
516 double numSpec = vecCalc.numNZ(vecCalc.addVecs(vec1,vec2));
517 return 2*numSpec/(vecCalc.numNZ(vec1)+vecCalc.numNZ(vec2))-1;
519 catch(exception& e) {
520 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getWhitt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
524 cout << "An unknown error has occurred in the BDiversity class function getWhitt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
528 /***********************************************************************/
529 double BDiversity::getMS(vector<double> vec1, vector<double> vec2){
532 double a = vecCalc.numNZ(vecCalc.multVecs(vec1, vec2));
533 double b = vecCalc.numPos(vecCalc.addVecs(vecCalc.addVecs(vec1, vecCalc.multiply(vec2, -1)), vecCalc.multiply(vecCalc.multVecs(vec1, vec2), -1)));
534 double c = vecCalc.numPos(vecCalc.addVecs(vecCalc.addVecs(vec2, vecCalc.multiply(vec1, -1)), vecCalc.multiply(vecCalc.multVecs(vec2, vec1), -1)));
537 catch(exception& e) {
538 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getMS. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
542 cout << "An unknown error has occurred in the BDiversity class function getMS. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
546 /***********************************************************************/
547 double BDiversity::getSor(vector<double> vec1, vector<double> vec2){
552 for(int i = 0; i < vec1.size(); i++)
556 if(vec1.at(i) >= vec2.at(i))
561 return 2*sum/(asum+bsum);
563 catch(exception& e) {
564 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getSor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
568 cout << "An unknown error has occurred in the BDiversity class function getSor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
572 /***********************************************************************/
573 double BDiversity::getMor(vector<double> vec1, vector<double> vec2){
580 for(int i = 0; i < vec1.size(); i++)
582 sum += vec1.at(i)*vec2.at(i);
583 asum += pow(vec1.at(i),2);
584 bsum += pow(vec2.at(i),2);
588 return 2*sum/((asum/pow(sum1,2)+bsum/pow(sum2,2))*(sum1*sum2));
590 catch(exception& e) {
591 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function getMor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
595 cout << "An unknown error has occurred in the BDiversity class function getMor. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
599 /***********************************************************************/
600 void BDiversity::printD(vector<vector<double> > columns, int type){
603 for(int i = 0; i < columns.size(); i++)
604 cout << "Column:" << i << " ";
606 for(int i = 0; i < columns.size(); i++)
608 cout << "Column " << i << ":";
609 for(int j = 0; j < columns.size(); j++)
615 B = getWhitt(columns.at(i), columns.at(j));
617 B = 1-getMS(columns.at(i), columns.at(j));
619 B = 1-getSor(columns.at(i), columns.at(j));
621 B = 1-getMor(columns.at(i), columns.at(j));
631 catch(exception& e) {
632 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function printD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
636 cout << "An unknown error has occurred in the BDiversity class function printD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
640 /***********************************************************************/
641 void BDiversity::doBD(vector<double> vec, double cols)//vec = The data vector if the data table was read left-to-right, top-to-bottom. cols = the number of columns in the data table.
644 vector<vector<double> > columns = vecCalc.gen2DVec(vec,cols,1);
646 cout.setf(ios_base::fixed, ios_base::floatfield);
647 cout.precision(6);//This formats the data so the tables look pretty.
649 //Whittaker's measure Bw (presence/absence data)
650 cout << "Whittaker's measure Bw (presence/absence data):\n";
653 for(int i = 0; i < cols; i++)
654 sum += vecCalc.numNZ(columns.at(i));
655 double meanRichness = sum/cols;
656 vector<double> totVec = vecCalc.genTotVec(columns);
657 double totRichness = vecCalc.numNZ(totVec);
658 cout << "\nOverall B Diversity = " << totRichness/meanRichness << "\n\n\n";//The overall B Diversity
660 //Marczewski-Steinhaus(MS) distance(Jaccard index) (presence/abscence data)
661 cout << "Marczewski-Steinhaus(MS) distance(Jaccard index) (presence/abscence data):\n";
664 for(int i = 0; i < cols; i++)
665 for(int j = 0; j < cols; j++)
667 sum += vecCalc.numNZ(columns.at(i))+vecCalc.numNZ(columns.at(j)) - 2*vecCalc.numNZ(vecCalc.multVecs(columns.at(i), columns.at(j)));
668 cout << "\nOverall B Diversity = " << sum/cols << "\n\n\n";//The overall B Diversity
670 //Sorensen quantitative index (abundance data)
671 cout << "Sorensen quantitative index (abundance data):\n";
675 //Sorensen quantitative index (abundance data)
676 cout << "Morisita-Horn index (abundance data):\n";
679 catch(exception& e) {
680 cout << "Standard Error: " << e.what() << " has occurred in the BDiversity class Function doBD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
684 cout << "An unknown error has occurred in the BDiversity class function doBD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
689 /***********************************************************************/
690 void BrokenStick::doBStick(vector<double> vec)//vec = The data vector.
693 vec = vecCalc.quicksort(vec);
694 double numInd = vecCalc.sumElements(vec);
695 double numSpec = vecCalc.numNZ(vec);
696 vector<double> cObsVec = vecCalc.genCVec(vec);
697 vector<double> cExpVec;
698 vec = vecCalc.power(vec, -1);
700 for(int i = 0; i < numSpec; i++)
703 for(int j = i; j < numSpec; j++)
704 n += 1/(double)(j+1);
706 sumExp += numInd/numSpec*n;
707 cExpVec.push_back(sumExp);
710 //Statistical analysis
711 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
712 double DStatistic = maxDiff/numInd;
714 cout << "D-Statistic = " << DStatistic << "\n";
716 cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
718 catch(exception& e) {
719 cout << "Standard Error: " << e.what() << " has occurred in the BrokenStick class Function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
723 cout << "An unknown error has occurred in the BrokenStick class function doBStick. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
728 /***********************************************************************/
729 double kEq(double k, double spec)
731 return k/(1-k)*pow(1-k, spec)/(1-pow(1-k, spec));
733 /***********************************************************************/
734 void GeometricSeries::doGeomTest(vector<double> vec)//vec = the data vector
737 vec = vecCalc.quicksort(vec);//Sorts vec
738 double numInd = vecCalc.sumElements(vec);//numInd = The total number of individuals in the data set.
739 double numSpec = vecCalc.numNZ(vec);//numSpec = the number of nonzero elements in the data set.
740 double min = -1.0*vecCalc.findMax(vecCalc.multiply(vec, -1));
742 double step = .49999;
743 while(fabs(min - numInd*kEq(k, numSpec)) > .0001)//This uses a binary search to find the value of k.
745 if(numInd*kEq(k, numSpec) > min)
751 cout << "k = " << k << "\n";
752 double cK = 1/(1-pow(1-k, numSpec));
754 vector<double> cObsVec = vecCalc.genCVec(vec);
755 vector<double> cExpVec;
757 for(int i = 0; i < vec.size(); i++)
759 sumExp += numInd*cK*k*pow(1-k, i);
760 cExpVec.push_back(sumExp);
762 double maxDiff = vecCalc.findMaxDiff(cObsVec, cExpVec);
763 double DStatistic = maxDiff/numInd;
765 //Statistical Analysis
766 cout << "D-Statistic = " << DStatistic << "\n";
768 cout << "Critical value for 95% confidence interval = " << .886/sqrt(vec.size()) << "\n";
770 catch(exception& e) {
771 cout << "Standard Error: " << e.what() << " has occurred in the GeometricSeries class Function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
775 cout << "An unknown error has occurred in the GeometricSeries class function doGeomTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
780 /***********************************************************************/
781 void Jackknifing::doJK(vector<double> vec, double cols)//vec = the data vector if the data table was read left-to-right, top-to-bottom. cols = The number of columns in the data table.
784 SSBPDiversityIndices ssbp;
785 double rows = vec.size()/cols;
786 vector<vector<double> > species = vecCalc.gen2DVec(vec, rows, 0);//converts vec into a 2 dimensional vector
787 vector<double> sumRows;
788 vector<double> pVals;
790 for(int i = 0; i < rows; i++)
791 sumRows.push_back(vecCalc.sumElements(species.at(i)));
792 double st = 1/ssbp.getSimp(sumRows);//The observed estimate using the Simpson Index. Can be changed to use other indexes of diversity.
793 for(int i = 0; i < cols; i++)
795 vector<double> newVec;
796 for(int j = 0; j < rows; j++)
797 newVec.push_back(vecCalc.sumElements(species.at(j))-species.at(j).at(i));
798 pVals.push_back(cols*st-((cols-1)/ssbp.getSimp(newVec)));
801 double mean = vecCalc.mean(pVals);
802 double stErr = vecCalc.stError(pVals);//stErr = standard error
804 double confidence = .95;
806 cout << "Confidence Level (.8, .9, .95, .98, .99, .998, .999): ";
808 double confLevels[] = {.80,.90,.95,.98,.99,.998,.999};
809 for(int i = 0; i < 7; i++)
810 if(confidence == confLevels[i])
812 confLimit = table.getConfLimit(cols-2,i);
816 //Statistical Analysis
817 cout << "Lower limit = " << mean - confLimit*stErr << "\n";
818 cout << "Upper limit = " << mean + confLimit*stErr << "\n";
819 cout << "Observed estimate = " << st << "\n\n";
820 cout << "Jackknifed estimate = " << mean << "\n";
822 catch(exception& e) {
823 cout << "Standard Error: " << e.what() << " has occurred in the Jackknifing class Function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
827 cout << "An unknown error has occurred in the Jackknifing class function doJK. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
831 /***********************************************************************/
832 void KS2SampleTest::doKSTest(vector<double> abun1, vector<double> abun2)//abun1 = 1st vector of abundances, abun2 = 2nd vector of abundances
835 double numNZ1 = vecCalc.numNZ(abun1);//Number of nonzero elements in abun1
836 double numNZ2 = vecCalc.numNZ(abun2);//Number of nonzero elements in abun2
837 abun1 = vecCalc.quicksort(abun1);//Sorts abun1
838 abun2 = vecCalc.quicksort(abun2);//Sorts abun2
839 abun1 = vecCalc.genRelVec(abun1);//Generates the relative vector for abun1
840 abun2 = vecCalc.genRelVec(abun2);//Generates the relative vector for abun2
842 abun1 = vecCalc.genCVec(abun1);//Generates the cumulative vector for abun1
843 abun2 = vecCalc.genCVec(abun2);//Generates the cumulative vector for abun2
845 double maxDiff = vecCalc.findMaxDiff(abun1, abun2);
846 double DStatistic = maxDiff*numNZ1*numNZ2;
848 cout << "Null Hypothesis = There is no difference.\n";
849 cout << "D-Statistic = " << DStatistic << "\n";
851 double a = pow((numNZ1 + numNZ2)/(numNZ1*numNZ2),.5);
852 double pVal = exp(-2*pow(maxDiff/a,2));
854 if(numNZ1 > 25 && numNZ2 > 25) //If the sample sizes are both bigger than 25.
855 cout << "P-Value = " << pVal << "\n\n";
858 cout << "D-Statistic must be higher than the critical value to reject the null hypothesis.\n" ;
859 cout << "90% Confidence Critical Value = " << 1.22*a*numNZ1*numNZ2 << "\n";
860 cout << "95% Confidence Critical Value = " << 1.36*a*numNZ1*numNZ2 << "\n";
863 catch(exception& e) {
864 cout << "Standard Error: " << e.what() << " has occurred in the KS2SampleTest class Function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
868 cout << "An unknown error has occurred in the KS2SampleTest class function doKSTest. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
873 /***********************************************************************/
874 double logS(double num)
876 return -(1-num)*log(1-num)/num;
879 /***********************************************************************/
880 void LogSD::doLogSD(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
883 double numSpec = vecCalc.sumElements(specVec);//numSpec = The total number of species
884 cout << "number of species = " << numSpec << "\n";
885 double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec, specVec));
886 double snRatio = numSpec/numInd;
888 double step = .4999999999;
889 while(fabs(snRatio - logS(x)) > .00001) //This uses a binary search to find the value of x.
891 if(logS(x) > snRatio)
897 double alpha = numInd*(1-x)/x;
900 cout << "Number of individuals:"; //Ask the user for the number of individuals.
902 double spec = alpha*pow(x, ind)/ind;
903 cout << "Number of species expected = " << spec << "\n" << "X value = " << x << "\n" << "Alpha value= " << alpha << "\n";//Outputs the number of species expected with the given number of individuals.
905 vector<double> obsSpec;
906 vector<double> cObsSpec;
907 vector<double> expSpec;
908 vector<double> cExpSpec;
909 vector<double> cDiff;
911 // Generates the cumulative observed species vector.
913 double octSumObs = 0;
914 for(int y = 0; y < specVec.size(); y++)
916 if(indVec.at(y) - .5 < pow(2.0, oct))
917 octSumObs += specVec.at(y);
920 obsSpec.push_back(octSumObs);
921 octSumObs = specVec.at(y);
924 if(y == specVec.size()-1)
925 obsSpec.push_back(octSumObs);
927 cObsSpec = vecCalc.genCVec(obsSpec);
928 cObsSpec = vecCalc.add(cObsSpec,-.5);
930 // Generates the cumulative expected species vector.
932 double octSumExp = 0;
933 for(int g = 1; g <= indVec.at(indVec.size()-1); g++)
935 if(g - .5 < pow(2.0, oct))
936 octSumExp += alpha*pow(x,g)/(g);
939 expSpec.push_back(octSumExp);
940 octSumExp = alpha*pow(x,g)/(g);
943 if(g == indVec.at(indVec.size()-1))
944 expSpec.push_back(octSumExp);
946 cExpSpec = vecCalc.genCVec(expSpec);
948 // Statistical Analysis
949 double dTStat = vecCalc.findDStat(cObsSpec, cExpSpec, numSpec);
950 cout << "D Test Statistic = " << dTStat << "\n";
951 cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
952 cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
954 catch(exception& e) {
955 cout << "Standard Error: " << e.what() << " has occurred in the LogSD class Function doLogSD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
959 cout << "An unknown error has occurred in the LogSD class function doLogSD. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
963 /***********************************************************************/
964 void QStatistic::doQStat(vector<double> vec)//vec = The data vector.
967 vector<double> cVec = vecCalc.genCSVec(vec);
968 vector<int> q = vecCalc.findQuartiles(cVec);//q.at(0) is the index of the first quartile; q.at(1) is the index of the second quartile.
969 vector<double> nDupVec = vecCalc.remDup(vec);//nDupVec only contains one of every unique element in cVec.
971 if(q.at(0) != 0)//The case if neither quartile is 0 or 1
972 Q = (.5*(cVec.at(q.at(0))-cVec.at(q.at(0)-1)) + (cVec.at(q.at(1)-1)-cVec.at(q.at(0))) + .5*(cVec.at(q.at(1))-cVec.at(q.at(1)-1)))/log(nDupVec.at(nDupVec.size()-1-q.at(1))/nDupVec.at(nDupVec.size()-1-q.at(0)));
973 else if(q.at(0) == 0 && (q.at(1) == 0 || q.at(1) == 1))//The case if the quartiles are both at index 0 or one is at 0 and the other at 1.
974 Q = (.5*cVec.at(0) + .5*(cVec.at(1)-cVec.at(0)))/log(nDupVec.at(nDupVec.size()-2)/nDupVec.at(nDupVec.size()-1));
975 else if(q.at(0) == 0 && q.at(1) > 1) //The case if the lower quartile is at index 0 and upper quartile index is above index 1.
976 Q = (.5*cVec.at(0) + (cVec.at(q.at(1)-1)-cVec.at(q.at(0))) + .5*(cVec.at(q.at(1))-cVec.at(q.at(1)-1)))/log(nDupVec.at(nDupVec.size()-1-q.at(1))/nDupVec.at(nDupVec.size()-1));
978 cout << "Q = " << Q << "\n";
980 catch(exception& e) {
981 cout << "Standard Error: " << e.what() << " has occurred in the QStatistic class Function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
985 cout << "An unknown error has occurred in the QStatistic class function doQStat. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
989 /***********************************************************************/
990 double SSBPDiversityIndices::getShan(vector<double> vec)//vec = The data vector.
993 double nz = vecCalc.numNZ(vec);
994 double nSum = vecCalc.sumElements(vec);
996 for(int i = 0; i < nz; i++)
997 H += vec.at(i)/nSum*log(vec.at(i)/nSum);
1001 catch(exception& e) {
1002 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1006 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getShan. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1010 /***********************************************************************/
1011 double SSBPDiversityIndices::getSimp(vector<double> vec)//vec = The data vector.
1014 double nSum = vecCalc.sumElements(vec);
1016 for(int j = 0; j < vec.size(); j++)
1017 D += vec.at(j)*(vec.at(j)-1)/(nSum*(nSum-1));
1020 catch(exception& e) {
1021 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1025 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getSimp. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1029 /***********************************************************************/
1030 double SSBPDiversityIndices::getBP(vector<double> vec)//vec = The data vector.
1033 double nSum = vecCalc.sumElements(vec);
1034 return vecCalc.findMax(vec)/nSum;
1036 catch(exception& e) {
1037 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1041 cout << "An unknown error has occurred in the SSBPDiversityIndices class function getBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1045 /***********************************************************************/
1046 void SSBPDiversityIndices::doSSBP(vector<double> vec)//vec = The data vector.
1049 double nz = vecCalc.numNZ(vec);
1052 double H = getShan(vec);
1053 cout << "H = " << H << "\n";
1054 cout << "Eveness = " << H/log(nz) << "\n\n";
1057 double D = getSimp(vec);
1058 cout << "D diversity = " << 1/D << "\n";
1059 cout << "Eveness = " << 1/D/nz << "\n\n";
1061 //Berger-Parker index
1062 double BP = getBP(vec);
1063 cout << "BP index = " << BP << "\n\n";
1065 catch(exception& e) {
1066 cout << "Standard Error: " << e.what() << " has occurred in the SSBPDiversityIndices class Function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1070 cout << "An unknown error has occurred in the SSBPDiversityIndices class function doSSBP. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1074 /***********************************************************************/
1075 double TDTable::getConfLimit(int row, int col) //Rows are the degrees of freedom
1077 //Confidence Level .80 .90 .95 .98 .99 .998 .999
1078 double values[33][7] = {{3.078, 6.314, 12.706, 31.821, 63.656, 318.289, 636.578},
1079 {1.886, 2.920, 4.303, 6.965, 9.925, 22.328, 31.600},
1080 {1.638, 2.353, 3.182, 4.541, 5.841, 10.214, 12.924},
1081 {1.533, 2.132, 2.776, 3.747, 4.604, 7.173, 8.610},
1082 {1.476, 2.015, 2.571, 3.365, 4.032, 5.894, 6.869},
1083 {1.440, 1.943, 2.447, 3.143, 3.707, 5.208, 5.959},
1084 {1.415, 1.895, 2.365, 2.998, 3.499, 4.785, 5.408},
1085 {1.397, 1.860, 2.306, 2.896, 3.355, 4.501, 5.041},
1086 {1.383, 1.833, 2.262, 2.821, 3.250, 4.297, 4.781},
1087 {1.372, 1.812, 2.228, 2.764, 3.169, 4.144, 4.587},
1088 {1.363, 1.796, 2.201, 2.718, 3.106, 4.025, 4.437},
1089 {1.356, 1.782, 2.179, 2.681, 3.055, 3.930, 4.318},
1090 {1.350, 1.771, 2.160, 2.650, 3.012, 3.852, 4.221},
1091 {1.345, 1.761, 2.145, 2.624, 2.977, 3.787, 4.140},
1092 {1.341, 1.753, 2.131, 2.602, 2.947, 3.733, 4.073},
1093 {1.337, 1.746, 2.120, 2.583, 2.921, 3.686, 4.015},
1094 {1.333, 1.740, 2.110, 2.567, 2.898, 3.646, 3.965},
1095 {1.330, 1.734, 2.101, 2.552, 2.878, 3.610, 3.922},
1096 {1.328, 1.729, 2.093, 2.539, 2.861, 3.579, 3.883},
1097 {1.325, 1.725, 2.086, 2.528, 2.845, 3.552, 3.850},
1098 {1.323, 1.721, 2.080, 2.518, 2.831, 3.527, 3.819},
1099 {1.321, 1.717, 2.074, 2.508, 2.819, 3.505, 3.792},
1100 {1.319, 1.714, 2.069, 2.500, 2.807, 3.485, 3.768},
1101 {1.318, 1.711, 2.064, 2.492, 2.797, 3.467, 3.745},
1102 {1.316, 1.708, 2.060, 2.485, 2.787, 3.450, 3.725},
1103 {1.315, 1.706, 2.056, 2.479, 2.779, 3.435, 3.707},
1104 {1.314, 1.703, 2.052, 2.473, 2.771, 3.421, 3.689},
1105 {1.313, 1.701, 2.048, 2.467, 2.763, 3.408, 3.674},
1106 {1.311, 1.699, 2.045, 2.462, 2.756, 3.396, 3.660},
1107 {1.310, 1.697, 2.042, 2.457, 2.750, 3.385, 3.646},
1108 {1.296, 1.671, 2.000, 2.390, 2.660, 3.232, 3.460},
1109 {1.289, 1.658, 1.980, 2.358, 2.617, 3.160, 3.373},
1110 {1.282, 1.645, 1.960, 2.326, 2.576, 3.091, 3.291}};
1111 return values[row][col];
1113 catch(exception& e) {
1114 cout << "Standard Error: " << e.what() << " has occurred in the TDTable class Function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1118 cout << "An unknown error has occurred in the TDTable class function getConfLimit. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
1122 /***********************************************************************
1123 void TrunLN::doTrunLN(vector<double> indVec, vector<double> specVec) //indVec = individuals vector, specVec = species vector
1126 double numSpec = vecCalc.sumElements(specVec); //numSpec = The total number of species
1127 double numInd = vecCalc.sumElements(vecCalc.multVecs(indVec,specVec)); //numInd = The total number of individuals
1130 for(int i = 0; i < indVec.size(); i++)
1131 obsMean += log10(indVec.at(i));
1132 obsMean /= numSpec; //obsMean = observed mean of the individuals vector
1133 cout << "obsMean = " << obsMean << "\n";
1134 double variance = 0;
1135 for(int t = 0; t < indVec.size(); t++)
1136 variance += pow(log10(indVec.at(t))-obsMean,2)/numSpec;
1139 for(int k = 0; k < indVec.size(); k++)
1140 rO += log10(indVec.at(k));
1141 rO /= indVec.size();
1142 double veilLine = .5;//The desired veil line.
1143 double auxFunc = -(obsMean-rO)/(obsMean-log10(veilLine));
1144 double uX = obsMean-auxFunc*(obsMean-log10(veilLine));
1145 double vX = variance + auxFunc*pow(obsMean-log10(veilLine),2);
1146 double z = (log10(veilLine)-uX)/pow(vX, .5);
1147 double p = .5*(erf(z)+1);
1148 double specRichness = numSpec/(1-p);
1150 double numUnseen = .5*(erf((log10(.5)-uX)/pow(vX,.5))+1)*specRichness;
1153 vector<double> cExp;
1154 for(int i = 1; i < 8; i++)
1156 double a = pow(10, i-1)+.5;
1157 double b = log10(a);
1158 double c = (b - uX)/pow(vX,.5);
1159 double d = .5*(erf(c)+1);
1160 double numS = d*specRichness;
1161 double toPush = numS - numUnseen;
1162 cExp.push_back(toPush);
1164 vector<double> cObs;
1166 for(int i = 0; i < 8; i++)
1169 for(int r = 0; r < indVec.size(); r++)
1171 if(indVec.at(r) < pow(10, i-1)+.5)
1172 sumOct += specVec.at(r);
1175 cObs.push_back(sumOct);
1176 sumOct = specVec.at(r);
1179 if(r == indVec.size()-1)
1180 cObs.push_back(sumOct);
1184 //Statistical Analysis
1185 double d = vecCalc.findDStat(cExp, cObs, numSpec);
1186 cout << "DStat = " << d << "\n";
1187 cout << ".05 confidence value = " << .89196/sqrt(numSpec) << "\n";
1188 cout << ".01 confidence value = " << 1.0471/sqrt(numSpec) << "\n\n";
1190 /***********************************************************************/