function compareTwoModels (mdl1, mdl2, rt1, rt2, fr1, fr2, useChiBar)
{
	idx1   = mdl1*8+2*rt1+fr1;
	idx2   = mdl2*8+2*rt2+fr2;
	
	if (cacheMatrix[idx1][0]==0)
	{
		res1 = runAModel (modelStrings[mdl1],rt1,fr1,1,idx1);
		cacheMatrix[idx1][0] = res1[1][0];
		cacheMatrix[idx1][1] = res1[1][1]+3*fr1;
	}
	
	if (cacheMatrix[idx2][0]==0)
	{
		res2 = runAModel (modelStrings[mdl2],rt2,fr2,1,idx2);
		cacheMatrix[idx2][0] = res2[1][0];
		cacheMatrix[idx2][1] = res2[1][1]+3*fr2;
	}
	
	LRT	   = 2*(cacheMatrix[idx2][0]-cacheMatrix[idx1][0]);
	DEGF   = cacheMatrix[idx2][1]-cacheMatrix[idx1][1];
	
	if (LRT<=0)
	{
		return 1;					
	}
	if (useChiBar)
	{
		return .5-.5*CChi2(LRT,DEGF);	
	}
	return 1-CChi2(LRT,DEGF);
}

/**************************************************************************/


function echoModelSpec (modelString, modelRate, modelFreq)
{
	fprintf (stdout, "\n\nTo define this model (if it isn't predefined) for analyses,\nUse 'Custom' model with the following options:\n",
					 "Model String:", modelString, "\n");
					 
	fprintf (stdout, "Model Options: ");
	if (modelRate==0)
	{
		fprintf (stdout, "Global");
	}
	else
	{
		if (modelRate==1)
		{
			fprintf (stdout, "Rate Variation, then choose the Gamma distribution");
		}
		else
		{
			if (modelRate==2)
			{
				fprintf (stdout, "Rate Variation, then choose the Gamma+Inv distribution");
			}
			else
			{
				if (modelRate==3)
				{
					fprintf (stdout, "Rate Variation, then choose the Inv Class distribution");
				}
			}
		}
	}
	fprintf (stdout, "\nEquilibrium Frequencies Option: ");
	if (modelFreq)
	{
		fprintf (stdout, "Observed\n\n");	
	}
	else
	{
		fprintf (stdout, "Equal\n\n");	
	}
	return 1;
}   

/**************************************************************************/


function reportTest (idx1, idx2, idx3, idx4, idx5, idx6, pValue, rt1, rt2)
{
	modelName1 = modelNameStrings[idx1][idx5];
	modelName2 = modelNameStrings[idx2][idx6];
	
	if (rt1==1)
	{
		modelName1 = modelName1+"+G";
	}
	if (rt1==2)
	{
		modelName1 = modelName1+"+G+I";
	}
	if (rt1==3)
	{
		modelName1 = modelName1+"+I";
	}
	
	if (rt2==1)
	{
		modelName2 = modelName2+"+G";
	}
	if (rt2==2)
	{
		modelName2 = modelName2+"+G+I";
	}
	if (rt2==3)
	{
		modelName2 = modelName2+"+I";
	}
	
	fprintf (stdout, "\n\tNull:", modelName1, "\n\t\tLog-likelihood = ", Format(cacheMatrix[idx3][0],12,4),
														", ", cacheMatrix[idx3][1], " parameters."); 
														
	fprintf (stdout, "\n\tAlt :", modelName2, "\n\t\tLog-likelihood = ", Format(cacheMatrix[idx4][0],12,4), 
														", ", cacheMatrix[idx4][1], " parameters."); 
														
	fprintf (stdout, "\n\n\tLR statistic       : ", 2*(cacheMatrix[idx4][0]-cacheMatrix[idx3][0]),
					 "\n\tDegrees of freedom : ",cacheMatrix[idx4][1]-cacheMatrix[idx3][1],
					 "\n\tP-Value            : ", pValue);
	
	if (pValue<rejectAt)
	{
		fprintf (stdout, "\n\n\tNull hypothesis rejected\n\t", modelName2, " chosen.\n");
		return modelName2;
	}
	else
	{
		fprintf (stdout, "\n\n\tNull hypothesis accepted\n\t", modelName1, " chosen.\n");
		return modelName1;
	}	
	return 1;
}                         

/**************************************************************************/

function doHierarchicalTest (echoFlag)
{
	if (echoFlag)
	{
		fprintf (stdout, "\n\n****** RUNNING HIERARCHICAL MODEL TESTING ******\n\n");
	}

	chosenMatrix = 0;
	chosenRates  = 0;
	chosenFreqs  = 0;
	chosen2		 = 0;
	chosen3      = 0;
	chosen4		 = 0;
	
	/* do the frequencies test */
	
	if (echoFlag)
	{
		fprintf (stdout,"\n1). Checking for equilibrium frequencies equality.\n");
	}
	pVal =  compareTwoModels (0,0,0,0,0,1,0);

	chosenFreqs = (pVal<rejectAt);

	if (echoFlag)
	{
		r = reportTest (0, 0, 0, 1, 0, 1, pVal,0,0);
	}
	
	/* do the tr=tv test */
	
	if (echoFlag)
	{
		fprintf (stdout,"\n2). Checking for equality of transition and transversion rates.\n");
	}
	
	pVal =  compareTwoModels (0,1,0,0,chosenFreqs,chosenFreqs,0);
	
	chosen2 = (pVal<rejectAt);
	
	if (echoFlag)
	{
		r = reportTest (0, 1, chosenFreqs, 8+chosenFreqs, chosenFreqs, chosenFreqs, pVal,0,0);
	}
	
	if (chosen2)
	/* fine tune the matrix */
	{
		if (echoFlag)
		{
			fprintf (stdout,"\n3). Checking for equality of A<->G and C<->T (transition) rates.\n");
		}
		pVal =  compareTwoModels (1,2,0,0,chosenFreqs,chosenFreqs,0);
		chosen3 = (pVal<rejectAt);
		
		if (echoFlag)
		{
			r = reportTest (1, 2, 8+chosenFreqs, 8*2+chosenFreqs, chosenFreqs, chosenFreqs, pVal,0,0);
		}
		
		if (echoFlag)
		{
			fprintf (stdout,"\n4). Checking whether there are 1 or 2 transversion rates.\n");
		}
		
		pVal =  compareTwoModels (1+chosen3,3+chosen3,0,0,chosenFreqs,chosenFreqs,0);
		
		chosen4 = (pVal<rejectAt);
		chosenMatrix = 1+chosen3+2*chosen4;
		if (echoFlag)
		{
			r = reportTest (1+chosen3, 3+chosen3, 8(1+chosen3)+chosenFreqs, 8(3+chosen3)+chosenFreqs, chosenFreqs, chosenFreqs, pVal,0,0);
		}
		if (chosen4)
		/* fine tune some more */
		{
			if (echoFlag)
			{
				fprintf (stdout,"\n5). Checking whether there are 2 or 4 transversion rates.\n");
			}
			pVal =  compareTwoModels (3+chosen3,5+chosen3,0,0,chosenFreqs,chosenFreqs,0);
			if (echoFlag)
			{
				r = reportTest (3+chosen3, 5+chosen3, 8(3+chosen3)+chosenFreqs, 8(5+chosen3)+chosenFreqs, chosenFreqs, chosenFreqs, pVal,0,0);
			}
			chosenMatrix = 3+chosen3+2*(pVal<rejectAt);
		}
		else
		{
			if (echoFlag)
			{
				fprintf (stdout,"\n...Skipping step 5...\n");
			}
		}
	}
	else
	{
		if (echoFlag)
		{
			fprintf (stdout,"\n...Skipping steps 3 through 6...\n");
		}
		chosenMatrix = 0;
	}
	
	if (echoFlag)
	{
		fprintf (stdout,"\n6). Checking for evidence of rate variation.\n");
	}
	pVal =  compareTwoModels (chosenMatrix,chosenMatrix,0,1,chosenFreqs,chosenFreqs,1);
	
	if (echoFlag)
	{
		r = reportTest (chosenMatrix, chosenMatrix, 8*chosenMatrix+chosenFreqs, 8*chosenMatrix+chosenFreqs+2, chosenFreqs, chosenFreqs, pVal,0,1);
		fprintf (stdout,"\n7). Checking for evidence of an invariant rate class.\n");
	}
	if (pVal<rejectAt)
	/* accept rate variation */
	{
		pVal =  compareTwoModels (chosenMatrix,chosenMatrix,1,2,chosenFreqs,chosenFreqs,1);
		if (echoFlag)
		{
			r = reportTest (chosenMatrix, chosenMatrix, 8*chosenMatrix+chosenFreqs+2, 8*chosenMatrix+chosenFreqs+4, chosenFreqs, chosenFreqs, pVal,1,2);
		}
		if (pVal<rejectAt)
		{
			chosenRates = 2;
		}
		else
		{
			chosenRates = 1;
		}
	}
	else
	{
		pVal =  compareTwoModels (chosenMatrix,chosenMatrix,0,3,chosenFreqs,chosenFreqs,1);
		if (echoFlag)
		{
			r = reportTest (chosenMatrix, chosenMatrix, 8*chosenMatrix+chosenFreqs, 8*chosenMatrix+chosenFreqs+6, chosenFreqs, chosenFreqs, pVal,0,3);
		}
		if (pVal<rejectAt)
		{
			chosenRates = 3;
		}
		else
		{
			chosenRates = 0;
		}	
	}
	return r;
}

/**************************************************************************/

function doAICTest (echoFlag)
{
	if (echoFlag)
	{
		fprintf (stdout, "\n\n****** RUNNING AIC BASED MODEL SELECTION ******\n\n");
	}
	aicRates = 0;
	aicFreqs = 0;
	aicType  = 0;
	aicScore = 1e100;
	cModelN  = "";
	
	if (echoFlag)
	{	
		fprintf (stdout,"\n|  Model     | # prm |    lnL    |    AIC     |");   
		fprintf (stdout,"\n|------------|-------|-----------|------------|");   
	}		
	for (k=0; k<Columns(modelStrings); k=k+1)
	{
		for (i=0; i<2; i=i+1)
		{
			for (j=0; j<4; j=j+1)
			{
				thisModelName = modelNameStrings[k][i];
				if (j==1)
				{
					thisModelName = thisModelName+"+G";
				}
				if (j==2)
				{
					thisModelName = thisModelName+"+G+I";
				}
				if (j==3)
				{
					thisModelName = thisModelName+"+I";
				}
				while (Abs(thisModelName)<11)
				{
					thisModelName = thisModelName+" ";						
				}
				idx1   = k*8+2*j+i;
				res = runAModel (modelStrings[k],j,i,1,idx1);
				if (i)
				{
					res[1][1] = res[1][1]+3;
				}
				modelAIC = 2*(res[1][1]-res[1][0]);
				
				cacheMatrix[idx1][0] = res[1][0];
				cacheMatrix[idx1][1] = res[1][1];
				
				if (modelAIC < aicScore)
				{
					aicFreqs = i;
					aicRates = j;
					aicType  = k;
					aicScore = modelAIC;
					cModelN  = thisModelName;
				}
				if (echoFlag)
				{
					fprintf (stdout, "\n| ", thisModelName, "| ", Format (res[1][1],5,0),
										   " | ", Format (res[1][0],9,3), " | ", Format (modelAIC,10,3), " |");
				}
			}
		}
	}
	if (echoFlag)
	{
		fprintf (stdout,"\n|------------|-------|-----------|------------|");   
	}
	
	return cModelN;
}

/**************************************************************************/

function printModelMatrix (modelString)
{
	mstrConv = "1 ";
	for (v2 = 1; v2 < 6; v2 = v2+1)
	{
		if (modelString[v2]=="0")
		{
			mstrConv = mstrConv+"1 ";
		}
		else
		{
			if (modelString[v2]=="1")
			{
				mstrConv = mstrConv+"R1";
			}
			else
			{
				if (modelString[v2]=="2")
				{
					mstrConv = mstrConv+"R2";
				}
				else
				{
					if (modelString[v2]=="3")
					{
						mstrConv = mstrConv+"R3";
					}
					else
					{
						if (modelString[v2]=="4")
						{
							mstrConv = mstrConv+"R4";
						}
						else
						{
							mstrConv = mstrConv+"R5";
						}
					}
				}
			}
		}
	}
	sep = "+---+-------+-------+-------+-------+\n";
	fprintf (stdout, sep,
					 "|   |   A   |   C   |   G   |   T   |\n",
					 sep,
					 "| A |   *   | ", mstrConv[0][1], "*t  | ", mstrConv[2][3], "*t  | ", mstrConv[4][5], "*t  |\n",
					 sep,
					 "| C | ", mstrConv[0][1], "*t  |   *   | ", mstrConv[6][7], "*t  | ", mstrConv[8][9], "*t  |\n",
					 sep,
					 "| G | ", mstrConv[2][3], "*t  | " , mstrConv[6][7], "*t  |   *   | ", mstrConv[10][11], "*t  |\n",
					 sep,
					 "| T | ", mstrConv[4][5], "*t  | " , mstrConv[8][9], "*t  | ", mstrConv[10][11], "*t  |   *   |\n",
					 sep);
						 
	
	return 1;
}

/**************************************************************************/

function  setElement (h,v,cell)
{
	if (modelType==0)
	{
		if (cell==0)
		{
			m[h][v]:=t;
			m[v][h]:=t;
		}
		else
		{
			if (cell==1)
			{
				m[h][v]:=R1*t;
				m[v][h]:=R1*t;
			}
			else
			{
				if (cell==2)
				{
					m[h][v]:=R2*t;
					m[v][h]:=R2*t;
				}
				else
				{
					if (cell==3)
					{
						m[h][v]:=R3*t;
						m[v][h]:=R3*t;
					}
					else
					{
						if (cell==4)
						{
							m[h][v]:=R4*t;
							m[v][h]:=R4*t;
						}
						else
						{
							m[h][v]:=R5*t;
							m[v][h]:=R5*t;
						}
					}
				}
			}
		}
	}
	else
	{
		if (cell==0)
		{
			m[h][v]:=c*t;
			m[v][h]:=c*t;
		}
		else
		{
			if (cell==1)
			{
				m[h][v]:=c*R1*t;
				m[v][h]:=c*R1*t;
			}
			else
			{
				if (cell==2)
				{
					m[h][v]:=c*R2*t;
					m[v][h]:=c*R2*t;
				}
				else
				{
					if (cell==3)
					{
						m[h][v]:=c*R3*t;
						m[v][h]:=c*R3*t;
					}
					else
					{
						if (cell==4)
						{
							m[h][v]:=c*R4*t;
							m[v][h]:=c*R4*t;
						}
						else
						{
							m[h][v]:=c*R5*t;
							m[v][h]:=c*R5*t;
						}
					}
				}
			}
		}
	}
	return 1;
}

/**************************************************************************/

function runAModel (modelSpecString, gammaType, freqType, doRun, mIdx)
{
	modelType = (gammaType>0);
	/* define the rate matrix */
	r = setElement  (0,1,modelSpecString[0]);
	r = setElement  (0,2,modelSpecString[1]);
	r = setElement  (0,3,modelSpecString[2]);
	r = setElement  (1,2,modelSpecString[3]);
	r = setElement  (1,3,modelSpecString[4]);
	r = setElement  (2,3,modelSpecString[5]);
	/* define rate variation */
	if (gammaType==1)
	{
		global alpha = .5;
		alpha:>0.01;alpha:<100;
		category c = (classCount, EQUAL, MEAN, 
						GammaDist(_x_,alpha,alpha), 
						CGammaDist(_x_,alpha,alpha), 
						0 , 
				  	    1e25,
				  	    CGammaDist(_x_,alpha+1,alpha)
				  	 );
	}
	else
	{
		if (gammaType==2)
		{
			global alpha   = .5;
			alpha:>0.01;alpha:<100;

			global P;
			P :< 1-(1e-10);
			P = 1/(classCount+1);
			
			
			catFreqMatrix = {classCount+1,1};
			catFreqMatrix [0] := P;
			
			for (h=1; h<=classCount;h=h+1)
			{
				catFreqMatrix [h] := (1-P)/classCount__;
			}
			
			category c = (classCount+1, catFreqMatrix, MEAN, 
							(1-P)*GammaDist(_x_,alpha,alpha)*(_x_>0), 
							(1-P)*CGammaDist(_x_,alpha,alpha)*(_x_>0)+P, 
							0, 
					  	    1e25,
					  	    (1-P)*CGammaDist(_x_,alpha+1,alpha)*(_x_>0)
					  	 );
		}	
		else
		{
			if (gammaType==3)
			{
				global P;
				P :< 1-(1e-10);
				P = 1/2;
				
				catFreqMatrix = {2,1};
				catFreqMatrix [0] := P;
				catFreqMatrix [1] := (1-P);
				
				catRateMatrix = {{0,1}};
				
				category c = (2, catFreqMatrix, MEAN, 
								, 
								catRateMatrix, 
								0, 
						  	    1e25,
						  	 );
			}	
		}
	}
	if (freqType)
	{
		Model testingModel = (m,PI_EFV);
	}
	else
	{
		Model testingModel = (m,quarters_EFV);
	}
	
	if (doRun)
	{
		Tree MTTree = treeString;
		LikelihoodFunction lf = (filteredData, MTTree);
		Optimize (res,lf);
		if (Rows (savedMLEs)==1)
		{
			savedMLEs     = {100,Columns (res)+10};
			savedParNames = {100,Columns (res)+10};
			
			for (kk=0; kk<100; kk=kk+1)
			{
				for (jj=0; jj<Columns (savedMLEs); jj=jj+1)
				{
					savedMLEs[kk][jj] = -1;
				}
			}
		}
		for (kk=0; kk<Columns (res); kk=kk+1)
		{
			savedMLEs[mIdx][kk] = res[0][kk];
			GetString (parName, lf, kk);
			savedParNames[mIdx][kk] = parName;
		}
		return   res;
	}
	return 0;
}

/**************************************************************************/

OPTIMIZE_SUMMATION_ORDER = 1;

m={4,4};
savedMLEs = {1,1};

global R1 = 1;
global R2 = 1;
global R3 = 1;
global R4 = 1;
global R5 = 1;

fprintf (stdout, "\n+--------------------------------+\n| RUNNING MODEL TESTING ANALYSIS |\n| Based on the program ModelTest |\n|              by                |\n|        David  Posada           |\n|             and                |\n|        Keith Crandall          |\n|                                |\n|    If you use this analysis,   |\n| be sure to cite the original   |\n| reference, which can be found  |\n| in the HyPhy documentation     |\n| page for this analysis.        |\n+--------------------------------+\n");

SetDialogPrompt   ("Please specify a nucleotide data file:");

DataSet ds    = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter filteredData = CreateFilter (ds,1);

fprintf (stdout, "\n\nData read from \"",LAST_FILE_PATH,"\"\n",ds,"\n");
HarvestFrequencies (PI_EFV,filteredData,1,1,1);

quarters_EFV = {{.25}{.25}{.25}{.25}};

_DO_TREE_REBALANCE_ = 1;

#include "queryTree.bf";

classCount = 0;

while (classCount<=0)
{
	fprintf (stdout, "\nNumber of rate classes in rate variation models (e.g. 4):");
	fscanf  (stdin,"Number", classCount);
}

KEEP_OPTIMAL_ORDER = 1;
MESSAGE_LOGGING = 0;

modelNum	= 0;
rejectCount = 0;


cacheMatrix     = {100,2};

modelStrings 	= {{"000000", /*F81,JC69*/
                 	"010010", /*HKY85, K80*/
                 	"010020", /*TrN*/
                 	"012210", /*K81*/
                 	"012230", /*TIM*/
                 	"012310", /*TVM*/
                	"012345" /*REV*/
                }};
                
 modelNameStrings = {{"JC69",
                      "F81"},
   				     {"K80",
   				      "HKY85"},
   				     {"TrNef",
   				      "TrN"},
   				     {"K81",
   				      "K81uf"},
   				     {"TIMef",
   				      "TIM"},
   				     {"TVMef",
   				      "TVM"},
   				     {"SYM",
   				      "GTR"}};
   				        				     
ChoiceList (runType,"Model Selection Method",1,SKIP_NONE,
			"Hierarchical Test","Perform a series of nested model comparisons to select the model.",
			"AIC Test","Obtain MLEs for each model, and select the best one using Akaike Information Criterion.",
			"Both","Run the AIC test, followed by the Hierarchical test");
			
if (runType<0)
{
	return;
}

if (runType!=1)
{
	rejectAt = 0;

	while ((rejectAt<=0)||(rejectAt>=1))
	{
		fprintf (stdout, "\nModel rejection level (e.g. 0.05):");
		fscanf  (stdin,"Number", rejectAt);
	}
}

if (runType)
{
	aicMdl = doAICTest (1);
	fprintf (stdout, "\n\nAIC based model: ",aicMdl,", AIC = ", aicScore);
}

if (runType!=1)
{
	hMdl = doHierarchicalTest (1);
}

OPTIMIZE_SUMMATION_ORDER = 0;

if (runType!=1)
{
	res = runAModel (modelStrings[chosenMatrix],chosenRates,chosenFreqs,0,0);
	Tree MTTree = treeString;
	LikelihoodFunction AIC_LF = (filteredData, MTTree);
	idx1 = 8*chosenMatrix+2*chosenRates+chosenFreqs;
	
	for (k=0; k<Columns (savedMLEs); k=k+1)
	{
		j = savedMLEs[idx1][k];
		if (j<0)
		{
			break;
		}
		
		GetString (parName, AIC_LF, k);
		
		for (i = 0; i< Columns (savedMLEs); i=i+1)
		{
			if (parName == savedParNames[idx1][i])
			{
				break;
			}
		}
		
		if (i==Columns (savedMLEs))
		{
			fprintf (stdout, "\nInternal Error in ModelTest.bf\n");
		}
		
		j = savedMLEs[idx1][i];
		SetParameter (AIC_LF,k,j);
	}
	
	
	fprintf (stdout, "\n\n ******* Hierarchical Testing based model (",hMdl,") ******** \n\nRate Matrix:\n");
	r = printModelMatrix (modelStrings[chosenMatrix]);
	r = echoModelSpec	 (modelStrings[chosenMatrix], chosenRates, chosenFreqs);
	fprintf (stdout, "\n\nMLEs for the model selected by Hierarchical Testing:\n", AIC_LF, "\n");
}

if (runType)
{
	res = runAModel (modelStrings[aicType],aicRates,aicFreqs,0,0);
	Tree MTTree = treeString;
	LikelihoodFunction AIC_LF = (filteredData, MTTree);

	for (k=0; k<Columns (savedMLEs); k=k+1)
	{
		j = savedMLEs[idx1][k];
		if (j<0)
		{
			break;
		}
		
		GetString (parName, AIC_LF, k);
		
		for (i = 0; i< Columns (savedMLEs); i=i+1)
		{
			if (parName == savedParNames[idx1][i])
			{
				break;
			}
		}
		
		if (i==Columns (savedMLEs))
		{
			fprintf (stdout, "\nInternal Error in ModelTest.bf\n");
		}
		
		j = savedMLEs[idx1][i];
		SetParameter (AIC_LF,k,j);
	}

	fprintf (stdout, "\n\n ******* AIC based model (", aicMdl,") ******** \n\nRate Matrix:\n");
	r = printModelMatrix (modelStrings[aicType]);
	r = echoModelSpec	 (modelStrings[aicType], aicRates, aicFreqs);
	fprintf (stdout, "\n\nMLEs for the model selected by AIC:\n", AIC_LF, "\n");
}

OPTIMIZE_SUMMATION_ORDER = 1;