X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=sharedjackknife.cpp;fp=sharedjackknife.cpp;h=0998db8fe3106bbb40c26dc8621d37030fb275c3;hb=0caf3fbabaa3ece404f8ce77f4c883dc5b1bf1dc;hp=0000000000000000000000000000000000000000;hpb=1b73ff67c83892a025e597dabd9df6fe7b58206a;p=mothur.git

diff --git a/sharedjackknife.cpp b/sharedjackknife.cpp
new file mode 100644
index 0000000..0998db8
--- /dev/null
+++ b/sharedjackknife.cpp
@@ -0,0 +1,144 @@
+/*
+ *  sharedjackknife.cpp
+ *  Mothur
+ *
+ *  Created by Thomas Ryabin on 3/30/09.
+ *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
+ *
+ */
+
+#include "sharedjackknife.h"
+
+/***************************************************************************************
+
+***************************************************************************************/
+double SharedJackknife::simpson(vector<int> abunds, double numInd, int numBins){
+	double denom = numInd*(numInd-1);
+	double sum = 0;
+	for(int i = 0; i < numBins; i++)
+		sum += (double)abunds[i]*((double)abunds[i]-1)/denom;
+	
+	return sum;
+}
+
+/*****************************************************************************************/
+
+double* SharedJackknife::jackknife(){		
+	int numBins = groups.at(0)->getNumBins()-1;
+	vector<int> cArray(numBins);
+	for(int i = 0; i < numBins; i++)
+		cArray[i] = 0;
+
+	double numInd = 0;
+	for(int i = 0; i < numGroups; i++)
+		for(int j = 0; j < numBins; j++) {
+			int curAbund = groups.at(i)->get(j+1).abundance;
+			cArray[j] += curAbund;
+			numInd += (double)curAbund;
+		}
+
+	double baseD = 1/simpson(cArray, numInd, numBins);
+	
+	vector<double> pseudoVals(numBins);
+	double jackknifeEstimate = 0;
+	for(int i = 0; i < numGroups; i++) {
+		for(int j = 0; j < numBins-1; j++) {
+			int abundDiff = -groups.at(i)->get(j+1).abundance;
+			if(i > 0)
+				abundDiff += groups.at(i-1)->get(j+1).abundance;
+
+			cArray[j] += abundDiff;
+			numInd += abundDiff;	
+		}
+		
+		double curD = 1/simpson(cArray, numInd, numBins);
+		pseudoVals[i] = (double)numGroups*(baseD - curD) + curD;
+		jackknifeEstimate += pseudoVals[i];
+	}
+	jackknifeEstimate /= (double)numGroups;
+		
+	double variance = 0;
+	for(int i = 0; i < numGroups; i++)
+		variance += pow(pseudoVals[i]-jackknifeEstimate, 2);
+		
+	variance /= (double)numGroups*((double)numGroups-1);
+	double stErr = sqrt(variance);
+	TDTable table;
+	double confLimit = 0;
+	if(numGroups <= 30)
+		confLimit = table.getConfLimit(numGroups-1, 1);
+	else
+		confLimit = 1.645;
+	
+	confLimit *= stErr;
+	
+	double* rdata = new double[3];
+	rdata[0] = baseD;
+	rdata[1] = jackknifeEstimate - confLimit;
+	rdata[2] = jackknifeEstimate + confLimit;
+	
+	return rdata;
+}
+
+/************************************************************************************************
+************************************************************************************************/
+
+EstOutput SharedJackknife::getValues(vector<SharedRAbundVector*> vectorShared){ //Fix this for collect, mistake was that it was made with summary in mind.
+	try {
+		SharedRAbundVector* shared1 = vectorShared[0];
+		SharedRAbundVector* shared2 = vectorShared[1];
+		if(numGroups == -1) {
+			numGroups = m->getNumGroups();
+		}
+
+		if(callCount == numGroups*(numGroups-1)/2) {
+			currentCallDone = true;
+			callCount = 0;
+		}
+		callCount++;
+
+		if(currentCallDone) {
+			groups.clear();	
+			currentCallDone = false;
+		}
+		
+		if(groups.size() != numGroups) {	
+			if(groups.size() == 0)
+				groups.push_back(shared1);
+			groups.push_back(shared2);
+		}
+		
+		if(groups.size() == numGroups && callCount < numGroups) {
+			data.resize(3,0);
+
+			double* rdata = jackknife();
+			data[0] = rdata[0];
+			data[1] = rdata[1];
+			data[2] = rdata[2];
+		
+			if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
+			if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
+			if (isnan(data[2]) || isinf(data[0])) { data[2] = 0; }
+			
+			return data;
+		}
+		
+		data.resize(3,0);
+		data[0] = 0;
+		data[1] = 0;
+		data[2] = 0;
+		
+		if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
+		if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
+		if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
+		return data;	
+	}
+		
+	catch(exception& e) {
+		m->errorOut(e, "SharedJackknife", "getValues");
+		exit(1);
+	}
+}
+
+/***********************************************************************/
+