/* 
	This file defines the transition matrix an arbitrary time reversible model.
	
	By Sergei L Kosakovsky Pond
*/

function  setElement (h,v,cc)
{	
	mSpecMatrix[h][v]=cc+1;
	mSpecMatrix[v][h]=cc+1;
	return 1;
}

done = 0;

while (!done)
{
	fprintf (stdout,"\nPlease enter a 6 character model designation (e.g:010010 defines HKY85):");
	fscanf  (stdin,"String", modelDesc);
	if (Abs(modelDesc)==6)
	{	
		done = 1;
	}
}

sharedFlag = 0;
categoriesUsed = 0;

#include "modelParameters.mdl";

if (modelType >0)
{
	sharedFlag = 1;
	if (modelType == 2)
	{
		categoriesUsed = 1;
		#include "defineGamma.mdl";
	}
	if (modelType == 3)
	{
		categoriesUsed = 1;
		#include "defineHM.mdl";
	}
}

if (sharedFlag)
{
	global B=1;
	global C=1;
	global D=1;
	global E=1;
	global F=1;
}

a = 1;
b = 2;
cc = 3;
d = 4;
e = 5;
f = 6;
t = 1;


mSpecMatrix = {{*,1,1,1}{1,*,1,1}{1,1,*,1}{1,1,1,*}};

elementAssignments = {{0,0,0,0,0,0}};

paramCount = 1;

r=setElement (0,1,0);

for (j=1; j<6; j=j+1)
{
	for (i=0; i<j; i=i+1)
	{
		if (modelDesc[j]==modelDesc[i])
		{
			elementAssignments[j]=elementAssignments[i];
			if (j<3)
			{
				r=setElement (0,j+1,elementAssignments[i]);
			}
			else
			{
				if (j<5)
				{
					r=setElement (1,j-1,elementAssignments[i]);
				}
				else
				{
					r=setElement (2,3,elementAssignments[i]);
				}
			}
			break;
		}
	}
	if (i==j)
	{
		if (j<3)
		{
			r=setElement (0,j+1,paramCount);
		}
		else
		{
			if (j<5)
			{
				r=setElement (1,j-1,paramCount);
			}
			else
			{
				r=setElement (2,3,paramCount);
			}
		}
		elementAssignments[j] = paramCount;
		paramCount = paramCount+1;
	}
	
}


ChoiceList  (freqResp,"Equilibrium Frequencies",1,NO_SKIP,
			 "Observed","Nucleotide frequencies collected from the data file will be used as equilibrium frequencies.",
			 "Equal","Equal (.25) frequencies are used as equilibrium frequencies",
			 "Custom","Define your own vector of frequencies.");
			 
if (freqResp<0) 
{
	return;
}
else
{
	 if (freqResp<1)
	 {
		HarvestFrequencies (vectorOfFrequencies, filteredData,1,1,0);
		FREQUENCY_SENSITIVE = 1;
	 }
	 else
	 {
	 	if (freqResp<2)
	 	{
	 		equalFreqs = {{.25}{.25}{.25}{.25}};
	 	}
	 	else
	 	{
	 		userFreqs = {4,1};
	 		while (1)
	 		{
	 			fprintf (stdout, "\nEnter the frequencies of A,C and G (T computed automatically). Use spaces to separate the entries:");
	 			fscanf  (stdin, "Number,Number,Number",i,j,r);
	 			if ((i>=0)&&(j>=0)&&(r>=0)&&(i+j+r<=1))
	 			{
	 				userFreqs[0][0]=i;
	 				userFreqs[1][0]=j;
	 				userFreqs[2][0]=r;
	 				userFreqs[3][0]=1-i-j-r;
	 				fprintf (stdout, "\nUsing the following frequencies:",userFreqs);
	 				break;
	 			}
	 			else
	 			{
	 				fprintf (stdout, "\nInvalid choice of frequencies. Please try again");
	 			}
	 		}
	 	}
	 	
	 }
	 
}

function PopulateModelMatrix (ModelMatrixName&, EFV)
{
	ModelMatrixName = {4,4};
	if (sharedFlag)
	{
		t = 1;
		B = 2;
		C = 3;
		D = 4;
		E = 5;
		F = 6;
		if (categoriesUsed)
		{
			for (customLoopCounter=0;customLoopCounter<4;customLoopCounter=customLoopCounter+1)
			{
				for (customLoopCounter_2=customLoopCounter+1;customLoopCounter_2<4;customLoopCounter_2=customLoopCounter_2+1)
				{
					customLoopCounter_r = mSpecMatrix[customLoopCounter][customLoopCounter_2];
					if (customLoopCounter_r==1)
					{
						ModelMatrixName[customLoopCounter][customLoopCounter_2]:=c*t;
						ModelMatrixName[customLoopCounter_2][customLoopCounter]:=c*t;
					}
					else
					{
						if (customLoopCounter_r==2)
						{
							ModelMatrixName[customLoopCounter][customLoopCounter_2]:=c*B*t;
							ModelMatrixName[customLoopCounter_2][customLoopCounter]:=c*B*t;
						}
						else
						{
							if (customLoopCounter_r==3)
							{
								ModelMatrixName[customLoopCounter][customLoopCounter_2]:=c*C*t;
								ModelMatrixName[customLoopCounter_2][customLoopCounter]:=c*C*t;
							}
							else
							{
								if (customLoopCounter_r==4)
								{
									ModelMatrixName[customLoopCounter][customLoopCounter_2]:=c*D*t;
									ModelMatrixName[customLoopCounter_2][customLoopCounter]:=c*D*t;
								}
								else
								{
									if (customLoopCounter_r==5)
									{
										ModelMatrixName[customLoopCounter][customLoopCounter_2]:=c*E*t;
										ModelMatrixName[customLoopCounter_2][customLoopCounter]:=c*E*t;
									}
									else
									{
										ModelMatrixName[customLoopCounter][customLoopCounter_2]:=c*F*t;
										ModelMatrixName[customLoopCounter_2][customLoopCounter]:=c*F*t;
									}
								}
							}
						}
					}
					
				}
			}
		}
		else
		{
			for (customLoopCounter=0;customLoopCounter<4;customLoopCounter=customLoopCounter+1)
			{
				for (customLoopCounter_2=customLoopCounter+1;customLoopCounter_2<4;customLoopCounter_2=customLoopCounter_2+1)
				{
					customLoopCounter_r = mSpecMatrix[customLoopCounter][customLoopCounter_2];
					if (customLoopCounter_r==1)
					{
						ModelMatrixName[customLoopCounter][customLoopCounter_2]:=t;
						ModelMatrixName[customLoopCounter_2][customLoopCounter]:=t;
					}
					else
					{
						if (customLoopCounter_r==2)
						{
							ModelMatrixName[customLoopCounter][customLoopCounter_2]:=B*t;
							ModelMatrixName[customLoopCounter_2][customLoopCounter]:=B*t;
						}
						else
						{
							if (customLoopCounter_r==3)
							{
								ModelMatrixName[customLoopCounter][customLoopCounter_2]:=C*t;
								ModelMatrixName[customLoopCounter_2][customLoopCounter]:=C*t;
							}
							else
							{
								if (customLoopCounter_r==4)
								{
									ModelMatrixName[customLoopCounter][customLoopCounter_2]:=D*t;
									ModelMatrixName[customLoopCounter_2][customLoopCounter]:=D*t;
								}
								else
								{
									if (customLoopCounter_r==5)
									{
										ModelMatrixName[customLoopCounter][customLoopCounter_2]:=E*t;
										ModelMatrixName[customLoopCounter_2][customLoopCounter]:=E*t;
									}
									else
									{
										ModelMatrixName[customLoopCounter][customLoopCounter_2]:=F*t;
										ModelMatrixName[customLoopCounter_2][customLoopCounter]:=F*t;
									}
								}
							}
						}
					}
					
				}
			}
		}
	}
	else
	{
		a = 1;
		b = 2;
		cc = 3;
		d = 4;
		e = 5;
		f = 6;
		for (customLoopCounter=0;customLoopCounter<4;customLoopCounter=customLoopCounter+1)
		{
			for (customLoopCounter_2=customLoopCounter+1;customLoopCounter_2<4;customLoopCounter_2=customLoopCounter_2+1)
			{
				customLoopCounter_r = mSpecMatrix[customLoopCounter][customLoopCounter_2];
				if (customLoopCounter_r==1)
				{
					ModelMatrixName[customLoopCounter][customLoopCounter_2]:=a;
					ModelMatrixName[customLoopCounter_2][customLoopCounter]:=a;
				}
				else
				{
					if (customLoopCounter_r==2)
					{
						ModelMatrixName[customLoopCounter][customLoopCounter_2]:=b;
						ModelMatrixName[customLoopCounter_2][customLoopCounter]:=b;
					}
					else
					{
						if (customLoopCounter_r==3)
						{
							ModelMatrixName[customLoopCounter][customLoopCounter_2]:=cc;
							ModelMatrixName[customLoopCounter_2][customLoopCounter]:=cc;
						}
						else
						{
							if (customLoopCounter_r==4)
							{
								ModelMatrixName[customLoopCounter][customLoopCounter_2]:=d;
								ModelMatrixName[customLoopCounter_2][customLoopCounter]:=d;
							}
							else
							{
								if (customLoopCounter_r==5)
								{
									ModelMatrixName[customLoopCounter][customLoopCounter_2]:=e;
									ModelMatrixName[customLoopCounter_2][customLoopCounter]:=e;
								}
								else
								{
									ModelMatrixName[customLoopCounter][customLoopCounter_2]:=f;
									ModelMatrixName[customLoopCounter_2][customLoopCounter]:=f;
								}
							}
						}
					}
				}
				
			}
		}
	}	
	return 1;
}

custom4x4matrix = 0;

if (freqResp>0)
{
	if (freqResp<2)
	{
		MULTIPLY_BY_FREQS = PopulateModelMatrix ("custom4x4matrix",equalFreqs);
		Model custom4x4 = (custom4x4matrix,equalFreqs,MULTIPLY_BY_FREQS);	
	}
	else
	{
		MULTIPLY_BY_FREQS = PopulateModelMatrix ("custom4x4matrix",userFreqs);
		Model custom4x4 = (custom4x4matrix,userFreqs,MULTIPLY_BY_FREQS);	
		equalFreqs = userFreqs; 
		/* Hack. Used in rebuilding models in NJ and CLUSTER bootstrap */
	}
}
else
{
	MULTIPLY_BY_FREQS = PopulateModelMatrix ("custom4x4matrix",vectorOfFrequencies);
	Model custom4x4 = (custom4x4matrix,vectorOfFrequencies,MULTIPLY_BY_FREQS);	
}