/***********************************************************************
* Filename: LOX.cpp
* Abstract: Level Of eXpression

* Version: 1.8
* Author: Zhang Zhang  (zhangzhang@big.ac.cn)
* Date: February 4, 2012
* Note: (1) Remove the bug for calculating p-values.
        (2) Change the first column name to be outputted, "ID" into "Gene".
        (3) Output p-values files for methodologies.  

* Version: 1.7
* Author: Zhang Zhang  (zhangzhang.cn@gmail.com)
* Date: July 26, 2011
* Note: Add an option to output posterior distribution for each gene.

* Version: 1.6
* Author: Zhang Zhang  (zhang.zhang@kaust.edu.sa)
* Date: Jun. 3, 2011
* Note: Remove the bug for calculating p-values.

* Version: 1.5
* Author: Zhang Zhang  (zhang.zhang@kaust.edu.sa)
* Date: Mar. 8, 2011
* Note: (1) Remove the bug when first column count=0, resulting from locating
        the lowest value and normarlization (getPQResults).
		(2) Provide more specifica error message for the input file (setHeader
		and Run). 

* Version: 1.4
* Author: Zhang Zhang  (zhang.zhang@kaust.edu.sa)
* Date: Nov. 9, 2010
* Note: Remove the bug for normalizing of p.

* Version: 1.3
* Author: Zhang Zhang  (zhang.zhang@kaust.edu.sa)
* Date: Nov. 1, 2010
* Note: Set q=1 for any input file with only one methodology.

* Version: 1.2
* Author: Zhang Zhang  (zhang.zhang@kaust.edu.sa)
* Date: July, 2010
* Note: (1) Accommodate input file with only one methodology;
        (2) Update the LOX reference published in Bioinformatics.

* Version: 1.1
* Author: Zhang Zhang  (zhang.zhang@kaust.edu.sa)
* Date: April, 2010
* Note: Output pvalue files for treatments

* Version: 1.0
* Author: Zhang Zhang  (zhang.zhang@yale.edu)
* Date: August 27, 2009

* Note: 
LOX is distributed as open-source software and licensed under the GNU 
General Public License (Version 3; http://www.gnu.org/licenses/gpl.txt), 
in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even 
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
See the GNU General Public License for more details. Commercial use of LOX 
requires a special contract.
************************************************************************/

#include "LOX.h"

LOX::LOX() {
	isOutputMethodsResults = false;
	isShowMLERatios = false;
	isOutputPosterior = false;
	InputFileName = "";
	LOXOutputFileName = "";
	MetOutputFileName = "";
	NumOfColumn = 0;

	ChainSampleSize = 1000;
	BurnIn = 2000;
	SamplingPeriod = 20;
	MaxTuningDepth = 8;

	LowerAcceptance = 0.15;
	UpperAcceptance = 0.50;
}

LOX::~LOX() {
}


/*****************************************************
* Name: setExpressionNodes
* Function: Parse the expression counts for each gene
* Input Parameter: 
* Output: 
* Return Value: int, 1=Success, 0=Fail
*****************************************************/
int LOX::setExpressionNodes() {

	int OK = 1;

	try {
		cout<<"Assigning expression counts...";

		unsigned long i, j, line_number = 3;	

		for (i=line_number; i<DatLines.size(); i++) {
			
			vector<string> cells = splitString(DatLines[i]);

			string id = "";
			string common_name = "";		
			vector<double> nodes;
			string msg = "";

			j = 0;
			if (cells.size()>2) {			
				id = cells[j++];
				common_name = cells[j++];
			}		
			for (; j<cells.size(); j++) {
				string t = cells[j];
				if (cells[j]=="") {
					cells[j] = "0";
					msg = "Count can not be empty.";
				}
				nodes.push_back(CONVERT<double>(cells[j]));
			}

			if(cells.size()!=NumOfColumn) msg = "Error input data";

			LOXNode lxn(id, common_name, nodes, msg);
			NodesVector.push_back(lxn);
		}

		//Free memory
		DatLines.clear();

		//Calculate total counts for each column
		for (i=0; i<HeadersVector.size(); i++) {
			double sum = 0.0;
			for (j=0; j<NodesVector.size(); j++) sum += NodesVector[j].Nodes[i];
			ColumnSumsVector.push_back(sum);		
		}
		
		//Calculate total counts for each treatment across all genes
		for (i=0; i<TreatmentsVector.size(); i++) TreatmentSumsVector.push_back(0);
		for (i=0; i<HeadersVector.size(); i++) TreatmentSumsVector[HeadersVector[i].treatmentNum] += ColumnSumsVector[i];	

		//Calculate total counts for each method across all genes
		for (i=0; i<MethodsVector.size(); i++) MethodSumsVector.push_back(0);
		for (i=0; i<HeadersVector.size(); i++) MethodSumsVector[HeadersVector[i].methodNum] += ColumnSumsVector[i];

		//Calculate total counts for each treatment/method under each gene; version 1.4
		for (i=0; i<NodesVector.size(); i++) {
			//initialize for treatments
			for (j=0; j<TreatmentsVector.size(); j++) NodesVector[i].NodesOrderByTreatment.push_back(0);
			//initialize for methods
			for (j=0; j<MethodsVector.size(); j++) NodesVector[i].NodesOrderByMethod.push_back(0);
			//calculate the total count
			for (j=0; j<HeadersVector.size(); j++) {
				NodesVector[i].NodesOrderByTreatment[HeadersVector[j].treatmentNum] += NodesVector[i].Nodes[j];
				NodesVector[i].NodesOrderByMethod[HeadersVector[j].methodNum] += NodesVector[i].Nodes[j];
			}
		}
	}
	catch (...) {
		OK = 0;
	}

	return OK;
}

int LOX::showParameterSettings() {
	cout<<ABBRNAME<<" ("<<FULLNAME<<"), Version: "<<VERSION<<" ["<<LASTUPDATE<<"]"<<endl;	
	cout<<"REFERENCE: "<<REFERENCE<<endl;
	cout<<endl;

	cout<<"Parameter settings:"<<endl;
	cout<<"Output methodology results = "<<toString(isOutputMethodsResults)<<endl;
	cout<<"Show results by maximum likelihood estimation = "<<toString(isShowMLERatios)<<endl;
	cout<<"Output posterior distribution for each gene = "<<toString(isOutputPosterior)<<endl;
	cout<<"Chain sample size = "<<ChainSampleSize<<endl;
	cout<<"Burnin size = "<<BurnIn<<endl;
	cout<<"Sampling period = "<<SamplingPeriod<<endl;
	cout<<"Tuning depth = "<<MaxTuningDepth<<endl;
	cout<<"Lower acceptance rate = "<<LowerAcceptance<<endl;
	cout<<"Upper acceptance rate = "<<UpperAcceptance<<endl;
	
	cout<<endl;

	return 1;
}

/*****************************************************
* Name: setHeader
* Function: Parse the header (the first 3 lines) and 
				set treatments and methods names
* Input Parameter: 
* Output: treatments and methods names
* Return Value: int, 1=Success, 0=Fail
*****************************************************/
int LOX::setHeader() {
	
	int OK = 1;	

	try {
		string firstline =  DatLines[0];
		string secondline = DatLines[1];
		string thirdline =  DatLines[2];
		string fourthline =  DatLines[3];

		//cout<<firstline.c_str()<<endl;
		
		//Read line NO.2 and NO.3
		//cout<<"Initializing ExpressionNodeNameList..."<<endl;	
		vector<string> cells_No2Line = splitString(secondline);
		vector<string> cells_No3Line = splitString(thirdline);
		vector<string> cells_No4Line = splitString(fourthline);

		if(cells_No2Line.size()!=cells_No3Line.size()) {
			cout<<"2nd line has "<<cells_No2Line.size()<<" columns, whereas 3rd line has "<<cells_No3Line.size()<<" columns. They should be equal in column count."<<endl;
			throw 1;
		}

		if(cells_No3Line.size()!=cells_No4Line.size()) {
			cout<<"3rd line has "<<cells_No3Line.size()<<" columns, whereas 4th line has "<<cells_No4Line.size()<<" columns. They should be equal in column count."<<endl;
			throw 1;
		}

		NumOfColumn = cells_No2Line.size();

		int count1 = 0, count2 = 0;
		string Notes1 = cells_No2Line[count1++];
		string Channel1 = cells_No2Line[count1++];

		string Notes2 = cells_No3Line[count2++];
		string Channel2 = cells_No3Line[count2++];

		cout<<"Assigning names to expression methodologies and treatments...";

		while(count1 < cells_No2Line.size()) {		
			
			//Methods
			string name_method = trim(cells_No2Line[count1++]);		
			cout<<".";		
			
			int num_method = 0;			
			bool MemorizedYet = false;
			while (num_method < MethodsVector.size() && MemorizedYet==false) {			
				if (MethodsVector[num_method]==name_method) MemorizedYet = true;
				else num_method++;
			}
					
			if (MemorizedYet==false) MethodsVector.push_back(name_method);

			//Treatments
			string name_treatment = trim(cells_No3Line[count2++]);
			cout<<".";
			
			int num_treatment = 0;
			MemorizedYet = false;

			while (num_treatment < TreatmentsVector.size() && MemorizedYet==false) {			
				if (TreatmentsVector[num_treatment]==name_treatment) MemorizedYet = true;
				else num_treatment++;
			}
			
			if (MemorizedYet==false) TreatmentsVector.push_back(name_treatment);

			//Put values to HeadersVector
			HeadersVector.push_back(LOXHeader(name_method, num_method, name_treatment, num_treatment));
		}

		int Node;
		cout<<endl<<"Expression methodologies:";
		for (Node = 0; Node < MethodsVector.size(); Node++) cout<<"\t"<<MethodsVector[Node].c_str();

		cout<<endl<<"Expression treatments:";
		for (Node = 0; Node < TreatmentsVector.size(); Node++) cout<<"\t"<<TreatmentsVector[Node].c_str();
		
		cout<<endl;
	}
	catch (...) {
		OK = 0;
	}

	return OK;	
}

/*****************************************************
* Name: calculateLikelihood
* Function: Calculate the log-likelihood for each gene
* Input Parameter: lox, p, q
* Output: 
* Return Value: double, log-likelihood value
* Note: input parameter vector p, q are changeable
*****************************************************/
double LOX::calculateLikelihood(LOXNode lox, vector<double>& p, vector<double>& q) {
	double L = 0.0;

	try {
		int i;
		for (i = 0; i < HeadersVector.size(); i++) {
			//set p=0 if the count is zero; version 1.3
			if (CONVERT<int>(lox.NodesOrderByTreatment[HeadersVector[i].treatmentNum])==0) {
				p[HeadersVector[i].treatmentNum] = 0.0;
			}
			double zz = p[HeadersVector[i].treatmentNum];
			
			//set q=1 in only one method available; version 1.3
			if (MethodSumsVector.size()==1) {
				q[HeadersVector[i].methodNum] = 1.0;
			}

			//zz = p*q
			zz *=  q[HeadersVector[i].methodNum];
			if(zz > SMALLVALUE && zz < 1.0) {
				L += lox.Nodes[i]*log(zz) + (ColumnSumsVector[i] - lox.Nodes[i])*log(1.0 - zz);
			}			
		}
	}
	catch(...) {
	}

	return L;
}

/*****************************************************
* Name: randomPQ
* Function: Randomly change either P(treatments) or Q(methods)
* Input Parameter: pORq, step
* Output: 
* Return Value: int, Success if Prior==1
* Note: pORq will be changed
*****************************************************/
int LOX::randomPQ(vector<double>& pORq, vector<double> step) {

	int Prior = 1;
	//Accept the input file with only one method "num=1"; version 1.2
	int num = pORq.size();

	int FirstNodeChanged = getRandLow2High(0, num);
	int SecondNodeChanged = getRandLow2High(0, num);
	while (SecondNodeChanged==FirstNodeChanged && num>1) {
		SecondNodeChanged = getRandLow2High(0, num);
	}

	double RandomWalk = getRandNeg1toPos1() * (step[FirstNodeChanged] + step[SecondNodeChanged]) / 2.0;
	pORq[FirstNodeChanged]  += RandomWalk;
	pORq[SecondNodeChanged] -= RandomWalk;
	if (pORq[FirstNodeChanged] < SMALLVALUE || pORq[SecondNodeChanged] < SMALLVALUE 
		|| pORq[FirstNodeChanged] > 1.0 || pORq[SecondNodeChanged] > 1.0) {
		Prior = 0;
	}
	//version 1.4: Remove the bug for normalization
	
	return Prior;
}

/*****************************************************
* Name: randomStepPQ
* Function: Randomly change steps for either P(treatments) or Q(methods)
* Input Parameter: MaxTuningDepth, acceptance_ratio, step, tuned
* Output: 
* Return Value: int
* Note: step and tuned will be changed
*****************************************************/
int LOX::randomStepPQ(int MaxTuningDepth, double acceptance_ratio, vector<double>& step, bool& tuned) {

	if (acceptance_ratio > LowerAcceptance && acceptance_ratio < UpperAcceptance)  {
		tuned = true;
	}		
	else {
		if (MaxTuningDepth > 1) {
			int i;
			if (acceptance_ratio < LowerAcceptance) {//make step smaller
				for (i=0; i<step.size(); i++) {
					step[i] /= (1 + (getRand0to1() + getRand0to1()) / 2 * ((1 - acceptance_ratio) / 0.5) * (1.5 - acceptance_ratio) * (1.5 - acceptance_ratio) * (1.5 - acceptance_ratio));
				}
			}
			else if (acceptance_ratio > UpperAcceptance) {//make step larger
				for (i=0; i<step.size(); i++) {
					step[i] *= (1 + (getRand0to1() + getRand0to1()) / 2 * (acceptance_ratio/0.15 - 1) * (acceptance_ratio/0.15 - 1));
				}
			}
		}
	}

	return 1;
}

/*****************************************************
* Name: MCMC
* Function: Calculate expression levels by MCMC
* Input Parameter: tally
* Output: Success=[OK], Fail=[Error]
* Return Value: int
*****************************************************/
int LOX::MCMC(unsigned long tally) {

	double DensityOfLikelihood, NewDensityOfLikelihood, BestDensityOfLikelihood;
	double PAcceptanceRatio, QAcceptanceRatio;
	long PStepsAccepted, QStepsAccepted;
	bool PTuned = false, QTuned = false;
	int TuningDepth = 0;
	int Prior = 1;
	int Generation, TotalGenerationsSampled, NumSampled;
	bool Burnt;
	unsigned long Node, Node1, Node2;

	int NumOfTreatments = TreatmentsVector.size();
	int NumOfMethods = MethodsVector.size();

	vector<DoubleVector> PTally(NumOfTreatments), QTally(NumOfMethods);

	//P-values for treatments, version 1.8
	vector<double> PMatrix(NumOfTreatments*NumOfTreatments, 0.0);
	//P-values for methodologies, version 1.8
	vector<double> PMethod(NumOfMethods*NumOfMethods, 0.0);

	DensityOfLikelihood = 0.0;
	
	try {
		LOXNode lox = NodesVector[tally];
		cout<<"["<<tally+1<<"] ";
		cout<<lox.ID.c_str()<<"  "<<lox.CommonName.c_str();

		if (lox.ErrorMsg!="") {
			cout<<"\t"<<lox.ErrorMsg.c_str();
			throw 1;
		}
		
		vector<double> best_p(NumOfTreatments, 0.0), best_q(NumOfMethods, 0.0);
		vector<double> step_p(NumOfTreatments, 0.0), step_q(NumOfMethods, 0.0);

		do {
			vector<double> p(NumOfTreatments, 0.0), q(NumOfMethods, 0.0);
			
			TuningDepth++;
			
			for (Node = 0; Node < HeadersVector.size(); Node++) {
				int treatmentNum = HeadersVector[Node].treatmentNum;
				double tmp_p = (double)lox.Nodes[Node] / TreatmentSumsVector[treatmentNum];
				p[treatmentNum] += tmp_p;

				int methodNum = HeadersVector[Node].methodNum;
				double tmp_q = (double)lox.Nodes[Node] / MethodSumsVector[methodNum];
				q[methodNum] += tmp_q;
			}

			if (TuningDepth==1) {
				for (Node=0; Node<NumOfTreatments; Node++) step_p[Node] = 0.5 * p[Node];
				for (Node=0; Node<NumOfMethods; Node++) step_q[Node] = 0.5 * q[Node];
			}

			DensityOfLikelihood = calculateLikelihood(lox, p, q);
			
			//Initial setting for the density of likelihood		
			NewDensityOfLikelihood = DensityOfLikelihood;
			BestDensityOfLikelihood = DensityOfLikelihood;
			
			TotalGenerationsSampled = 0;
			NumSampled = 0;
			Burnt = false;
			PStepsAccepted = 0;
			QStepsAccepted = 0;
			
			string posterior = "Sample #\tLog-Density of Likelihood";
			for(Node=0; Node<NumOfTreatments; Node++) posterior += "\tTreatment("+TreatmentsVector[Node]+")";
			for(Node=0; Node<NumOfMethods; Node++) posterior += "\tMethod("+MethodsVector[Node]+")";
			posterior += "\n";

			vector<double> new_p(NumOfTreatments, 0.0),  new_q(NumOfMethods, 0.0);
			
			//Init setting for p-values of treatments and methodologies
			for (Node=0; Node<NumOfTreatments*NumOfTreatments; Node++) PMatrix[Node] = 0.0;
			for (Node=0; Node<NumOfMethods*NumOfMethods;       Node++) PMethod[Node] = 0.0;

			do {
				Generation = 0;
				
				do {
					Generation++;
					TotalGenerationsSampled++;
					
					Prior = 1;
					copyVector(p, best_p);
					copyVector(q, best_q);
					
					copyVector(p, new_p);
					copyVector(q, new_q);
					
					//Change P, if odd values or input data has only one method
					if (Generation % 2 == 1 || MethodSumsVector.size()==1) {
						Prior = randomPQ(new_p, step_p);
					}
					//Change Q
					else {
						Prior = randomPQ(new_q, step_q);
					}

					//new_q[0] = 1;
					if (Prior==1) {
						NewDensityOfLikelihood = calculateLikelihood(lox, new_p, new_q);				
						if (log(getRand0to1()) < NewDensityOfLikelihood - DensityOfLikelihood ) {
							if (Generation % 2 == 1) PStepsAccepted++;
							else QStepsAccepted++;

							DensityOfLikelihood = NewDensityOfLikelihood;
							
							copyVector(new_p, p);
							copyVector(new_q, q);					

							if (DensityOfLikelihood >= BestDensityOfLikelihood) {
								BestDensityOfLikelihood = DensityOfLikelihood;
								copyVector(new_p, best_p);
								copyVector(new_q, best_q);
							}
						}
					}

				} while ((Generation!=BurnIn) && (!Burnt || Generation!=SamplingPeriod));

				if (!Burnt) {
					PStepsAccepted = 0;
					QStepsAccepted = 0;
					TotalGenerationsSampled = 0;
				}

				Burnt = true;								
				
				//Record posterior distribution
				posterior += toString(NumSampled) + "\t" + toString(DensityOfLikelihood); 
				for (Node=0; Node<NumOfTreatments; Node++) {
					posterior += "\t" + toString(p[Node]);
				}
				for (Node=0; Node<NumOfMethods; Node++) {
					posterior += "\t" + toString(q[Node]);
				}
				posterior += "\n";

				//Record p and q into PTally and QTally, respectively
				for (Node = 0; Node < NumOfTreatments; Node++) PTally[Node].push_back(p[Node]);			
				for (Node = 0; Node < NumOfMethods; Node++) QTally[Node].push_back(q[Node]);

				//Record P-values for treatments
				for (Node1 = 0; Node1 < NumOfTreatments; Node1++) {
					for (Node2 = 0; Node2 < NumOfTreatments; Node2++) {
						//Version 1.6: ">" is changed as ">=", for the case with p[node1]=p[node2].
						if (p[Node1] >= p[Node2]) PMatrix[Node1 * NumOfTreatments + Node2]++;
					}
				}				
				
				//Record P-values for methodologies, version 1.8
				for (Node1 = 0; Node1 < NumOfMethods; Node1++) {
					for (Node2 = 0; Node2 < NumOfMethods; Node2++) {
						if (q[Node1] >= q[Node2]) PMethod[Node1 * NumOfMethods + Node2]++;
					}
				}
				
				NumSampled++;

			} while (NumSampled!=ChainSampleSize);

			//Change step_p according to PAcceptanceRatio
			PAcceptanceRatio = ((double)PStepsAccepted * 2.0) / TotalGenerationsSampled;
			randomStepPQ(MaxTuningDepth, PAcceptanceRatio, step_p, PTuned);
			
			//Change step_q according to QAcceptanceRatio
			QAcceptanceRatio = ((double)QStepsAccepted * 2.0) / TotalGenerationsSampled;
			randomStepPQ(MaxTuningDepth, QAcceptanceRatio, step_q, QTuned);

			//Output posterior distribution
			if (isOutputPosterior) {
				string folder = "Posterior_" + InputFileName;
				#ifdef WIN32
					mkdir(folder.c_str());
				#else
					mkdir(folder.c_str(), 0755);
				#endif
				chdir(folder.c_str());
				writeFile(lox.ID, posterior.c_str(), 0);
				chdir("..");
			}

		} while ((TuningDepth!=MaxTuningDepth) && (!PTuned || !QTuned));

		//Put results into NodesVector
		getPQResults(NodesVector[tally], NodesVector[tally].NodesOrderByTreatment, PTally, best_p, PAcceptanceRatio, PTuned);
		getPQResults(NodesVector[tally], NodesVector[tally].NodesOrderByMethod,    QTally, best_q, QAcceptanceRatio, QTuned);
		NodesVector[tally].BestLogLikelihood = BestDensityOfLikelihood;
		NodesVector[tally].PValuesTreatment = PMatrix;
		NodesVector[tally].PValuesMethodology = PMethod;

		//Output information on screen
		cout<<"\t[OK]";
	} 
	catch(...) {
		cout<<"\t[Error]";
	}
	
	cout<<endl;	
	
	return 1;	
}

/*****************************************************
* Name: getPQResults
* Function: Record treatments and methods' results
* Input Parameter: lox, tally, best, AcceptanceRatio, Tuned
* Output: 
* Return Value: int
* Note: lox will be changed
*****************************************************/
int LOX::getPQResults(LOXNode& lox, vector<double> CountByTreatmentOrMethod, vector<DoubleVector> tally, vector<double> best, double AcceptanceRatio, bool Tuned) {
	
	int i, Node;	
	double sum_p = 0.0;

	int num = best.size();

	//Normalization
	for (Node=0; Node < num; Node++) {
		for (i=0; i < tally[Node].size(); i++) sum_p += tally[Node][i];
	}
	for (Node=0; Node < num; Node++) {
		if (CONVERT<int>(CountByTreatmentOrMethod[Node])!=0) {
			for (i=0; i < tally[Node].size(); i++) tally[Node][i] /= sum_p;
		}
		else {
			for (i=0; i < tally[Node].size(); i++) tally[Node][i] = 0;
		}
	}
	
	//Expression levels and confidence intervals
	vector<double> Low, Middle, High;	
		
	for (Node=0; Node < num; Node++) {			
		vector<double> Posterior(tally[Node].begin(), tally[Node].end());
		sort(Posterior.begin(), Posterior.end());
		
		//(-)97.5% CI
		Low.push_back(Posterior[(int)(0.025 * Posterior.size())]);		
		//Median Value 50%
		Middle.push_back(Posterior[(int)(0.5 * Posterior.size())]);
		//(+)97.5% CI
		High.push_back(Posterior[(int)(0.975 * Posterior.size())]);
	}

	//Find the lowest value and fix bug for first column count=0; version 1.5
	double NormalizerToOne = 1;//1 is a big number for this case
	for (Node=0; Node < num; Node++) {
		if (Middle[Node] < NormalizerToOne && CONVERT<int>(CountByTreatmentOrMethod[Node])!=0) NormalizerToOne = Middle[Node];
	}
	
	vector<double> tmp;

	//Expression levels
	for (Node=0; Node < num; Node++) tmp.push_back(Middle[Node] / NormalizerToOne);
	lox.Ratios.push_back(tmp);

	//(-)97.5% Expression levels	
	for (tmp.clear(), Node=0; Node < num; Node++) tmp.push_back((Middle[Node] - Low[Node]) / NormalizerToOne);
	lox.RatiosLowCI.push_back(tmp);

	//(+)97.5% Expression levels	
	for (tmp.clear(), Node=0; Node < num; Node++) tmp.push_back((High[Node] - Middle[Node])/ NormalizerToOne);
	lox.RatiosHighCI.push_back(tmp);

	//Record acceptance ratio
	lox.AcceptanceRatio.push_back(AcceptanceRatio);

	//Record tuned
	lox.Tuned.push_back(Tuned);
	
	//Record ratios by maximum likelihood estimation
	//Bug removing when finding the lowest ratio; version 1.5
	double lowest_best = 1;
	for (Node=0; Node < num; Node++) {
		if (CONVERT<int>(CountByTreatmentOrMethod[Node])!=0) { 
			if (lowest_best > best[Node]) lowest_best = best[Node];
		}
		else {
			best[Node] = 0.0;
		}
	}
	for (tmp.clear(), Node=0; Node < num; Node++) tmp.push_back(best[Node]/lowest_best);
	lox.MLERatios.push_back(tmp);
		
	return 1;
}


/**************************************************
* Function: parseParameter
* Input Parameter: int, const char* []
* Output: Parse the input parameters
* Return Value: bool, success=true, fail=false
***************************************************/
bool LOX::parseParameter(int argc, const char* argv[]) {
	
	bool flag = true;
	string temp = "";
	string msg = "";

	try {
		//Input sequence
		int inputfile_flag=0;
		int i;

		temp=stringtoUpper(argv[1]);
		if (argc==2) {//help information			
			if (temp=="-H") throw 1;
			else throw "no input file";
		}
		else {
			//parse parameters			
			for (i=1; i<argc; i++) {				
				temp = stringtoUpper(argv[i]);
				//Input file
				if (temp=="-I" && (i+1)<argc && inputfile_flag==0) {
					InputFileName = argv[++i];						
					inputfile_flag++;
				}
				//Output methods' results
				else if (temp=="-M" && (i+1)<argc) {					
					int num = CONVERT<int>(argv[++i]);
					if (num==0) isOutputMethodsResults = false;
					else if (num==1) isOutputMethodsResults = true;
					else throw (msg=temp + " " + argv[i]).c_str();
				}				
				//Output MLE ratios
				else if (temp=="-R" && (i+1)<argc) {
					int num = CONVERT<int>(argv[++i]);
					if (num==0) isShowMLERatios = false;
					else if (num==1) isShowMLERatios = true;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Output posterior distribution
				else if (temp=="-G" && (i+1)<argc) {
					int num = CONVERT<int>(argv[++i]);
					if (num==0) isOutputPosterior = false;
					else if (num==1) isOutputPosterior = true;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Sample Size
				else if (temp=="-S" && (i+1)<argc) {
					int num = CONVERT<int>(argv[++i]);
					if (num>0) ChainSampleSize = num;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Burnin
				else if (temp=="-B" && (i+1)<argc) {
					int num = CONVERT<int>(argv[++i]);
					if (num>0) BurnIn = num;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Sampling Period
				else if (temp=="-P" && (i+1)<argc) {
					int num = CONVERT<int>(argv[++i]);
					if (num>0) SamplingPeriod = num;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Tuning Depth
				else if (temp=="-D" && (i+1)<argc) {
					int num = CONVERT<int>(argv[++i]);
					if (num>0) MaxTuningDepth = num;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Lower acceptance rate
				else if (temp=="-L" && (i+1)<argc) {
					double num = CONVERT<double>(argv[++i]);
					if (num>0.0 && num<1.0) LowerAcceptance = num;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				//Upper acceptance rate
				else if (temp=="-U" && (i+1)<argc) {
					double num = CONVERT<double>(argv[++i]);
					if (num>0.0 && num<1.0) UpperAcceptance = num;
					else throw (msg=temp + " " + argv[i]).c_str();
				}
				else {
					throw temp.c_str();
				}				
			}
			//no input file
			if (inputfile_flag==0) throw "no input file";
			if (LowerAcceptance > UpperAcceptance) throw "Lower acceptance rate > higher acceptance rate";
		}		
	}
	catch (const char* e) {
		cout<<ABBRNAME<<": incorrect option, '"<<e<<"'"<<endl;
		cout<<"For more help information, please type '"ABBRNAME<<" -h'."<<endl;
		flag = false;		
	}
	catch(int i) {
		showHelpInfo();
		flag = false;
	}
	catch (...) {		
		cout<<"For more help information, please type '"ABBRNAME<<" -h'."<<endl;
		flag = false;		
	}
	
	return flag;
}

void LOX::showHelpInfo() {
	cout<<"***********************************************************************"<<endl;
	cout<<ABBRNAME<<" ("<<FULLNAME<<"), Version: "<<VERSION<<" ["<<LASTUPDATE<<"]"<<endl;	
	cout<<"Function: "<<FUNCTION<<endl;
	cout<<"Usage: "<<ABBRNAME<<" [OPTIONS]"<<endl;
	cout<<"***********************************************************************"<<endl<<endl;

	cout<<"OPTIONS:"<<endl;	
	cout<<"  -i\tinput file name [string, required]"<<endl;			
	cout<<"  -m\twhether to output methodology results [integer, optional], {0:No || 1:Yes}, default = 0"<<endl;
	cout<<"  -r\twhether to show results by maximum likelihood estimation [integer, optional], {0:No || 1:Yes], default = 0"<<endl;
	cout<<"  -g\twhether to output posterior distribution for each gene [integer, optional], {0:No || 1:Yes], default = 0"<<endl;
	cout<<"  -s\tset chain sample size [integer, optional], default = 1000"<<endl;
	cout<<"  -b\tset burnin size [integer, optional], default = 2000"<<endl;
	cout<<"  -p\tset sampling period [integer, optional], default = 20"<<endl;
	cout<<"  -d\tset tuning depth [integer, optional], default = 8"<<endl;
	cout<<"  -l\tset lower acceptance rate [double, optional], default = 0.15"<<endl;
	cout<<"  -u\tset upper acceptance rate [double, optional], default = 0.5"<<endl;
	cout<<endl;
	
	cout<<"COPYRIGHT & LICENSE:"<<endl;
	cout<<"LOX is distributed as open-source software and licensed under the GNU General Public License "<<endl;
	cout<<"(Version 3; http://www.gnu.org/licenses/gpl.txt), in the hope that it will be useful, but "<<endl;
	cout<<"WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A "<<endl;
	cout<<"PARTICULAR PURPOSE. See the GNU General Public License for more details."<<endl;
	cout<<endl;

	cout<<"REFERENCE:"<<endl;	
	cout<<REFERENCE<<endl;
	cout<<endl;

	cout<<"SEE ALSO:"<<endl;
	cout<<"For more information, please see <http://www.yale.edu/townsend/software.html>."<<endl;
	cout<<endl;

	cout<<"CONTACT:"<<endl;
	cout<<"Please send suggestions or bugs to Zhang Zhang at zhangzhang@big.ac.cn or zhangzhang.cn@gmail.com."<<endl;
	cout<<endl;
}

/*****************************************************
* Name: readDatFile
* Function: Read input file and save as a vector DatLines
* Input Parameter: FileName
* Output: information when reading file
* Return Value: int, success=1, fail=0
* Note: DatLines stores all data
*****************************************************/
int LOX::readDatFile(string FileName) {

	int OKFileName = 0;
	
	try {
		cout<<"Checking input file '"<<FileName.c_str()<<"'..."<<endl;

		ifstream is(trim(FileName).c_str());
		
		if (!is || is.eof())  {
			cout<<endl<<"Error in opening '"<<FileName.c_str()<<"' or no data. Please try again."<<endl;
		}
		else {			
			cout<<"Reading data from the input file..."<<endl;
			while (!is.eof()) {
				string tmp = "";
				getline(is, tmp, '\n');				
				if(tmp!="") DatLines.push_back(tmp);
			}

			//cout<<"Input file '"<<FileName.c_str()<<"' has "<<DatLines.size()<<" lines."<<endl;

			if (DatLines.size() < 4 ) {
				cout<<"Not enough data in '"<<FileName.c_str()<<"'. Please check the input file and try again."<<endl;
			}
			else {
				//cout<<"Input file OK."<<endl;
				OKFileName = 1;
			}
		}
	}
	catch (...) { 
		OKFileName = 0;
	}

	return OKFileName;
}

/*****************************************************
* Name: Run
* Function: calcuate expression ratios for each gene
* Input Parameter: argc, argv
* Output: information for each operation
* Return Value: int
*****************************************************/
int LOX::Run(int argc, const char* argv[]) {
	
	try {
		time_t timeused = time(NULL);

		//Parse input parameters
		if (!parseParameter(argc, argv)) throw "";	
		
		//Output parameter settings
		showParameterSettings();

		//Read input file
		if (!readDatFile(InputFileName)) throw "Error in reading input file.";

		//Set the header from input file
		if (!setHeader()) throw "Error in parsing the header.";
		
		//Set the expression nodes from input file
		if (!setExpressionNodes()) throw "Error in preprocessing expression data.";	

		cout<<endl<<endl<<"Please wait while calculating..."<<endl;
		
		unsigned long num = 0;		
		for(num=0; num < NodesVector.size(); num++) MCMC(num);
		
		cout<<endl<<"Outputing results:"<<endl;
		outputResult();
		outputPValues();

		cout<<endl<<"Mission accomplished. (Time used: "<<formatDateTime(time(NULL) - timeused)<<")"<<endl;
	}
	catch(const char* e) {
		cout<<e<<endl;
	}	
	catch(...) {
		cout<<UNEXPECTED_ERROR<<endl;
	}

	return 1;
}

/*****************************************************
* Name: outputResult
* Function: Output results
* Input Parameter:
* Output: Two files generated: one is for treatments, 
			and the other is for methods.
* Return Value: int
*****************************************************/
int LOX::outputResult() {

	int i, j, k;

	//output columns' titles: 0=treatments, 1=methods
	vector<string> results(2, "");
	
	for (i=0; i<2; i++) {
		results[i] += "Gene";	
		results[i] += "\t";
		results[i] += "Common Name";

		vector<string> tmp;
		if (i==0) tmp = TreatmentsVector;
		else tmp = MethodsVector;

		for(j=0; j<tmp.size(); j++) {
			results[i] += "\t";		results[i] += tmp[j];		
		}

		for(j=0; j<tmp.size(); j++) {
			results[i] += "\t";		results[i] += "(-)97.5%[";
			results[i] += tmp[j];
			results[i] += "]";
		}

		for(j=0; j<tmp.size(); j++) {
			results[i] += "\t";		results[i] += "(+)97.5%[";
			results[i] += tmp[j];
			results[i] += "]";
		}

		results[i] += "\t"; results[i] += "Acceptance Ratio";
		results[i] += "\t"; results[i] += "Density of Best Log-Likelihood";

		if (isShowMLERatios) {
			for(j=0; j<tmp.size(); j++) {
				results[i] += "\t";		results[i] += "MLE(" + tmp[j];		results[i] += ")";
			}
		}

		results[i] += "\n";

		if (!isOutputMethodsResults) break;
	}	
	
	//output results: 0=treatments, 1=methods
	for (i=0; i<2; i++) {	
		for (j=0; j<NodesVector.size(); j++) {
			results[i] += NodesVector[j].ID;
			results[i] += "\t";
			results[i] += NodesVector[j].CommonName;

			if (NodesVector[j].ErrorMsg=="") {
				//Ratio
				for (k=0; k<NodesVector[j].Ratios[i].size(); k++) {
					results[i] += "\t"; results[i] += toString(NodesVector[j].Ratios[i][k]);
				}

				//Lower CI
				for (k=0; k<NodesVector[j].RatiosLowCI[i].size(); k++) {
					results[i] += "\t"; results[i] += toString(NodesVector[j].RatiosLowCI[i][k]);
				}

				//Higher CI
				for (k=0; k<NodesVector[j].RatiosHighCI[i].size(); k++) {
					results[i] += "\t"; results[i] += toString(NodesVector[j].RatiosHighCI[i][k]);
				}

				//Acceptance ratio
				results[i] += "\t"; results[i] += toString(NodesVector[j].AcceptanceRatio[i]);

				//Boolean value
				results[i] += "(";  results[i] += toString(NodesVector[j].Tuned[i]); results[i] += ")";

				//Best Log-likelihood
				results[i] += "\t"; results[i] += toString(NodesVector[j].BestLogLikelihood);

				//MLE Ratio
				for (k=0; k<NodesVector[j].MLERatios[i].size() && isShowMLERatios; k++) {
					results[i] += "\t"; results[i] += toString(NodesVector[j].MLERatios[i][k]);
				}
			}

			results[i] += "\n";
		}
		if (isOutputMethodsResults) continue;
	}

	LOXOutputFileName = InputFileName + ".lox";
	MetOutputFileName = InputFileName + ".met";

	//write treatments' results into file (default)
	writeFile(LOXOutputFileName, results[0].c_str());

	//write methods' results into file (not default)
	if (isOutputMethodsResults) writeFile(MetOutputFileName, results[1].c_str());	

	return 1;
}

/*****************************************************
* Name: outputPValues
* Function: Output PValues
* Input Parameter:
* Output: Each file for every treatment.
* Return Value: int
*****************************************************/
int LOX::outputPValues() {

	int i, Node1, Node2;
	
	//output p-values for treatments
	int NumOfTreatments = TreatmentsVector.size();
	for (Node1 = 0; Node1 < NumOfTreatments; Node1++) {
	
		//output columns' header
		string header = "";	
		header += "Gene";	
		header += "\t";
		header += "Common Name";

		for (Node2 = 0; Node2 < NumOfTreatments; Node2++) {
			if (Node1!=Node2) {
				header += "\t";
				header += "P(" + TreatmentsVector[Node1];
				header += ">" + TreatmentsVector[Node2] +")";
			}
		}

		header += "\n";
		
		//output pvalue results
		string results = "";
		for (i=0; i<NodesVector.size(); i++) {
			results += NodesVector[i].ID;
			results += "\t";
			results += NodesVector[i].CommonName;

			if (NodesVector[i].ErrorMsg=="") {
				for (Node2 = 0; Node2 < NumOfTreatments; Node2++) {
					if (Node1!=Node2) {
						double pvalue = NodesVector[i].PValuesTreatment[Node1*NumOfTreatments + Node2] / (NodesVector[i].PValuesTreatment[Node1*NumOfTreatments + Node2] + NodesVector[i].PValuesTreatment[Node2*NumOfTreatments + Node1]);
						results += "\t";
						results += toString(pvalue);						
					}
				}
			}
			
			results += "\n";
		}
		
		//write treatments' pvalue into file
		string PValueFileName = InputFileName + "." + TreatmentsVector[Node1] + ".pvalue";		
		writeFile(PValueFileName, (header + results).c_str());

	}
	
	//output p-values for methodologies


	int NumOfMethods = MethodsVector.size();
	for (Node1 = 0; Node1 < NumOfMethods; Node1++) {
	
		//output columns' header
		string header = "";	
		header += "Gene";	
		header += "\t";
		header += "Common Name";

		for (Node2 = 0; Node2 < NumOfMethods; Node2++) {
			if (Node1!=Node2) {
				header += "\t";
				header += "P(" + MethodsVector[Node1];
				header += ">" + MethodsVector[Node2] +")";
			}
		}

		header += "\n";
		
		//output pvalue results
		string results = "";
		for (i=0; i<NodesVector.size(); i++) {
			results += NodesVector[i].ID;
			results += "\t";
			results += NodesVector[i].CommonName;

			if (NodesVector[i].ErrorMsg=="") {
				for (Node2 = 0; Node2 < NumOfMethods; Node2++) {
					if (Node1!=Node2) {
						double pvalue = NodesVector[i].PValuesMethodology[Node1*NumOfMethods + Node2] / (NodesVector[i].PValuesMethodology[Node1*NumOfMethods + Node2] + NodesVector[i].PValuesMethodology[Node2*NumOfMethods + Node1]);
						results += "\t";
						results += toString(pvalue);						
					}
				}
			}
			
			results += "\n";
		}
		
		//write treatments' pvalue into file
		string PValueFileName = InputFileName + "." + MethodsVector[Node1] + ".pvalue";		
		writeFile(PValueFileName, (header + results).c_str());

	}

	return 1;
}

/*****************************************************
* Name: main
* Function: Main function to run LOX
* Input Parameter: argc, argv
* Output: Expression levels by LOX
* Return Value: int
*****************************************************/
int main(int argc, const char* argv[]) {
	try {
		LOX zz;
		zz.Run(argc, argv);
	}
	catch (...) {		
	}

	return 1;
}
