/*
 *  ccode.cpp
 *  Mothur
 *
 *  Created by westcott on 8/24/09.
 *  Copyright 2009 Schloss Lab. All rights reserved.
 *
 */

#include "ccode.h"
#include "ignoregaps.h"
#include "eachgapdist.h"

//***************************************************************************************************************
Ccode::Ccode(string filename, string temp, bool f, string mask, int win, int numW, string o) : Chimera() {  
 try {	
	
	fastafile = filename;  
	outputDir = o; 
	templateFileName = temp;  templateSeqs = readSeqs(temp); 
	setMask(mask);
	filter = f;
	window = win;
	numWanted = numW;
	
	distCalc = new eachGapDist();
	decalc = new DeCalculator();
	
	mapInfo = outputDir + m->getRootName(m->getSimpleName(fastafile)) + "mapinfo";
	
	#ifdef USE_MPI
		
		//char* inFileName = new char[mapInfo.length()];
		//memcpy(inFileName, mapInfo.c_str(), mapInfo.length());
		
		char inFileName[1024];
		strcpy(inFileName, mapInfo.c_str());
		
		int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;

		MPI_File_open(MPI_COMM_WORLD, inFileName, outMode, MPI_INFO_NULL, &outMap);  //comm, filename, mode, info, filepointer
		
		//delete inFileName;

		int pid;
		MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
		
		if (pid == 0) {
			string outString = "Place in masked, filtered and trimmed sequence\tPlace in original alignment\n";
			
			MPI_Status status;
			int length = outString.length();
			char* buf2 = new char[length];
			memcpy(buf2, outString.c_str(), length);
				
			MPI_File_write_shared(outMap, buf2, length, MPI_CHAR, &status);
			delete buf2;
		}
	#else

		ofstream out2;
		m->openOutputFile(mapInfo, out2);
		
		out2 << "Place in masked, filtered and trimmed sequence\tPlace in original alignment" << endl;
		out2.close();
	#endif
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "Ccode");
		exit(1);
	}
}
//***************************************************************************************************************
Ccode::~Ccode() {
	delete distCalc;
	delete decalc;
	
	#ifdef USE_MPI
		MPI_File_close(&outMap);
	#endif
}	
//***************************************************************************************************************
Sequence Ccode::print(ostream& out, ostream& outAcc) {
	try {
		
		ofstream out2;
		m->openOutputFileAppend(mapInfo, out2);
		
		out2 << querySeq->getName() << endl;
		for (it = spotMap.begin(); it!= spotMap.end(); it++) {
			out2 << it->first << '\t' << it->second << endl;
		}
		out2.close();
		out << querySeq->getName() << endl << endl << "Reference sequences used and distance to query:" << endl;
			
		for (int j = 0; j < closest.size(); j++) {
			out << closest[j].seq->getName() << '\t' << closest[j].dist << endl;
		}
		out << endl << endl;
		
		//for each window
		//window mapping info.
		out << "Mapping information: ";
		//you mask and did not filter
		if ((seqMask != "") && (!filter)) { out << "mask and trim."; }
				
		//you filtered and did not mask
		if ((seqMask == "") && (filter)) { out << "filter and trim."; }
				
		//you masked and filtered
		if ((seqMask != "") && (filter)) { out << "mask, filter and trim."; }
			
		out << endl << "Window\tStartPos\tEndPos" << endl;
		it = trim.begin();
		for (int k = 0; k < windows.size()-1; k++) {
			out << k+1 << '\t' << spotMap[windows[k]-it->first] << '\t' << spotMap[windows[k]-it->first+windowSizes] << endl;
		}
			
		out << windows.size() << '\t' << spotMap[windows[windows.size()-1]-it->first] << '\t' << spotMap[it->second-it->first-1] << endl;
		out << endl;
		out << "Window\tAvgQ\t(sdQ)\tAvgR\t(sdR)\tRatio\tAnova" << endl;
		for (int k = 0; k < windows.size(); k++) {
			float ds = averageQuery[k] / averageRef[k]; 
			out << k+1 << '\t' << averageQuery[k] << '\t' << sdQuery[k] << '\t' << averageRef[k] << '\t'<< sdRef[k] << '\t' << ds << '\t' << anova[k] << endl;
		}
		out << endl;
			
		//varRef
		//varQuery
		/* F test for differences among variances.
		* varQuery is expected to be higher or similar than varRef */
		//float fs = varQuery[query] / varRef[query];	/* F-Snedecor, test for differences of variances */
			
		bool results = false;	
					
		//confidence limit, t - Student, anova
		out << "Window\tConfidenceLimit\tt-Student\tAnova" << endl;
			
		for (int k = 0; k < windows.size(); k++) {
			string temp = "";
			if (isChimericConfidence[k]) {  temp += "*\t"; }
			else { temp += "\t"; }
				
			if (isChimericTStudent[k]) {  temp += "*\t"; }
			else { temp += "\t"; }
				
			if (isChimericANOVA[k]) {  temp += "*\t"; }
			else { temp += "\t"; }
			
			out << k+1 << '\t' << temp << endl;
				
			if (temp == "*\t*\t*\t") {  results = true;  }
		}
		out << endl;	
			
		if (results) {
			m->mothurOut(querySeq->getName() + " was found have at least one chimeric window."); m->mothurOutEndLine();
			outAcc << querySeq->getName() << endl;
		}

		//free memory
		for (int i = 0; i < closest.size(); i++) {  delete closest[i].seq;  }

		return *querySeq;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "print");
		exit(1);
	}
}
#ifdef USE_MPI
//***************************************************************************************************************
Sequence Ccode::print(MPI_File& out, MPI_File& outAcc) {
	try {
		
		string outMapString = "";
		
		outMapString += querySeq->getName() + "\n";
		for (it = spotMap.begin(); it!= spotMap.end(); it++) {
			outMapString += toString(it->first)  + "\t" + toString(it->second)  + "\n";
		}
		printMapping(outMapString);
		outMapString = "";
		
		string outString = "";
		string outAccString = "";
		
		outString +=  querySeq->getName() + "\n\nReference sequences used and distance to query:\n";
			
		for (int j = 0; j < closest.size(); j++) {
			outString += closest[j].seq->getName() + "\t" + toString(closest[j].dist) + "\n";
		}
		outString += "\n\nMapping information: ";
		
		//for each window
		//window mapping info.
		//you mask and did not filter
		if ((seqMask != "") && (!filter)) { outString += "mask and trim."; }
				
		//you filtered and did not mask
		if ((seqMask == "") && (filter)) { outString += "filter and trim."; }
				
		//you masked and filtered
		if ((seqMask != "") && (filter)) { outString += "mask, filter and trim."; }
			
		outString += "\nWindow\tStartPos\tEndPos\n";
		it = trim.begin();
		for (int k = 0; k < windows.size()-1; k++) {
			outString += toString(k+1) + "\t" + toString(spotMap[windows[k]-it->first]) + "\t" + toString(spotMap[windows[k]-it->first+windowSizes]) + "\n";
		}
			
		outString += toString(windows.size()) + "\t" + toString(spotMap[windows[windows.size()-1]-it->first]) + "\t" + toString(spotMap[it->second-it->first-1]) + "\n\n";
		
		outString += "Window\tAvgQ\t(sdQ)\tAvgR\t(sdR)\tRatio\tAnova\n";
		for (int k = 0; k < windows.size(); k++) {
			float ds = averageQuery[k] / averageRef[k]; 
			outString += toString(k+1) + "\t" + toString(averageQuery[k]) + "\t" + toString(sdQuery[k]) + "\t" + toString(averageRef[k]) + "\t" + toString(sdRef[k]) + "\t" + toString(ds) + "\t" + toString(anova[k]) + "\n";
		}
			
		//varRef
		//varQuery
		/* F test for differences among variances.
		* varQuery is expected to be higher or similar than varRef */
		//float fs = varQuery[query] / varRef[query];	/* F-Snedecor, test for differences of variances */
			
		bool results = false;	
					
		//confidence limit, t - Student, anova
		outString += "\nWindow\tConfidenceLimit\tt-Student\tAnova\n";
			
		for (int k = 0; k < windows.size(); k++) {
			string temp = "";
			if (isChimericConfidence[k]) {  temp += "*\t"; }
			else { temp += "\t"; }
				
			if (isChimericTStudent[k]) {  temp += "*\t"; }
			else { temp += "\t"; }
				
			if (isChimericANOVA[k]) {  temp += "*\t"; }
			else { temp += "\t"; }
			
			outString += toString(k+1) + "\t" + temp + "\n";
				
			if (temp == "*\t*\t*\t") {  results = true;  }
		}
		outString += "\n";	
		
		MPI_Status status;
		int length = outString.length();
		char* buf2 = new char[length];
		memcpy(buf2, outString.c_str(), length);
				
		MPI_File_write_shared(out, buf2, length, MPI_CHAR, &status);
		delete buf2;

		if (results) {
			m->mothurOut(querySeq->getName() + " was found have at least one chimeric window."); m->mothurOutEndLine();
			outAccString += querySeq->getName() + "\n";
			
			MPI_Status statusAcc;
			length = outAccString.length();
			char* buf = new char[length];
			memcpy(buf, outAccString.c_str(), length);
				
			MPI_File_write_shared(outAcc, buf, length, MPI_CHAR, &statusAcc);
			delete buf;
		}

		//free memory
		for (int i = 0; i < closest.size(); i++) {  delete closest[i].seq;  }

		return *querySeq;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "print");
		exit(1);
	}
}
//***************************************************************************************************************
int Ccode::printMapping(string& output) {
	try {
			MPI_Status status;
			int length = output.length();
			char* buf = new char[length];
			memcpy(buf, output.c_str(), length);
				
			MPI_File_write_shared(outMap, buf, length, MPI_CHAR, &status);
			delete buf;
			
			return 0;

	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "printMapping");
		exit(1);
	}
}
#endif
//***************************************************************************************************************
int Ccode::getChimeras(Sequence* query) {
	try {
	
		closest.clear();
		refCombo = 0;
		sumRef.clear(); 
		varRef.clear(); 
		varQuery.clear(); 
		sdRef.clear(); 
		sdQuery.clear();     
		sumQuery.clear();
		sumSquaredRef.clear(); 
		sumSquaredQuery.clear(); 
		averageRef.clear();
		averageQuery.clear();
		anova.clear();
		isChimericConfidence.clear();
		isChimericTStudent.clear();
		isChimericANOVA.clear();
		trim.clear();
		spotMap.clear();
		windowSizes = window;
		windows.clear();

	
		querySeq = query;
		
		//find closest matches to query
		closest = findClosest(query, numWanted);
		
		if (m->control_pressed) {  return 0;  }
		
		//initialize spotMap
		for (int i = 0; i < query->getAligned().length(); i++) {	spotMap[i] = i;		}
	
		//mask sequences if the user wants to 
		if (seqMask != "") {
			decalc->setMask(seqMask);
			
			decalc->runMask(query);
			
			//mask closest
			for (int i = 0; i < closest.size(); i++) {	decalc->runMask(closest[i].seq);	}
			
			spotMap = decalc->getMaskMap();
		}
		
		if (filter) {
			vector<Sequence*> temp;
			for (int i = 0; i < closest.size(); i++) { temp.push_back(closest[i].seq);  }
			temp.push_back(query);  
			
			createFilter(temp, 0.5);
		
			for (int i = 0; i < temp.size(); i++) { 
				if (m->control_pressed) {  return 0;  }
				runFilter(temp[i]);  
			}
			
			//update spotMap
			map<int, int> newMap;
			int spot = 0;
			
			for (int i = 0; i < filterString.length(); i++) {
				if (filterString[i] == '1') {
					//add to newMap
					newMap[spot] = spotMap[i];
					spot++;  
				}
			}
			spotMap = newMap;
		}

		//trim sequences - this follows ccodes remove_extra_gaps 
		trimSequences(query);
		if (m->control_pressed) {  return 0;  }
		
		//windows are equivalent to words - ccode paper recommends windows are between 5% and 20% on alignment length().  
		//Our default will be 10% and we will warn if user tries to use a window above or below these recommendations
		windows = findWindows();  
		if (m->control_pressed) {  return 0;  }

		//remove sequences that are more than 20% different and less than 0.5% different - may want to allow user to specify this later 
		removeBadReferenceSeqs(closest);
		if (m->control_pressed) {  return 0;  }
		
		//find the averages for each querys references
		getAverageRef(closest);  //fills sumRef, averageRef, sumSquaredRef and refCombo.
		getAverageQuery(closest, query);  //fills sumQuery, averageQuery, sumSquaredQuery.
		if (m->control_pressed) {  return 0;  }			
		
		//find the averages for each querys references 
		findVarianceRef();  //fills varRef and sdRef also sets minimum error rate to 0.001 to avoid divide by 0.
		if (m->control_pressed) {  return 0;  }	
			
		//find the averages for the query 
		findVarianceQuery();  //fills varQuery and sdQuery also sets minimum error rate to 0.001 to avoid divide by 0.
		if (m->control_pressed) {  return 0;  }
					
		determineChimeras();  //fills anova, isChimericConfidence, isChimericTStudent and isChimericANOVA. 
		if (m->control_pressed) {  return 0;  }
		
		return 0;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "getChimeras");
		exit(1);
	}
}
/***************************************************************************************************************/
//ccode algo says it does this to "Removes the initial and final gaps to avoid biases due to incomplete sequences."
void Ccode::trimSequences(Sequence* query) {
	try {
		
		int frontPos = 0;  //should contain first position in all seqs that is not a gap character
		int rearPos = query->getAligned().length();
		
		//********find first position in closest seqs that is a non gap character***********//
		//find first position all query seqs that is a non gap character
		for (int i = 0; i < closest.size(); i++) {
			
			string aligned = closest[i].seq->getAligned();
			int pos = 0;
			
			//find first spot in this seq
			for (int j = 0; j < aligned.length(); j++) {
				if (isalpha(aligned[j])) {
					pos = j;
					break;
				}
			}
			
			//save this spot if it is the farthest
			if (pos > frontPos) { frontPos = pos; }
		}
		
		//find first position all querySeq[query] that is a non gap character
		string aligned = query->getAligned();
		int pos = 0;
			
		//find first spot in this seq
		for (int j = 0; j < aligned.length(); j++) {
			if (isalpha(aligned[j])) {
				pos = j;
				break;
			}
		}
		
		//save this spot if it is the farthest
		if (pos > frontPos) { frontPos = pos; }
		
		
		//********find last position in closest seqs that is a non gap character***********//
		for (int i = 0; i < closest.size(); i++) {
			
			string aligned = closest[i].seq->getAligned();
			int pos = aligned.length();
			
			//find first spot in this seq
			for (int j = aligned.length()-1; j >= 0; j--) {
				if (isalpha(aligned[j])) {
					pos = j;
					break;
				}
			}
			
			//save this spot if it is the farthest
			if (pos < rearPos) { rearPos = pos; }
		}
		
		//find last position all querySeqs[query] that is a non gap character
		aligned = query->getAligned();
		pos = aligned.length();
		
		//find first spot in this seq
		for (int j = aligned.length()-1; j >= 0; j--) {
			if (isalpha(aligned[j])) {
				pos = j;
				break;
			}
		}
		
		//save this spot if it is the farthest
		if (pos < rearPos) { rearPos = pos; }

		
		//check to make sure that is not whole seq
		if ((rearPos - frontPos - 1) <= 0) {  m->mothurOut("Error, when I trim your sequences, the entire sequence is trimmed."); m->mothurOutEndLine(); exit(1);  }
		
		map<int, int> tempTrim;
		tempTrim[frontPos] = rearPos;
		
		//save trimmed locations
		trim = tempTrim;
						
		//update spotMask
		map<int, int> newMap;
		int spot = 0;
		
		for (int i = frontPos; i < rearPos; i++) {
			//add to newMap
			newMap[spot] = spotMap[i];
			spot++;  
		}
		spotMap = newMap;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "trimSequences");
		exit(1);
	}
}
/***************************************************************************************************************/
vector<int> Ccode::findWindows() {
	try {
		
		vector<int> win; 
		it = trim.begin();
		
		int length = it->second - it->first;
	
		//default is wanted = 10% of total length
		if (windowSizes > length) { 
			m->mothurOut("You have slected a window larger than your sequence length after all filters, masks and trims have been done. I will use the default 10% of sequence length.");
			windowSizes = length / 10;
		}else if (windowSizes == 0) { windowSizes = length / 10;  }
		else if (windowSizes > (length * 0.20)) {
			m->mothurOut("You have selected a window that is larger than 20% of your sequence length.  This is not recommended, but I will continue anyway."); m->mothurOutEndLine();
		}else if (windowSizes < (length * 0.05)) {
			m->mothurOut("You have selected a window that is smaller than 5% of your sequence length.  This is not recommended, but I will continue anyway."); m->mothurOutEndLine();
		}
		
		//save starting points of each window
		for (int m = it->first;  m < (it->second-windowSizes); m+=windowSizes) {  win.push_back(m);  }
		
		//save last window
		if (win[win.size()-1] < (it->first+length)) {
			win.push_back(win[win.size()-1]+windowSizes); // ex. string length is 115, window is 25, without this you would get 0, 25, 50, 75
		}																									//with this you would get 1,25,50,75,100
		
		return win;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "findWindows");
		exit(1);
	}
}
//***************************************************************************************************************
int Ccode::getDiff(string seqA, string seqB) {
	try {
		
		int numDiff = 0;
		
		for (int i = 0; i < seqA.length(); i++) {
			//if you are both not gaps
			//if (isalpha(seqA[i]) && isalpha(seqA[i])) {
				//are you different
				if (seqA[i] != seqB[i]) { 
					 int ok; /* ok=1 means equivalent base. Checks for degenerate bases */

					/* the char in base_a and base_b have been checked and they are different */
					if ((seqA[i] == 'N') && (seqB[i] != '-')) ok = 1;
					else if ((seqB[i] == 'N') && (seqA[i] != '-')) ok = 1;
					else if ((seqA[i] == 'Y') && ((seqB[i] == 'C') || (seqB[i] == 'T'))) ok = 1;
					else if ((seqB[i] == 'Y') && ((seqA[i] == 'C') || (seqA[i] == 'T'))) ok = 1;
					else if ((seqA[i] == 'R') && ((seqB[i] == 'G') || (seqB[i] == 'A'))) ok = 1;
					else if ((seqB[i] == 'R') && ((seqA[i] == 'G') || (seqA[i] == 'A'))) ok = 1;
					else if ((seqA[i] == 'S') && ((seqB[i] == 'C') || (seqB[i] == 'G'))) ok = 1;
					else if ((seqB[i] == 'S') && ((seqA[i] == 'C') || (seqA[i] == 'G'))) ok = 1;
					else if ((seqA[i] == 'W') && ((seqB[i] == 'T') || (seqB[i] == 'A'))) ok = 1;
					else if ((seqB[i] == 'W') && ((seqA[i] == 'T') || (seqA[i] == 'A'))) ok = 1;
					else if ((seqA[i] == 'M') && ((seqB[i] == 'A') || (seqB[i] == 'C'))) ok = 1;
					else if ((seqB[i] == 'M') && ((seqA[i] == 'A') || (seqA[i] == 'C'))) ok = 1;
					else if ((seqA[i] == 'K') && ((seqB[i] == 'T') || (seqB[i] == 'G'))) ok = 1;
					else if ((seqB[i] == 'K') && ((seqA[i] == 'T') || (seqA[i] == 'G'))) ok = 1;
					else if ((seqA[i] == 'V') && ((seqB[i] == 'C') || (seqB[i] == 'A') || (seqB[i] == 'G'))) ok = 1;
					else if ((seqB[i] == 'V') && ((seqA[i] == 'C') || (seqA[i] == 'A') || (seqA[i] == 'G'))) ok = 1;
					else if ((seqA[i] == 'H') && ((seqB[i] == 'T') || (seqB[i] == 'A') || (seqB[i] == 'C'))) ok = 1;
					else if ((seqB[i] == 'H') && ((seqA[i] == 'T') || (seqA[i] == 'A') || (seqA[i] == 'C'))) ok = 1;
					else if ((seqA[i] == 'D') && ((seqB[i] == 'T') || (seqB[i] == 'A') || (seqB[i] == 'G'))) ok = 1;
					else if ((seqB[i] == 'D') && ((seqA[i] == 'T') || (seqA[i] == 'A') || (seqA[i] == 'G'))) ok = 1;
					else if ((seqA[i] == 'B') && ((seqB[i] == 'C') || (seqB[i] == 'T') || (seqB[i] == 'G'))) ok = 1;
					else if ((seqB[i] == 'B') && ((seqA[i] == 'C') || (seqA[i] == 'T') || (seqA[i] == 'G'))) ok = 1;
					else ok = 0;  /* the bases are different and not equivalent */
					
					//check if they are both blanks
					if ((seqA[i] == '.') && (seqB[i] == '-')) ok = 1;
					else if ((seqB[i] == '.') && (seqA[i] == '-')) ok = 1;
					
					if (ok == 0) {  numDiff++;  }
				}
			//}
		}
		
		return numDiff;
	
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "getDiff");
		exit(1);
	}
}
//***************************************************************************************************************
//tried to make this look most like ccode original implementation
void Ccode::removeBadReferenceSeqs(vector<SeqDist>& seqs) {
	try {
		
		vector< vector<int> > numDiffBases;
		numDiffBases.resize(seqs.size());
		//initialize to 0
		for (int i = 0; i < numDiffBases.size(); i++) { numDiffBases[i].resize(seqs.size(),0); }
		
		it = trim.begin();
		int length = it->second - it->first;
		
		//calc differences from each sequence to everyother seq in the set
		for (int i = 0; i < seqs.size(); i++) {
			
			string seqA = seqs[i].seq->getAligned().substr(it->first, length);
			
			//so you don't calc i to j and j to i since they are the same
			for (int j = 0; j < i; j++) {
				
				string seqB = seqs[j].seq->getAligned().substr(it->first, length);
				
				//compare strings
				int numDiff = getDiff(seqA, seqB);
				
				numDiffBases[i][j] = numDiff;
				numDiffBases[j][i] = numDiff;
			}
		}
		
		//initailize remove to 0
		vector<int> remove;  remove.resize(seqs.size(), 0);
		float top = ((20*length) / (float) 100);
		float bottom = ((0.5*length) / (float) 100);
		
		//check each numDiffBases and if any are higher than threshold set remove to 1 so you can remove those seqs from the closest set
		for (int i = 0; i < numDiffBases.size(); i++) {
			for (int j = 0; j < i; j++) {	
				//are you more than 20% different
				if (numDiffBases[i][j] > top)		{  remove[j] = 1;  }
				//are you less than 0.5% different
				if (numDiffBases[i][j] < bottom)	{  remove[j] = 1;  }
			}
		}
		
		int numSeqsLeft = 0;
		
		//count seqs that are not going to be removed
		for (int i = 0; i < remove.size(); i++) {  
			if (remove[i] == 0)  { numSeqsLeft++;  }
		}
		
		//if you have enough then remove bad ones
		if (numSeqsLeft >= 3) {
			vector<SeqDist> goodSeqs;
			//remove bad seqs
			for (int i = 0; i < remove.size(); i++) {
				if (remove[i] == 0) { 
					goodSeqs.push_back(seqs[i]);
				}
			}
			
			seqs = goodSeqs;
			
		}else { //warn, but dont remove any
			m->mothurOut(querySeq->getName() + " does not have an adaquate number of reference sequences that are within 20% and 0.5% similarity.  I will continue, but please check."); m->mothurOutEndLine();  
		}

	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "removeBadReferenceSeqs");
		exit(1);
	}
}
//***************************************************************************************************************
//makes copy of templateseq for filter
vector<SeqDist>  Ccode::findClosest(Sequence* q, int numWanted) {
	try{
	
		vector<SeqDist>  topMatches;  
		
		Sequence query = *(q);
			
		//calc distance to each sequence in template seqs
		for (int i = 0; i < templateSeqs.size(); i++) {
			
			Sequence ref = *(templateSeqs[i]); 
					
			//find overall dist
			distCalc->calcDist(query, ref);
			float dist = distCalc->getDist();	
				
			//save distance
			SeqDist temp;
			temp.seq = new Sequence(templateSeqs[i]->getName(), templateSeqs[i]->getAligned());
			temp.dist = dist;

			topMatches.push_back(temp);
		}
			
		sort(topMatches.begin(), topMatches.end(), compareSeqDist);
		
		for (int i = numWanted; i < topMatches.size(); i++) {  delete topMatches[i].seq;  }
		
		topMatches.resize(numWanted);
			
		return topMatches;

	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "findClosestSides");
		exit(1);
	}
}
/**************************************************************************************************/
//find the distances from each reference sequence to every other reference sequence for each window for this query
void Ccode::getAverageRef(vector<SeqDist> ref) {
	try {
		
		vector< vector< vector<int> > >  diffs;  //diffs[0][1][2] is the number of differences between ref seq 0 and ref seq 1 at window 2.
		
		//initialize diffs vector
		diffs.resize(ref.size());
		for (int i = 0; i < diffs.size(); i++) {  
			diffs[i].resize(ref.size());
			for (int j = 0; j < diffs[i].size(); j++) {
				diffs[i][j].resize(windows.size(), 0);
			}
		}
		
		it = trim.begin();
				
		//find the distances from each reference sequence to every other reference sequence for each window for this query		
		for (int i = 0; i < ref.size(); i++) {
			
			string refI = ref[i].seq->getAligned();
			
			//j<i, so you don't find distances from i to j and then j to i.
			for (int j = 0; j < i; j++) {
			
				string refJ = ref[j].seq->getAligned();
			
				for (int k = 0; k < windows.size(); k++) {
					
					string refIWindowk, refJWindowk;
					
					if (k < windows.size()-1) {
						//get window strings
						refIWindowk = refI.substr(windows[k], windowSizes);
						refJWindowk = refJ.substr(windows[k], windowSizes);
					}else { //last window may be smaller than rest - see findwindows
						//get window strings
						refIWindowk = refI.substr(windows[k], (it->second-windows[k]));
						refJWindowk = refJ.substr(windows[k], (it->second-windows[k]));
					}
					
					//find differences
					int diff = getDiff(refIWindowk, refJWindowk);

					//save differences in [i][j][k] and [j][i][k] since they are the same
					diffs[i][j][k] = diff;
					diffs[j][i][k] = diff;

				}//k
				
			}//j
		
		}//i
		
		//initialize sumRef for this query 
		sumRef.resize(windows.size(), 0);
		sumSquaredRef.resize(windows.size(), 0);
		averageRef.resize(windows.size(), 0);
		
		//find the sum of the differences for hte reference sequences
		for (int i = 0; i < diffs.size(); i++) {
			for (int j = 0; j < i; j++) {
			
				//increment this querys reference sequences combos
				refCombo++;
				
				for (int k = 0; k < diffs[i][j].size(); k++) {
					sumRef[k] += diffs[i][j][k];
					sumSquaredRef[k] += (diffs[i][j][k]*diffs[i][j][k]);
				}//k
				
			}//j
		}//i

		
		//find the average of the differences for the references for each window
		for (int i = 0; i < windows.size(); i++) {
			averageRef[i] = sumRef[i] / (float) refCombo;
		}
		
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "getAverageRef");
		exit(1);
	}
}
/**************************************************************************************************/
void Ccode::getAverageQuery (vector<SeqDist> ref, Sequence* query) {
	try {
	
		vector< vector<int> >  diffs;  //diffs[1][2] is the number of differences between querySeqs[query] and ref seq 1 at window 2.
		
		//initialize diffs vector
		diffs.resize(ref.size());
		for (int j = 0; j < diffs.size(); j++) {
			diffs[j].resize(windows.size(), 0);
		}
		
		it = trim.begin();
							
		string refQuery = query->getAligned();
			 
		//j<i, so you don't find distances from i to j and then j to i.
		for (int j = 0; j < ref.size(); j++) {
			 
			 string refJ = ref[j].seq->getAligned();
			 
			 for (int k = 0; k < windows.size(); k++) {
					
					string QueryWindowk, refJWindowk;
					
					if (k < windows.size()-1) {
						//get window strings
						QueryWindowk = refQuery.substr(windows[k], windowSizes);
						refJWindowk = refJ.substr(windows[k], windowSizes);					
					}else { //last window may be smaller than rest - see findwindows
						//get window strings
						QueryWindowk = refQuery.substr(windows[k], (it->second-windows[k]));
						refJWindowk = refJ.substr(windows[k], (it->second-windows[k]));
					}
					
					//find differences
					int diff = getDiff(QueryWindowk, refJWindowk);
					
					//save differences 
					diffs[j][k] = diff;
			 
			 }//k
		}//j
			 
		
		//initialize sumRef for this query 
		sumQuery.resize(windows.size(), 0);
		sumSquaredQuery.resize(windows.size(), 0);
		averageQuery.resize(windows.size(), 0);
		
		//find the sum of the differences 
		for (int j = 0; j < diffs.size(); j++) {
			for (int k = 0; k < diffs[j].size(); k++) {
				sumQuery[k] += diffs[j][k];
				sumSquaredQuery[k] += (diffs[j][k]*diffs[j][k]);
			}//k
		}//j
		
		
		//find the average of the differences for the references for each window
		for (int i = 0; i < windows.size(); i++) {
			averageQuery[i] = sumQuery[i] / (float) ref.size();
		}
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "getAverageQuery");
		exit(1);
	}
}
/**************************************************************************************************/
void Ccode::findVarianceRef() {
	try {
		
		varRef.resize(windows.size(), 0);
		sdRef.resize(windows.size(), 0);
		
		//for each window
		for (int i = 0; i < windows.size(); i++) {
			varRef[i] = (sumSquaredRef[i] - ((sumRef[i]*sumRef[i])/(float)refCombo)) / (float)(refCombo-1);
			sdRef[i] = sqrt(varRef[i]);
			
			//set minimum error rate to 0.001 - to avoid potential divide by zero - not sure if this is necessary but it follows ccode implementation
			if (averageRef[i] < 0.001)			{	averageRef[i] = 0.001;		}
			if (sumRef[i] < 0.001)				{	sumRef[i] = 0.001;			}
			if (varRef[i] < 0.001)				{	varRef[i] = 0.001;			}
			if (sumSquaredRef[i] < 0.001)		{	sumSquaredRef[i] = 0.001;	}
			if (sdRef[i] < 0.001)				{	sdRef[i] = 0.001;			}
				
		}
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "findVarianceRef");
		exit(1);
	}
}
/**************************************************************************************************/
void Ccode::findVarianceQuery() {
	try {
		varQuery.resize(windows.size(), 0);
		sdQuery.resize(windows.size(), 0);
		
		//for each window
		for (int i = 0; i < windows.size(); i++) {
			varQuery[i] = (sumSquaredQuery[i] - ((sumQuery[i]*sumQuery[i])/(float) closest.size())) / (float) (closest.size()-1);
			sdQuery[i] = sqrt(varQuery[i]);
			
			//set minimum error rate to 0.001 - to avoid potential divide by zero - not sure if this is necessary but it follows ccode implementation
			if (averageQuery[i] < 0.001)			{	averageQuery[i] = 0.001;		}
			if (sumQuery[i] < 0.001)				{	sumQuery[i] = 0.001;			}
			if (varQuery[i] < 0.001)				{	varQuery[i] = 0.001;			}
			if (sumSquaredQuery[i] < 0.001)		{	sumSquaredQuery[i] = 0.001;	}
			if (sdQuery[i] < 0.001)				{	sdQuery[i] = 0.001;			}
		}

	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "findVarianceQuery");
		exit(1);
	}
}
/**************************************************************************************************/
void Ccode::determineChimeras() {
	try {
		
		isChimericConfidence.resize(windows.size(), false);
		isChimericTStudent.resize(windows.size(), false);
		isChimericANOVA.resize(windows.size(), false);
		anova.resize(windows.size());

		
		//for each window
		for (int i = 0; i < windows.size(); i++) {
		
			//get confidence limits
			float t = getT(closest.size()-1);  //how many seqs you are comparing to this querySeq
			float dsUpper = (averageQuery[i] + (t * sdQuery[i])) / averageRef[i]; 
			float dsLower = (averageQuery[i] - (t * sdQuery[i])) / averageRef[i]; 
		
			if ((dsUpper > 1.0) && (dsLower > 1.0) && (averageQuery[i] > averageRef[i])) {  /* range does not include 1 */
					isChimericConfidence[i] = true;   /* significantly higher at P<0.05 */ 

			}
			
			//student t test
			int degreeOfFreedom = refCombo + closest.size() - 2;
			float denomForT = (((refCombo-1) * varQuery[i] + (closest.size() - 1) * varRef[i]) / (float) degreeOfFreedom) * ((refCombo + closest.size()) / (float) (refCombo * closest.size()));	/* denominator, without sqrt(), for ts calculations */
				
			float ts = fabs((averageQuery[i] - averageRef[i]) / (sqrt(denomForT)));  /* value of ts for t-student test */	
			t = getT(degreeOfFreedom);
	
			if ((ts >= t) && (averageQuery[i] > averageRef[i])) {   
				isChimericTStudent[i] = true;   /* significantly higher at P<0.05 */ 
			}
		
			//anova test
			float value1 = sumQuery[i] + sumRef[i];
			float value2 = sumSquaredQuery[i] + sumSquaredRef[i];
			float value3 = ((sumQuery[i]*sumQuery[i]) / (float) (closest.size())) + ((sumRef[i] * sumRef[i]) / (float) refCombo);
			float value4 = (value1 * value1) / ( (float) (closest.size() + refCombo) );
			float value5 = value2 - value4;
			float value6 = value3 - value4;
			float value7 = value5 - value6;
			float value8 = value7 / ((float) degreeOfFreedom);
			float anovaValue = value6 / value8;
			
			float f = getF(degreeOfFreedom);
			
			 if ((anovaValue >= f) && (averageQuery[i] > averageRef[i]))  {
	               isChimericANOVA[i] = true;   /* significant P<0.05 */
	        }
			
			if (isnan(anovaValue) || isinf(anovaValue)) { anovaValue = 0.0; }
			
			anova[i] = anovaValue;
		}
		
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "determineChimeras");
		exit(1);
	}
}
/**************************************************************************************************/	
float Ccode::getT(int numseq) {
	try {
	
		float tvalue = 0;
		
		/* t-student critical values for different degrees of freedom and alpha 0.1 in one-tail tests (equivalent to 0.05) */
		if (numseq > 120) tvalue = 1.645;
		else if (numseq > 60) tvalue = 1.658;
        else if (numseq > 40) tvalue = 1.671;
        else if (numseq > 30) tvalue = 1.684;
        else if (numseq > 29) tvalue = 1.697;
        else if (numseq > 28) tvalue = 1.699;
        else if (numseq > 27) tvalue = 1.701;
        else if (numseq > 26) tvalue = 1.703;
        else if (numseq > 25) tvalue = 1.706;
        else if (numseq > 24) tvalue = 1.708;
        else if (numseq > 23) tvalue = 1.711;
        else if (numseq > 22) tvalue = 1.714;
        else if (numseq > 21) tvalue = 1.717;
        else if (numseq > 20) tvalue = 1.721;
        else if (numseq > 19) tvalue = 1.725;
        else if (numseq > 18) tvalue = 1.729;
        else if (numseq > 17) tvalue = 1.734;
        else if (numseq > 16) tvalue = 1.740;
        else if (numseq > 15) tvalue = 1.746;
        else if (numseq > 14) tvalue = 1.753;
        else if (numseq > 13) tvalue = 1.761;
        else if (numseq > 12) tvalue = 1.771;
        else if (numseq > 11) tvalue = 1.782;
        else if (numseq > 10) tvalue = 1.796;
        else if (numseq > 9) tvalue = 1.812;
        else if (numseq > 8) tvalue = 1.833;
        else if (numseq > 7) tvalue = 1.860;
        else if (numseq > 6) tvalue = 1.895;
        else if (numseq > 5) tvalue = 1.943;
        else if (numseq > 4) tvalue = 2.015;
        else if (numseq > 3) tvalue = 2.132;
        else if (numseq > 2) tvalue = 2.353;
        else if (numseq > 1) tvalue = 2.920;
		else if (numseq <= 1) {
			m->mothurOut("Two or more reference sequences are required, your data will be flawed.\n"); m->mothurOutEndLine();
		}
		
		return tvalue;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "getT");
		exit(1);
	}
}
/**************************************************************************************************/	
float Ccode::getF(int numseq) {
	try {
	
		float fvalue = 0;
		
		 /* F-Snedecor critical values for v1=1 and different degrees of freedom v2 and alpha 0.05 */
        if (numseq > 120) fvalue = 3.84;
        else if (numseq > 60) fvalue = 3.92;
        else if (numseq > 40) fvalue = 4.00;
        else if (numseq > 30) fvalue = 4.08;
        else if (numseq > 29) fvalue = 4.17;
        else if (numseq > 28) fvalue = 4.18;
        else if (numseq > 27) fvalue = 4.20;
        else if (numseq > 26) fvalue = 4.21;
        else if (numseq > 25) fvalue = 4.23;
        else if (numseq > 24) fvalue = 4.24;
        else if (numseq > 23) fvalue = 4.26;
        else if (numseq > 22) fvalue = 4.28;
        else if (numseq > 21) fvalue = 4.30;
        else if (numseq > 20) fvalue = 4.32;
        else if (numseq > 19) fvalue = 4.35;
        else if (numseq > 18) fvalue = 4.38;
        else if (numseq > 17) fvalue = 4.41;
        else if (numseq > 16) fvalue = 4.45;
        else if (numseq > 15) fvalue = 4.49;
        else if (numseq > 14) fvalue = 4.54;
        else if (numseq > 13) fvalue = 4.60;
        else if (numseq > 12) fvalue = 4.67;
        else if (numseq > 11) fvalue = 4.75;
        else if (numseq > 10) fvalue = 4.84;
        else if (numseq > 9) fvalue = 4.96;
        else if (numseq > 8) fvalue = 5.12;
        else if (numseq > 7) fvalue = 5.32;
        else if (numseq > 6) fvalue = 5.59;
        else if (numseq > 5) fvalue = 5.99;
        else if (numseq > 4) fvalue = 6.61;
        else if (numseq > 3) fvalue = 7.71;
        else if (numseq > 2) fvalue = 10.1;
        else if (numseq > 1) fvalue = 18.5;
        else if (numseq > 0) fvalue = 161;
		else if (numseq <= 0) {
			m->mothurOut("Two or more reference sequences are required, your data will be flawed.\n"); m->mothurOutEndLine();
        }
		
		return fvalue;
	}
	catch(exception& e) {
		m->errorOut(e, "Ccode", "getF");
		exit(1);
	}
}
//***************************************************************************************************************