_S_NS_POSITIONS_ = {2,64};

hShift = 0;

for (h=0; h<64; h=h+1)
{
	myAA = _Genetic_Code[h];
	
	if (myAA != 10) /* not a stop codon */
	{
		norm_factor  = 0.0;
		sSites	 	 = 0.0;
		nsSites 	 = 0.0;
		
		/* first position change */
		/* actual first position */
		
		p1 = h$16; /* 0->A, 1->C, 2->G, 3->T */
		p2 = h%16; /* remainder - i.e. positions 2 and 3*/
				
		for (n1 = 0; n1 < 4; n1=n1+1)
		{
			if (n1 != p1) /* a change */
			{
				/* new codon */
				nc = n1*16 + p2;
				newAA = _Genetic_Code[nc];
				
				if (newAA!=10) /* not a stop codon */
				{
					if (newAA == myAA) /* syn. change */
					{
						sSites = sSites + _EFV_MATRIX0_[p1][n1];
					}
					else
					{
						nsSites = nsSites + _EFV_MATRIX0_[p1][n1];
					}
				}
				norm_factor = norm_factor + _EFV_MATRIX0_[p1][n1];
			}
		}
		
		if (norm_factor)
		{
			_S_NS_POSITIONS_[0][h] = _S_NS_POSITIONS_[0][h] + sSites/norm_factor;
			_S_NS_POSITIONS_[1][h] = _S_NS_POSITIONS_[1][h] + nsSites/norm_factor;
		}
		
		norm_factor  = 0.0;
		sSites	 	 = 0.0;
		nsSites 	 = 0.0;

		/* second position change */
		/* actual second position */
		
		p1 = (h%16)$4;
		p2 = (h$16)*16+h%4; /* remainder - i.e. positions 1 and 3*/
		
		for (n1 = 0; n1 < 4; n1=n1+1)
		{
			if (n1 != p1) /* a change */
			{
				/* new codon */
				nc = n1*4 + p2;
				newAA = _Genetic_Code[nc];
				
				if (newAA!=10) /* not a stop codon */
				{
					if (newAA == myAA) /* syn. change */
					{
						sSites = sSites + _EFV_MATRIX1_[p1][n1];
					}
					else
					{
						nsSites = nsSites + _EFV_MATRIX1_[p1][n1];
					}
				}
				norm_factor = norm_factor + _EFV_MATRIX1_[p1][n1];
			}
		}

		/* 3rd position change */
		/* actual 3rd position */
		
		if (norm_factor)
		{
			_S_NS_POSITIONS_[0][h] = _S_NS_POSITIONS_[0][h] + sSites/norm_factor;
			_S_NS_POSITIONS_[1][h] = _S_NS_POSITIONS_[1][h] + nsSites/norm_factor;
		}
		
		norm_factor  = 0.0;
		sSites	 	 = 0.0;
		nsSites 	 = 0.0;

		p1 = h%4;
		p2 = (h$4)*4; /* remainder - i.e. positions 1 and 2*/
		
		for (n1 = 0; n1 < 4; n1=n1+1)
		{
			if (n1 != p1) /* a change */
			{
				/* new codon */
				nc = n1 + p2;
				newAA = _Genetic_Code[nc];
				
				if (newAA!=10) /* not a stop codon */
				{
					if (newAA == myAA) /* syn. change */
					{
						sSites = sSites + _EFV_MATRIX2_[p1][n1];
					}
					else
					{
						nsSites = nsSites + _EFV_MATRIX2_[p1][n1];
					}
				}
				norm_factor = norm_factor + _EFV_MATRIX2_[p1][n1];
			}
		}
		
		if (norm_factor)
		{
			_S_NS_POSITIONS_[0][h] = _S_NS_POSITIONS_[0][h] + sSites/norm_factor;
			_S_NS_POSITIONS_[1][h] = _S_NS_POSITIONS_[1][h] + nsSites/norm_factor;
		}
		
	}
	else
	{
		hShift = hShift+1;
	}
}

stateCharCount = 64-hShift;


/* now compute pairwise codon distance matrix */

nuc_split 	  = {2,3};
_PAIRWISE_S_  = {stateCharCount,stateCharCount};
_PAIRWISE_NS_ = {stateCharCount,stateCharCount};
_OBSERVED_S_  = {stateCharCount,stateCharCount};
_OBSERVED_NS_ = {stateCharCount,stateCharCount};
_NUC_SUB_TYPE_= {stateCharCount,stateCharCount};

_NUC_TEMPL_	  = {{0,0,1,2}
                 {0,0,3,4}
                 {1,3,0,5}
                 {2,4,5,0}};

hShift = 0;

for (h = 0; h<64; h=h+1)
{
	myAA = _Genetic_Code[h];
	if (myAA == 10)
	{
		hShift = hShift+1;
	}
	else
	{
		nuc_split[0][0] = h$16;
		nuc_split[0][1] = (h%16)$4;
		nuc_split[0][2] = h%4;
		
		_PAIRWISE_S_  [h-hShift][h-hShift] = _S_NS_POSITIONS_[0][h];
		_PAIRWISE_NS_ [h-hShift][h-hShift] = _S_NS_POSITIONS_[1][h];
		
		vShift = hShift;
		
		for (v = h+1; v<64; v=v+1)
		{
			newAA = _Genetic_Code[v];
			if (newAA == 10)
			{
				vShift = vShift+1;
			}
			else
			{
				/*fprintf ("echo", codonString (h), "->", codonString(v), "\n");*/
				nuc_split[1][0] = v$16;
				nuc_split[1][1] = (v%16)$4;
				nuc_split[1][2] = v%4;
				
				p1 = (nuc_split[1][0]!=nuc_split[0][0])+
					 (nuc_split[1][1]!=nuc_split[0][1])+
					 (nuc_split[1][2]!=nuc_split[0][2]);
					 
				if (p1==1) 
				{
					_PAIRWISE_S_ [h-hShift][v-vShift] = (_S_NS_POSITIONS_[0][h]+_S_NS_POSITIONS_[0][v])/2;
					_PAIRWISE_NS_ [h-hShift][v-vShift] = (_S_NS_POSITIONS_[1][h]+_S_NS_POSITIONS_[1][v])/2;
					
					if (myAA == newAA)
					{
						_OBSERVED_S_[h-hShift][v-vShift] = 1;
					}
					else
					{
						_OBSERVED_NS_[h-hShift][v-vShift] = 1;					
					}
					
					/*fprintf ("echo", "\tOne change:", _PAIRWISE_S_ [h][v], "\t", _PAIRWISE_NS_ [h][v], "\n");*/
					
					for (p1 = 0; p1 < 3; p1=p1+1)
					{
						if (nuc_split[0][p1] != nuc_split[1][p1])
						{
							_NUC_SUB_TYPE_[h-hShift][v-vShift] = _NUC_TEMPL_[nuc_split[0][p1]][nuc_split[1][p1]];
							break;
						}
					}
				}
				else
				{
					if (p1==2) 
					{
						/*fprintf ("echo", "\tTwo changes:\n");*/
						if (nuc_split[1][0]==nuc_split[0][0]) 
						{
							pc1 = 16*nuc_split[0][0]+4*nuc_split[0][1]+nuc_split[1][2];
							pc2 = 16*nuc_split[0][0]+4*nuc_split[1][1]+nuc_split[0][2];
						}
						else
						{
							if (nuc_split[1][1]==nuc_split[0][1])
							{
								pc1 = 16*nuc_split[0][0]+4*nuc_split[0][1]+nuc_split[1][2];
								pc2 = 16*nuc_split[1][0]+4*nuc_split[0][1]+nuc_split[0][2];
							}
							else 
							{
								pc1 = 16*nuc_split[1][0]+4*nuc_split[0][1]+nuc_split[0][2];
								pc2 = 16*nuc_split[0][0]+4*nuc_split[1][1]+nuc_split[0][2];
							}
						}
						
						pc = 0;
						
						pc1AA = _Genetic_Code[pc1];
						if (pc1AA != 10)
						{	
							_OBSERVED_S_   [h-hShift][v-vShift]  = (pc1AA == myAA) + (pc1AA == newAA);
							_OBSERVED_NS_  [h-hShift][v-vShift]  = (pc1AA != myAA) + (pc1AA != newAA);
							
							_PAIRWISE_S_  [h-hShift][v-vShift] = (_S_NS_POSITIONS_[0][h]+_S_NS_POSITIONS_[0][v]+_S_NS_POSITIONS_[0][pc1])/3;
							_PAIRWISE_NS_ [h-hShift][v-vShift] = (_S_NS_POSITIONS_[1][h]+_S_NS_POSITIONS_[1][v]+_S_NS_POSITIONS_[1][pc1])/3;
							/*fprintf ("echo", "\t", codonString (h), "->", codonString(pc1), "->", codonString(v), "\t", _PAIRWISE_S_ [h][v], "\t", _PAIRWISE_NS_ [h][v], "\n");*/
							pc = 1;
						}
						
						pc1AA = _Genetic_Code[pc2];
						if (pc1AA != 10)
						{
							_OBSERVED_S_  [h-hShift][v-vShift]  = _OBSERVED_S_ [h-hShift][v-vShift] + (pc1AA == myAA) + (pc1AA == newAA);
							_OBSERVED_NS_ [h-hShift][v-vShift]  = _OBSERVED_NS_[h-hShift][v-vShift] + (pc1AA != myAA) + (pc1AA != newAA);

							_PAIRWISE_S_  [h-hShift][v-vShift] = _PAIRWISE_S_ [h-hShift][v-vShift]+(_S_NS_POSITIONS_[0][h]+_S_NS_POSITIONS_[0][v]+_S_NS_POSITIONS_[0][pc2])/3;
							_PAIRWISE_NS_ [h-hShift][v-vShift] = _PAIRWISE_NS_ [h-hShift][v-vShift]+(_S_NS_POSITIONS_[1][h]+_S_NS_POSITIONS_[1][v]+_S_NS_POSITIONS_[1][pc2])/3;
							/*fprintf ("echo", "\t", codonString (h), "->", codonString(pc2), "->", codonString(v), "\t", _PAIRWISE_S_ [h][v], "\t", _PAIRWISE_NS_ [h][v], "\n");*/
							pc = pc+1;
						}
						
						if (pc == 0)
						{
							_PAIRWISE_S_  [h-hShift][v-vShift] = 0;
							_PAIRWISE_NS_ [h-hShift][v-vShift] = 0;
							_OBSERVED_S_  [h-hShift][v-vShift] = 0;
							_OBSERVED_NS_ [h-hShift][v-vShift] = 0;
						}
						else
						{
							_PAIRWISE_S_  [h-hShift][v-vShift] = _PAIRWISE_S_  [h-hShift][v-vShift]/pc;
							_PAIRWISE_NS_ [h-hShift][v-vShift] = _PAIRWISE_NS_ [h-hShift][v-vShift]/pc;
							_OBSERVED_S_  [h-hShift][v-vShift] = _OBSERVED_S_  [h-hShift][v-vShift]/pc;
							_OBSERVED_NS_ [h-hShift][v-vShift] = _OBSERVED_NS_ [h-hShift][v-vShift]/pc;						
						}
					}
					else 
					{
						pc = 0;
						
						/*
						fprintf ("echo", "\tThree changes:\n");
						*/
						
						for (p1=0;p1<3;p1=p1+1)
						{
							pc1   = (h+4^(2-p1)*(nuc_split[1][p1]-nuc_split[0][p1])+0.5)$1;
							pc1AA = _Genetic_Code[pc1];
							if (pc1AA != 10)
							{
								for (p2=0;p2<3;p2=p2+1)
								{
									if (p2 != p1)
									{
										pc2 = (pc1+4^(2-p2)*(nuc_split[1][p2]-nuc_split[0][p2]) + 0.5)$1;
										pc2AA = _Genetic_Code[pc2];
										if (pc2AA != 10)
										{
											_OBSERVED_S_   [h-hShift][v-vShift]  = _OBSERVED_S_ [h-hShift][v-vShift] + (pc1AA == myAA) + (pc1AA == pc2AA) + (pc2AA == newAA);
											_OBSERVED_NS_  [h-hShift][v-vShift]  = _OBSERVED_NS_[h-hShift][v-vShift] + (pc1AA != myAA) + (pc1AA != pc2AA) + (pc2AA != newAA);
											
											_PAIRWISE_S_ [h-hShift][v-vShift] = _PAIRWISE_S_ [h-hShift][v-vShift]+
																(_S_NS_POSITIONS_[0][h]+_S_NS_POSITIONS_[0][v]+_S_NS_POSITIONS_[0][pc1]+_S_NS_POSITIONS_[0][pc2])/4;
											_PAIRWISE_NS_ [h-hShift][v-vShift] = _PAIRWISE_NS_ [h-hShift][v-vShift]+
																(_S_NS_POSITIONS_[1][h]+_S_NS_POSITIONS_[1][v]+_S_NS_POSITIONS_[1][pc1]+_S_NS_POSITIONS_[1][pc2])/4;
											pc = pc+1;
										
											/*
											fprintf (stdout, "\t", codonString (h), "->", codonString(pc1), "->", codonString(pc2), "->",codonString(v), "\t", _PAIRWISE_S_ [h][v], "\t", _PAIRWISE_NS_ [h][v], "\n");
											*/
										}
									}
								}
							}
						}
						if (pc == 0)
						{
							_PAIRWISE_S_  [h-hShift][v-vShift] = 0;
							_PAIRWISE_NS_ [h-hShift][v-vShift] = 0;
							_OBSERVED_S_  [h-hShift][v-vShift] = 0;
							_OBSERVED_NS_ [h-hShift][v-vShift] = 0;
						}
						else
						{
							_PAIRWISE_S_  [h-hShift][v-vShift] = _PAIRWISE_S_  [h-hShift][v-vShift]/pc;
							_PAIRWISE_NS_ [h-hShift][v-vShift] = _PAIRWISE_NS_ [h-hShift][v-vShift]/pc;
							_OBSERVED_S_  [h-hShift][v-vShift] = _OBSERVED_S_  [h-hShift][v-vShift]/pc;
							_OBSERVED_NS_ [h-hShift][v-vShift] = _OBSERVED_NS_ [h-hShift][v-vShift]/pc;						
						}
					}
				}
				_PAIRWISE_S_  [v-vShift][h-hShift] = _PAIRWISE_S_  [h-hShift][v-vShift];
				_PAIRWISE_NS_ [v-vShift][h-hShift] = _PAIRWISE_NS_ [h-hShift][v-vShift];
				_OBSERVED_S_  [v-vShift][h-hShift] = _OBSERVED_S_  [h-hShift][v-vShift];
				_OBSERVED_NS_ [v-vShift][h-hShift] = _OBSERVED_NS_ [h-hShift][v-vShift];
				
				
			}
		}
	}
}