X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=binaries%2Fsrc%2Fmafft%2Fcore%2Fconstants.c;h=080ce137a36a9fe30d9fab77e967dd83530f2bc0;hb=711ff37b582bf25ba17d553008c4bf06cbdb8c84;hp=efc06dc39cf8bccf0943692a42c934545ea85867;hpb=7c6018cf87e67d7ac21f5230284ef2162a19320f;p=jabaws.git

diff --git a/binaries/src/mafft/core/constants.c b/binaries/src/mafft/core/constants.c
index efc06dc..080ce13 100644
--- a/binaries/src/mafft/core/constants.c
+++ b/binaries/src/mafft/core/constants.c
@@ -11,6 +11,7 @@
 
 #define NORMALIZE1 1
 
+
 static int shishagonyuu( double in )
 {
 	int out;
@@ -21,6 +22,71 @@ static int shishagonyuu( double in )
 	return( out );
 }
 
+static void nscore( int *amino_n, int **n_dis )
+{
+	int i;
+	for( i=0; i<26; i++ )
+	{
+//		reporterr( "i=%d (%c), n_dis[%d][%d] = %d\n", i, amino[i], i, amino_n['n'], n_dis[i][amino_n['n']] );
+		n_dis[i][amino_n['n']] = shishagonyuu( (double)0.25 * n_dis[i][i] );
+//		reporterr( "-> i=%d, n_dis[%d][%d] = %d\n", i, i, amino_n['n'], n_dis[i][amino_n['n']] );
+		n_dis[amino_n['n']][i] = n_dis[i][amino_n['n']];
+	}
+//	n_dis[amino_n['n']][amino_n['n']] = shishagonyuu( (double)0.25 * 0.25 * ( n_dis[0][0] + n_dis[1][1] + n_dis[2][2] + n_dis[3][3] ) );
+	n_dis[amino_n['n']][amino_n['n']] = shishagonyuu( (double)0.25 * ( n_dis[0][0] + n_dis[1][1] + n_dis[2][2] + n_dis[3][3] ) ); // 2017/Jan/2
+
+#if 0 // Ato de kakunin
+	for( i=0; i<26; i++ )
+	{
+		n_dis[i][amino_n['-']] = shishagonyuu( (double)0.25 * n_dis[i][i] );
+		n_dis[amino_n['-']][i] = n_dis[i][amino_n['-']];
+	}
+//	n_dis[amino_n['-']][amino_n['-']] = shishagonyuu( (double)0.25 * 0.25 * ( n_dis[0][0] + n_dis[1][1] + n_dis[2][2] + n_dis[3][3] ) ); // DAME!
+#endif
+}
+
+
+static void ambiguousscore( int *amino_n, int **n_dis )
+{
+	int i;
+	for( i=0; i<26; i++ )
+	{
+		n_dis[i][amino_n['r']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['g']][i] ) );
+		n_dis[i][amino_n['y']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['c']][i] + n_dis[amino_n['t']][i] ) );
+		n_dis[i][amino_n['k']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][i] + n_dis[amino_n['t']][i] ) );
+		n_dis[i][amino_n['m']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['c']][i] ) );
+		n_dis[i][amino_n['s']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][i] + n_dis[amino_n['c']][i] ) );
+		n_dis[i][amino_n['w']] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['t']][i] ) );
+		n_dis[i][amino_n['b']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['c']][i] + n_dis[amino_n['g']][i] + n_dis[amino_n['t']][i] ) );
+		n_dis[i][amino_n['d']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['g']][i] + n_dis[amino_n['t']][i] ) );
+		n_dis[i][amino_n['h']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['c']][i] + n_dis[amino_n['t']][i] ) );
+		n_dis[i][amino_n['v']] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][i] + n_dis[amino_n['c']][i] + n_dis[amino_n['g']][i] ) );
+
+		n_dis[amino_n['r']][i] = n_dis[i][amino_n['r']];
+		n_dis[amino_n['y']][i] = n_dis[i][amino_n['y']];
+		n_dis[amino_n['k']][i] = n_dis[i][amino_n['k']];
+		n_dis[amino_n['m']][i] = n_dis[i][amino_n['m']];
+		n_dis[amino_n['s']][i] = n_dis[i][amino_n['s']];
+		n_dis[amino_n['w']][i] = n_dis[i][amino_n['w']];
+		n_dis[amino_n['b']][i] = n_dis[i][amino_n['b']];
+		n_dis[amino_n['d']][i] = n_dis[i][amino_n['d']];
+		n_dis[amino_n['h']][i] = n_dis[i][amino_n['h']];
+		n_dis[amino_n['v']][i] = n_dis[i][amino_n['v']];
+	}
+
+	i = amino_n['r']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['g']][amino_n['g']] ) );
+	i = amino_n['y']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['t']][amino_n['t']] ) );
+	i = amino_n['k']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['t']][amino_n['t']] ) );
+	i = amino_n['m']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['c']][amino_n['c']] ) );
+	i = amino_n['s']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['c']][amino_n['c']] ) );
+	i = amino_n['w']; n_dis[i][i] = shishagonyuu( (double)1/2 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['t']][amino_n['t']] ) );
+	i = amino_n['b']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['t']][amino_n['t']] ) );
+	i = amino_n['d']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['g']][amino_n['g']] + n_dis[amino_n['t']][amino_n['t']] ) );
+	i = amino_n['h']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['t']][amino_n['t']] ) );
+	i = amino_n['v']; n_dis[i][i] = shishagonyuu( (double)1/3 * ( n_dis[amino_n['a']][amino_n['a']] + n_dis[amino_n['c']][amino_n['c']] + n_dis[amino_n['g']][amino_n['g']] ) );
+}
+
+
 static void calcfreq_nuc( int nseq, char **seq, double *datafreq )
 {
 	int i, j, l;
@@ -43,19 +109,20 @@ static void calcfreq_nuc( int nseq, char **seq, double *datafreq )
 			}
 		}
 	}
-	for( i=0; i<4; i++ )
-		if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
+	total = 0.0; for( i=0; i<4; i++ ) total += datafreq[i];
+	for( i=0; i<4; i++ ) datafreq[i] /= (double)total;
+	for( i=0; i<4; i++ ) if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
 
 
 	total = 0.0; for( i=0; i<4; i++ ) total += datafreq[i];
-//	fprintf( stderr, "total = %f\n", total );
+//	reporterr(       "total = %f\n", total );
 	for( i=0; i<4; i++ ) datafreq[i] /= (double)total;
 
 #if 0
-	fprintf( stderr, "\ndatafreq = " );
+	reporterr(       "\ndatafreq = " );
 	for( i=0; i<4; i++ )
-		fprintf( stderr, "%10.5f ", datafreq[i] );
-	fprintf( stderr, "\n" );
+		reporterr(       "%10.5f ", datafreq[i] );
+	reporterr(       "\n" );
 	exit( 1 );
 #endif
 }
@@ -65,7 +132,7 @@ static void calcfreq( int nseq, char **seq, double *datafreq )
 	int i, j, l;
 	int aan;
 	double total;
-	for( i=0; i<20; i++ )
+	for( i=0; i<nscoredalphabets; i++ )
 		datafreq[i] = 0.0;
 	total = 0.0;
 	for( i=0; i<nseq; i++ )
@@ -74,29 +141,100 @@ static void calcfreq( int nseq, char **seq, double *datafreq )
 		for( j=0; j<l; j++ )
 		{
 			aan = amino_n[(int)seq[i][j]];
-			if( aan >= 0 && aan < 20 )
+			if( aan >= 0 && aan < nscoredalphabets && seq[i][j] != '-' )
 			{
 				datafreq[aan] += 1.0;
 				total += 1.0;
 			}
 		}
 	}
-	for( i=0; i<20; i++ )
-		if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
+	total = 0.0; for( i=0; i<nscoredalphabets; i++ ) total += datafreq[i];
+	for( i=0; i<nscoredalphabets; i++ ) datafreq[i] /= (double)total;
+	for( i=0; i<nscoredalphabets; i++ ) if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
 
-	fprintf( stderr, "datafreq = \n" );
-	for( i=0; i<20; i++ )
-		fprintf( stderr, "%f\n", datafreq[i] );
+//	reporterr(       "datafreq = \n" );
+//	for( i=0; i<nscoredalphabets; i++ )
+//		reporterr(       "%f\n", datafreq[i] );
 
-	total = 0.0; for( i=0; i<20; i++ ) total += datafreq[i];
-	fprintf( stderr, "total = %f\n", total );
-	for( i=0; i<20; i++ ) datafreq[i] /= (double)total;
+	total = 0.0; for( i=0; i<nscoredalphabets; i++ ) total += datafreq[i];
+//	reporterr(       "total = %f\n", total );
+	for( i=0; i<nscoredalphabets; i++ ) datafreq[i] /= (double)total;
 }
 
+static void calcfreq_extended( int nseq, char **seq, double *datafreq )
+{
+	int i, j, l;
+	int aan;
+	double total;
+	for( i=0; i<nscoredalphabets; i++ )
+		datafreq[i] = 0.0;
+	total = 0.0;
+	for( i=0; i<nseq; i++ )
+	{
+		l = strlen( seq[i] );
+		for( j=0; j<l; j++ )
+		{
+			aan = amino_n[(unsigned char)seq[i][j]];
+			if( aan >= 0 && aan < nscoredalphabets && seq[i][j] != '-' )
+			{
+				datafreq[aan] += 1.0;
+				total += 1.0;
+			}
+		}
+	}
+	total = 0.0; for( i=0; i<nscoredalphabets; i++ ) total += datafreq[i];
+	for( i=0; i<nscoredalphabets; i++ ) datafreq[i] /= (double)total;
+//	for( i=0; i<nscoredalphabets; i++ ) if( datafreq[i] < 0.0001 ) datafreq[i] = 0.0001;
+
+#if 0
+	reporterr(       "datafreq = \n" );
+	for( i=0; i<nscoredalphabets; i++ )
+		reporterr(       "%d %c %f\n", i, amino[i], datafreq[i] );
+#endif
+
+	total = 0.0; for( i=0; i<nscoredalphabets; i++ ) total += datafreq[i];
+//	reporterr(       "total = %f\n", total );
+	for( i=0; i<nscoredalphabets; i++ ) datafreq[i] /= (double)total;
+}
+
+static void generatenuc1pam( double **pam1, int kimuraR, double *freq )
+{
+	int i, j;
+	double R[4][4], mut[4], total, tmp;
+
+	R[0][0] = 0.0;     R[0][1] = kimuraR; R[0][2] = 1.0;     R[0][3] = 1.0;
+	R[1][0] = kimuraR; R[1][1] = 0.0;     R[1][2] = 1.0;     R[1][3] = 1.0;
+	R[2][0] = 1.0;     R[2][1] = 1.0;     R[2][2] = 0.0;     R[2][3] = kimuraR;
+	R[3][0] = 1.0;     R[3][1] = 1.0;     R[3][2] = kimuraR; R[3][3] = 0.0;
+
+
+	total = 0.0;
+	for( i=0; i<4; i++ )
+	{
+		tmp = 0.0;
+		for( j=0; j<4; j++ ) tmp += R[i][j] * freq[j];
+		mut[i] = tmp;
+		total += tmp * freq[i];
+	}
+	for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
+	{
+		if( i != j ) pam1[i][j] = 0.01 / total * R[i][j] * freq[j];
+		else pam1[i][j] = 1.0 - 0.01 / total * mut[i];
+	}
+}
+
+
 void constants( int nseq, char **seq )
 {
 	int i, j, x;
 //	double tmp;
+	char shiftmodel[100];
+	int charsize;
+
+	if( nblosum < 0 ) dorp = 'p';
+
+	if( penalty_shift_factor >= 10 ) trywarp = 0;
+	else trywarp = 1;
 
 	if( dorp == 'd' )  /* DNA */
 	{
@@ -106,11 +244,18 @@ void constants( int nseq, char **seq )
 		double **pam1 = AllocateDoubleMtx( 4, 4 );
 		double *freq = AllocateDoubleVec( 4 );
 
+		nalphabets = 26;
+		nscoredalphabets = 10;
+		charsize = 0x80;
+
+		n_dis = AllocateIntMtx( nalphabets, nalphabets );
+		n_disLN = AllocateDoubleMtx( nalphabets, nalphabets );
 
 		scoremtx = -1;
 		if( RNAppenalty == NOTSPECIFIED ) RNAppenalty = DEFAULTRNAGOP_N;
 		if( RNAppenalty_ex == NOTSPECIFIED ) RNAppenalty_ex = DEFAULTRNAGEP_N;
 		if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_N;
+		if( ppenalty_dist == NOTSPECIFIED ) ppenalty_dist = ppenalty;
 		if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_N;
 		if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_N;
 		if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_N;
@@ -121,23 +266,48 @@ void constants( int nseq, char **seq )
 
 		RNApenalty = (int)( 3 * 600.0 / 1000.0 * RNAppenalty + 0.5 );
 		RNApenalty_ex = (int)( 3 * 600.0 / 1000.0 * RNAppenalty_ex + 0.5 );
-//		fprintf( stderr, "DEFAULTRNAGOP_N = %d\n", DEFAULTRNAGOP_N );
-//		fprintf( stderr, "RNAppenalty = %d\n", RNAppenalty );
-//		fprintf( stderr, "RNApenalty = %d\n", RNApenalty );
+//		reporterr(       "DEFAULTRNAGOP_N = %d\n", DEFAULTRNAGOP_N );
+//		reporterr(       "RNAppenalty = %d\n", RNAppenalty );
+//		reporterr(       "RNApenalty = %d\n", RNApenalty );
 
 
 		RNAthr = (int)( 3 * 600.0 / 1000.0 * RNApthr + 0.5 );
 		penalty = (int)( 3 * 600.0 / 1000.0 * ppenalty + 0.5);
+		penalty_dist = (int)( 3 * 600.0 / 1000.0 * ppenalty_dist + 0.5);
+		penalty_shift = (int)( penalty_shift_factor * penalty );
 		penalty_OP = (int)( 3 * 600.0 / 1000.0 * ppenalty_OP + 0.5);
 		penalty_ex = (int)( 3 * 600.0 / 1000.0 * ppenalty_ex + 0.5);
 		penalty_EX = (int)( 3 * 600.0 / 1000.0 * ppenalty_EX + 0.5);
-		offset = (int)( 3 * 600.0 / 1000.0 * poffset + 0.5);
-		offsetFFT = (int)( 3 * 600.0 / 1000.0 * (-0) + 0.5);
-		offsetLN = (int)( 3 * 600.0 / 1000.0 * 100 + 0.5);
+		offset = (int)( 1 * 600.0 / 1000.0 * poffset + 0.5);
+		offsetFFT = (int)( 1 * 600.0 / 1000.0 * (-0) + 0.5);
+		offsetLN = (int)( 1 * 600.0 / 1000.0 * 100 + 0.5);
 		penaltyLN = (int)( 3 * 600.0 / 1000.0 * -2000 + 0.5);
 		penalty_exLN = (int)( 3 * 600.0 / 1000.0 * -100 + 0.5);
-		sprintf( modelname, "%s%d (%d), %6.3f (%6.3f), %6.3f (%6.3f)", rnakozo?"RNA":"DNA", pamN, kimuraR,
-        -(double)ppenalty*0.001, -(double)ppenalty*0.003, -(double)poffset*0.001, -(double)poffset*0.003 );
+
+		if( trywarp ) sprintf( shiftmodel, "%4.2f (%4.2f)", -(double)penalty_shift/1800, -(double)penalty_shift/600 );
+		else sprintf( shiftmodel, "noshift" );
+
+		sprintf( modelname, "%s%d (%d), %4.2f (%4.2f), %4.2f (%4.2f), %s", rnakozo?"RNA":"DNA", pamN, kimuraR, -(double)ppenalty*0.001, -(double)ppenalty*0.003, -(double)poffset*0.001, -(double)poffset*0.003, shiftmodel );
+
+		for( i=0; i<26; i++ ) amino[i] = locaminon[i];
+		for( i=0; i<0x80; i++ ) amino_n[i] = -1;
+		for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
+		if( fmodel == 1 )
+		{
+			calcfreq_nuc( nseq, seq, freq );
+			reporterr(       "a, freq[0] = %f\n", freq[0] );
+			reporterr(       "g, freq[1] = %f\n", freq[1] );
+			reporterr(       "c, freq[2] = %f\n", freq[2] );
+			reporterr(       "t, freq[3] = %f\n", freq[3] );
+		}
+		else
+		{
+			freq[0] = 0.25;
+			freq[1] = 0.25;
+			freq[2] = 0.25;
+			freq[3] = 0.25;
+		}
+
 
 		if( kimuraR == 9999 ) 
 		{
@@ -150,7 +320,7 @@ void constants( int nseq, char **seq )
 			average /= 16.0;
 	
    	     if( disp )
-				fprintf( stderr, "average = %f\n", average );
+				reporterr(       "average = %f\n", average );
 	
 			for( i=0; i<4; i++ ) for( j=0; j<4; j++ ) 
 				pamx[i][j] -= average;
@@ -164,6 +334,7 @@ void constants( int nseq, char **seq )
 		}
 		else
 		{
+#if 0
 				double f = 0.99;
 				double s = (double)kimuraR / ( 2 + kimuraR ) * 0.01;
 				double v = (double)1       / ( 2 + kimuraR ) * 0.01;
@@ -171,32 +342,49 @@ void constants( int nseq, char **seq )
 				pam1[1][0] = s; pam1[1][1] = f; pam1[1][2] = v; pam1[1][3] = v;
 				pam1[2][0] = v; pam1[2][1] = v; pam1[2][2] = f; pam1[2][3] = s;
 				pam1[3][0] = v; pam1[3][1] = v; pam1[3][2] = s; pam1[3][3] = f;
+#else
+				generatenuc1pam( pam1, kimuraR, freq );
+#endif
 	
-				fprintf( stderr, "generating %dPAM scoring matrix for nucleotides ... ", pamN );
+				reporterr(       "generating a scoring matrix for nucleotide (dist=%d) ... ", pamN );
 	
 		       	if( disp )
    		    	{
-   		     		fprintf( stderr, " TPM \n" );
+   		     		reporterr(       " TPM \n" );
    		        	for( i=0; i<4; i++ )
    			       	{
    		            	for( j=0; j<4; j++ )
-   		                	fprintf( stderr, "%+#6.10f", pam1[i][j] );
-   		            	fprintf( stderr, "\n" );
+   		                	reporterr(       "%+#6.10f", pam1[i][j] );
+   		            	reporterr(       "\n" );
    		        	}
-   		        	fprintf( stderr, "\n" );
+   		        	reporterr(       "\n" );
    		     	}
 	
 	
 				MtxuntDouble( pamx, 4 );
 				for( x=0; x < pamN; x++ ) MtxmltDouble( pamx, pam1, 4 );
+
+		       	if( disp )
+   		    	{
+   		     		reporterr(       " TPM \n" );
+   		        	for( i=0; i<4; i++ )
+   			       	{
+   		            	for( j=0; j<4; j++ )
+   		                	reporterr(       "%+#6.10f", pamx[i][j] );
+   		            	reporterr(       "\n" );
+   		        	}
+   		        	reporterr(       "\n" );
+   		     	}
+
 				for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
-					pamx[i][j] /= 1.0 / 4.0;
+					pamx[i][j] /= freq[j];
+//					pamx[i][j] /= 0.25;
 	
 				for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
 				{
 					if( pamx[i][j] == 0.0 ) 
 					{
-						fprintf( stderr, "WARNING: pamx[i][j] = 0.0 ?\n" );
+						reporterr(       "WARNING: pamx[i][j] = 0.0 ?\n" );
 						pamx[i][j] = 0.00001; /* by J. Thompson */
 					}
 					pamx[i][j] = log10( pamx[i][j] ) * 1000.0;
@@ -204,35 +392,18 @@ void constants( int nseq, char **seq )
 	
    	    		if( disp )
    	    		{
-   		     		fprintf( stderr, " after log\n" );
+   		     		reporterr(       " after log\n" );
    	        		for( i=0; i<4; i++ )
    		       		{
    	        	    	for( j=0; j<4; j++ )
-   	        	        	fprintf( stderr, "%+#6.10f", pamx[i][j] );
-   	        	    	fprintf( stderr, "\n" );
+   	        	        	reporterr(       "%+10.6f ", pamx[i][j] );
+   	        	    	reporterr(       "\n" );
    	        		}
-   	        		fprintf( stderr, "\n" );
+   	        		reporterr(       "\n" );
    		     	}
 
 
 // ?????
-			for( i=0; i<26; i++ ) amino[i] = locaminon[i];
-			for( i=0; i<0x80; i++ ) amino_n[i] = -1;
-			for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
-			if( fmodel == 1 )
-				calcfreq_nuc( nseq, seq, freq );
-			else
-			{
-				freq[0] = 0.25;
-				freq[1] = 0.25;
-				freq[2] = 0.25;
-				freq[3] = 0.25;
-			}
-//			fprintf( stderr, "a, freq[0] = %f\n", freq[0] );
-//			fprintf( stderr, "g, freq[1] = %f\n", freq[1] );
-//			fprintf( stderr, "c, freq[2] = %f\n", freq[2] );
-//			fprintf( stderr, "t, freq[3] = %f\n", freq[3] );
-
 			
 			average = 0.0;
 			for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
@@ -249,23 +420,23 @@ void constants( int nseq, char **seq )
 
 
 			for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
-				pamx[i][j] -= offset;        /* extending gap cost */
+				pamx[i][j] -= offset;
 
 			for( i=0; i<4; i++ ) for( j=0; j<4; j++ )
 				pamx[i][j] = shishagonyuu( pamx[i][j] );
 
        		if( disp )
        		{
-        		fprintf( stderr, " after shishagonyuu\n" );
+        		reporterr(       " after shishagonyuu\n" );
            		for( i=0; i<4; i++ )
    	       		{
            	    	for( j=0; j<4; j++ )
-           	        	fprintf( stderr, "%+#6.10f", pamx[i][j] );
-           	    	fprintf( stderr, "\n" );
+           	        	reporterr(       "%+#6.10f", pamx[i][j] );
+           	    	reporterr(       "\n" );
            		}
-           		fprintf( stderr, "\n" );
+           		reporterr(       "\n" );
         	}
-			fprintf( stderr, "done\n" );
+			reporterr(       "done\n" );
 		}
 	
 		for( i=0; i<5; i++ ) 
@@ -281,42 +452,49 @@ void constants( int nseq, char **seq )
 	
        	if( disp )
        	{
-       		fprintf( stderr, " before dis\n" );
+       		reporterr(       " before dis\n" );
           	for( i=0; i<4; i++ )
    	       	{
            	   	for( j=0; j<4; j++ )
-           	       	fprintf( stderr, "%+#6.10f", pamx[i][j] );
-           	   	fprintf( stderr, "\n" );
+           	       	reporterr(       "%+#6.10f", pamx[i][j] );
+           	   	reporterr(       "\n" );
            	}
-           	fprintf( stderr, "\n" );
+           	reporterr(       "\n" );
         }
 
        	if( disp )
        	{
-        	fprintf( stderr, " score matrix  \n" );
+        	reporterr(       " score matrix  \n" );
            	for( i=0; i<4; i++ )
    	       	{
                	for( j=0; j<4; j++ )
-                   	fprintf( stderr, "%+#6.10f", pamx[i][j] );
-               	fprintf( stderr, "\n" );
+                   	reporterr(       "%+#6.10f", pamx[i][j] );
+               	reporterr(       "\n" );
            	}
-           	fprintf( stderr, "\n" );
+           	reporterr(       "\n" );
+					exit( 1 );
         }
 
 		for( i=0; i<26; i++ ) amino[i] = locaminon[i];
 		for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpn[i];
 		for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = 0;
 		for( i=0; i<10; i++ ) for( j=0; j<10; j++ ) n_dis[i][j] = shishagonyuu( pamx[i][j] );
+
+		ambiguousscore( amino_n, n_dis );
+		if( nwildcard ) nscore( amino_n, n_dis );
+
         if( disp )
         {
-            fprintf( stderr, " score matrix  \n" );
+            reporterr(       " score matrix  \n" );
             for( i=0; i<26; i++ )
             {
                 for( j=0; j<26; j++ )
-                    fprintf( stderr, "%+6d", n_dis[i][j] );
-                fprintf( stderr, "\n" );
+                    reporterr(       "%+6d", n_dis[i][j] );
+                reporterr(       "\n" );
             }
-            fprintf( stderr, "\n" );
+            reporterr(       "\n" );
+			reporterr(       "penalty = %d, penalty_ex = %d\n", penalty, penalty_ex );
+//exit( 1 );
         }
 
 // RIBOSUM
@@ -368,24 +546,24 @@ void constants( int nseq, char **seq )
 
   		if( disp )
    		{
-     		fprintf( stderr, "ribosum after shishagonyuu\n" );
+     		reporterr(       "ribosum after shishagonyuu\n" );
        		for( i=0; i<4; i++ )
        		{
        	    	for( j=0; j<4; j++ )
-       	        	fprintf( stderr, "%+#6.10f", ribosum4[i][j] );
-       	    	fprintf( stderr, "\n" );
+       	        	reporterr(       "%+#6.10f", ribosum4[i][j] );
+       	    	reporterr(       "\n" );
        		}
-       		fprintf( stderr, "\n" );
-     		fprintf( stderr, "ribosum16 after shishagonyuu\n" );
+       		reporterr(       "\n" );
+     		reporterr(       "ribosum16 after shishagonyuu\n" );
        		for( i=0; i<16; i++ )
        		{
        	    	for( j=0; j<16; j++ )
-       	        	fprintf( stderr, "%+#7.0f", ribosum16[i][j] );
-       	    	fprintf( stderr, "\n" );
+       	        	reporterr(       "%+#7.0f", ribosum16[i][j] );
+       	    	reporterr(       "\n" );
        		}
-       		fprintf( stderr, "\n" );
+       		reporterr(       "\n" );
       	}
-		fprintf( stderr, "done\n" );
+//		reporterr(       "done\n" );
 
 #if 1
 		for( i=0; i<37; i++ ) for( j=0; j<37; j++ ) ribosumdis[i][j] = 0.0; //iru
@@ -403,16 +581,16 @@ void constants( int nseq, char **seq )
 
   		if( disp )
    		{
-     		fprintf( stderr, "ribosumdis\n" );
+     		reporterr(       "ribosumdis\n" );
        		for( i=0; i<37; i++ )
        		{
        	    	for( j=0; j<37; j++ )
-       	        	fprintf( stderr, "%+5d", ribosumdis[i][j] );
-       	    	fprintf( stderr, "\n" );
+       	        	reporterr(       "%+5d", ribosumdis[i][j] );
+       	    	reporterr(       "\n" );
        		}
-       		fprintf( stderr, "\n" );
+       		reporterr(       "\n" );
       	}
-		fprintf( stderr, "done\n" );
+//		reporterr(       "done\n" );
 #endif
 
 		FreeDoubleMtx( pam1 );
@@ -420,7 +598,7 @@ void constants( int nseq, char **seq )
 		free( freq );
 
 	}
-	else if( dorp == 'p' && scoremtx == 1 )  /* Blosum */
+	else if( dorp == 'p' && scoremtx == 1 && nblosum == -2 )  /* extended */
 	{
 		double *freq;
 		double *freq1;
@@ -429,11 +607,230 @@ void constants( int nseq, char **seq )
 //		double tmp;
 		double **n_distmp;
 
+		nalphabets = 0x100;
+		nscoredalphabets = 0x100;
+		charsize = 0x100;
+
+		reporterr(       "nalphabets = %d\n", nalphabets );
+
+		n_dis = AllocateIntMtx( nalphabets, nalphabets );
+		n_disLN = AllocateDoubleMtx( nalphabets, nalphabets );
+		n_distmp = AllocateDoubleMtx( nalphabets, nalphabets );
+		datafreq = AllocateDoubleVec( nalphabets );
+		freq = AllocateDoubleVec( nalphabets );
+
+		if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_B;
+		if( ppenalty_dist == NOTSPECIFIED ) ppenalty_dist = ppenalty;
+		if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_B;
+		if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_B;
+		if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_B;
+		if( poffset == NOTSPECIFIED ) poffset = DEFAULTOFS_B;
+		if( pamN == NOTSPECIFIED ) pamN = 0;
+		if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
+		penalty = (int)( 600.0 / 1000.0 * ppenalty + 0.5 );
+		penalty_dist = (int)( 600.0 / 1000.0 * ppenalty_dist + 0.5 );
+		penalty_shift = (int)( penalty_shift_factor * penalty );
+		penalty_OP = (int)( 600.0 / 1000.0 * ppenalty_OP + 0.5 );
+		penalty_ex = (int)( 600.0 / 1000.0 * ppenalty_ex + 0.5 );
+		penalty_EX = (int)( 600.0 / 1000.0 * ppenalty_EX + 0.5 );
+		offset = (int)( 600.0 / 1000.0 * poffset + 0.5 );
+		offsetFFT = (int)( 600.0 / 1000.0 * (-0) + 0.5);
+		offsetLN = (int)( 600.0 / 1000.0 * 100 + 0.5);
+		penaltyLN = (int)( 600.0 / 1000.0 * -2000 + 0.5);
+		penalty_exLN = (int)( 600.0 / 1000.0 * -100 + 0.5);
+
+		extendedmtx( n_distmp, freq, amino, amino_grp );
+
+		if( trywarp ) sprintf( shiftmodel, "%4.2f", -(double)penalty_shift/600 );
+		else sprintf( shiftmodel, "noshift" );
+
+		sprintf( modelname, "Extended, %4.2f, %+4.2f, %+4.2f, %s", -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000, shiftmodel );
+#if 0
+		for( i=0; i<26; i++ ) amino[i] = locaminod[i];
+		for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpd[i];
+		for( i=0; i<0x80; i++ ) amino_n[i] = 0;
+		for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
+#endif
+		for( i=0; i<0x100; i++ )amino_n[i] = -1;
+		for( i=0; i<nalphabets; i++) 
+		{
+			amino_n[(unsigned char)amino[i]] = i;
+//			reporterr(       "i=%d, amino = %c, amino_n = %d\n", i, amino[i], amino_n[amino[i]] );
+		}
+		if( fmodel == 1 )
+		{
+			calcfreq_extended( nseq, seq, datafreq );
+			freq1 = datafreq;
+		}
+		else
+			freq1 = freq;
+
+#if TEST
+		reporterr(       "raw scoreing matrix : \n" );
+		for( i=0; i<nalphabets; i++ )
+		{
+			for( j=0; j<nalphabets; j++ ) 
+			{
+				fprintf( stdout, "%6.2f", n_distmp[i][j] );
+			}
+			fprintf( stdout, "\n" );
+		}
+#endif
+		if( fmodel == -1 )
+			average = 0.0;
+		else
+		{
+			for( i=0; i<nalphabets; i++ )
+#if TEST 
+				fprintf( stdout, "freq[%c] = %f, datafreq[%c] = %f, freq1[] = %f\n", amino[i], freq[i], amino[i], datafreq[i], freq1[i] );
+#endif
+			average = 0.0;
+			for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ )
+				average += n_distmp[i][j] * freq1[i] * freq1[j];
+		}
+#if TEST
+		fprintf( stdout, "####### average2  = %f\n", average );
+#endif
+
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) 
+			n_distmp[i][j] -= average;
+#if TEST
+		fprintf( stdout, "average2 = %f\n", average );
+		fprintf( stdout, "after average substruction : \n" );
+		for( i=0; i<nalphabets; i++ )
+		{
+			for( j=0; j<nalphabets; j++ ) 
+			{
+				fprintf( stdout, "%6.2f", n_distmp[i][j] );
+			}
+			fprintf( stdout, "\n" );
+		}
+#endif
+		
+		average = 0.0;
+		for( i=0; i<nalphabets; i++ ) 
+			average += n_distmp[i][i] * freq1[i];
+#if TEST
+		fprintf( stdout, "####### average1  = %f\n", average );
+#endif
+
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) 
+			n_distmp[i][j] *= 600.0 / average;
+#if TEST
+        fprintf( stdout, "after average division : \n" );
+        for( i=0; i<nalphabets; i++ )
+        {
+            for( j=0; j<=i; j++ )
+            {
+                fprintf( stdout, "%7.1f", n_distmp[i][j] );
+            }
+            fprintf( stdout, "\n" );
+        }
+#endif
+
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) 
+			n_distmp[i][j] -= offset;  
+#if TEST
+		fprintf( stdout, "after offset substruction (offset = %d): \n", offset );
+		for( i=0; i<nalphabets; i++ )
+		{
+			for( j=0; j<=i; j++ ) 
+			{
+				fprintf( stdout, "%7.1f", n_distmp[i][j] );
+			}
+			fprintf( stdout, "\n" );
+		}
+#endif
+#if 0
+/* Ãí°Õ ¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª¡ª */
+			penalty -= offset;
+#endif
+
+
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) 
+			n_distmp[i][j] = shishagonyuu( n_distmp[i][j] );
+
+        if( disp )
+        {
+			fprintf( stdout, "freq = \n" );
+			for( i=0; i<nalphabets; i++ ) fprintf( stdout, "%c %f\n", amino[i], freq1[i] );
+            fprintf( stdout, " scoring matrix  \n" );
+            for( i=0; i<nalphabets; i++ )
+            {
+				fprintf( stdout, "%c    ", amino[i] );
+                for( j=0; j<nalphabets; j++ )
+                    fprintf( stdout, "%5.0f", n_distmp[i][j] );
+                fprintf( stdout, "\n" );
+            }
+			fprintf( stdout, "     " );
+            for( i=0; i<nalphabets; i++ )
+				fprintf( stdout, "    %c", amino[i] );
+
+			average = 0.0;
+        	for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ )
+				average += n_distmp[i][j] * freq1[i] * freq1[j];
+			fprintf( stdout, "average = %f\n", average );
+
+			average = 0.0;
+        	for( i=0; i<nalphabets; i++ )
+				average += n_distmp[i][i] * freq1[i];
+			fprintf( stdout, "itch average = %f\n", average );
+			reporterr(       "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
+
+			
+  			exit( 1 );
+        }
+
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) n_dis[i][j] = 0;
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) n_dis[i][j] = (int)n_distmp[i][j];
+		for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) n_dis[i][amino_n['-']] = n_dis[amino_n['-']][i] = 0.0;
+
+		FreeDoubleMtx( n_distmp );
+		FreeDoubleVec( datafreq );
+		FreeDoubleVec( freq );
+
+//		reporterr(       "done.\n" );
+
+	}
+	else if( dorp == 'p' && scoremtx == 1 )  /* Blosum, user-defined */
+	{
+		double *freq;
+		double *freq1;
+		double *datafreq;
+		double average;
+		double iaverage;
+//		double tmp;
+		double **n_distmp;
+		int makeaverage0;
+
+
+		if( nblosum == 0 )
+		{
+			reporterr( "nblosum=%d??\n", nblosum );
+			exit( 1 );
+		}
+//		if( nblosum < 0 )
+//		{
+//			nblosum *= -1;
+//			makeaverage0 = 0;
+//		}
+		else
+		{
+			makeaverage0 = 1;
+		}
+
+		nalphabets = 26;
+		nscoredalphabets = 20;
+		charsize = 0x80;
+
+		n_dis = AllocateIntMtx( nalphabets, nalphabets );
+		n_disLN = AllocateDoubleMtx( nalphabets, nalphabets );
 		n_distmp = AllocateDoubleMtx( 20, 20 );
 		datafreq = AllocateDoubleVec( 20 );
 		freq = AllocateDoubleVec( 20 );
 
 		if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_B;
+		if( ppenalty_dist == NOTSPECIFIED ) ppenalty_dist = ppenalty;
 		if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_B;
 		if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_B;
 		if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_B;
@@ -441,6 +838,8 @@ void constants( int nseq, char **seq )
 		if( pamN == NOTSPECIFIED ) pamN = 0;
 		if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
 		penalty = (int)( 600.0 / 1000.0 * ppenalty + 0.5 );
+		penalty_dist = (int)( 600.0 / 1000.0 * ppenalty_dist + 0.5 );
+		penalty_shift = (int)( penalty_shift_factor * penalty );
 		penalty_OP = (int)( 600.0 / 1000.0 * ppenalty_OP + 0.5 );
 		penalty_ex = (int)( 600.0 / 1000.0 * ppenalty_ex + 0.5 );
 		penalty_EX = (int)( 600.0 / 1000.0 * ppenalty_EX + 0.5 );
@@ -451,10 +850,14 @@ void constants( int nseq, char **seq )
 		penalty_exLN = (int)( 600.0 / 1000.0 * -100 + 0.5);
 
 		BLOSUMmtx( nblosum, n_distmp, freq, amino, amino_grp );
+
+		if( trywarp ) sprintf( shiftmodel, "%4.2f", -(double)penalty_shift/600 );
+		else sprintf( shiftmodel, "noshift" );
+
 		if( nblosum == -1 )
-			sprintf( modelname, "User-defined, %6.3f, %+6.3f, %+6.3f", -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000 );
+			sprintf( modelname, "User-defined, %4.2f, %+4.2f, %+4.2f, %s", -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000, shiftmodel );
 		else
-			sprintf( modelname, "BLOSUM%d, %6.3f, %+6.3f, %+6.3f", nblosum, -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000 );
+			sprintf( modelname, "BLOSUM%d, %4.2f, %+4.2f, %+4.2f, %s", nblosum, -(double)ppenalty/1000, -(double)poffset/1000, -(double)ppenalty_ex/1000, shiftmodel );
 #if 0
 		for( i=0; i<26; i++ ) amino[i] = locaminod[i];
 		for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpd[i];
@@ -471,7 +874,7 @@ void constants( int nseq, char **seq )
 		else
 			freq1 = freq;
 #if TEST
-		fprintf( stderr, "raw scoreing matrix : \n" );
+		reporterr(       "raw scoreing matrix : \n" );
 		for( i=0; i<20; i++ )
 		{
 			for( j=0; j<20; j++ ) 
@@ -497,8 +900,11 @@ void constants( int nseq, char **seq )
 		fprintf( stdout, "####### average2  = %f\n", average );
 #endif
 
-		for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) 
-			n_distmp[i][j] -= average;
+		if( makeaverage0 )
+		{
+			for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) 
+				n_distmp[i][j] -= average;
+		}
 #if TEST
 		fprintf( stdout, "average2 = %f\n", average );
 		fprintf( stdout, "after average substruction : \n" );
@@ -572,13 +978,13 @@ void constants( int nseq, char **seq )
 			average = 0.0;
         	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
 				average += n_distmp[i][j] * freq1[i] * freq1[j];
-			fprintf( stdout, "average = %f\n", average );
+			fprintf( stdout, "\naverage = %f\n", average );
 
-			average = 0.0;
+			iaverage = 0.0;
         	for( i=0; i<20; i++ )
-				average += n_distmp[i][i] * freq1[i];
-			fprintf( stdout, "itch average = %f\n", average );
-			fprintf( stderr, "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
+				iaverage += n_distmp[i][i] * freq1[i];
+			fprintf( stdout, "itch average = %f, E=%f\n", iaverage, average/iaverage );
+			reporterr(       "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
 
 			
   			exit( 1 );
@@ -591,38 +997,13 @@ void constants( int nseq, char **seq )
 		FreeDoubleVec( datafreq );
 		FreeDoubleVec( freq );
 
-		fprintf( stderr, "done.\n" );
+//		reporterr(       "done.\n" );
 
 	}
 	else if( dorp == 'p' && scoremtx == 2 ) /* Miyata-Yasunaga */
 	{
-		fprintf( stderr, "Not supported\n" );
+		reporterr(       "Not supported\n" );
 		exit( 1 );
-		for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = locn_dism[i][j];
-		for( i=0; i<26; i++ ) if( i != 24 ) n_dis[i][24] = n_dis[24][i] = exgpm;
-		n_dis[24][24] = 0;
-		if( ppenalty == NOTSPECIFIED ) ppenalty = locpenaltym;
-		if( poffset == NOTSPECIFIED ) poffset = -20;
-		if( pamN == NOTSPECIFIED ) pamN = 0;
-		if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
-
-		penalty = ppenalty;
-		offset = poffset;
-
-		sprintf( modelname, "Miyata-Yasunaga, %6.3f, %6.3f", -(double)ppenalty/1000, -(double)poffset/1000 );
-		for( i=0; i<26; i++ ) amino[i] = locaminom[i];
-		for( i=0; i<26; i++ ) amino_grp[(int)amino[i]] = locgrpm[i];
-#if DEBUG
-		fprintf( stdout, "scoreing matrix : \n" );
-		for( i=0; i<26; i++ )
-		{
-			for( j=0; j<26; j++ ) 
-			{
-				fprintf( stdout, "%#5d", n_dis[i][j] );
-			}
-			fprintf( stdout, "\n" );
-		}
-#endif
 	}
 	else         /* JTT */
 	{
@@ -634,9 +1015,17 @@ void constants( int nseq, char **seq )
 		double *mutab;
 		double *datafreq;
 		double average;
+		double iaverage;
 		double tmp;
 		double delta;
+		int makeaverage0;
+
+		nalphabets = 26;
+		nscoredalphabets = 20;
+		charsize = 0x80;
 
+		n_dis = AllocateIntMtx( nalphabets, nalphabets );
+		n_disLN = AllocateDoubleMtx( nalphabets, nalphabets );
 		rsr = AllocateDoubleMtx( 20, 20 );
 		pam1 = AllocateDoubleMtx( 20, 20 );
 		pamx = AllocateDoubleMtx( 20, 20 );
@@ -645,6 +1034,7 @@ void constants( int nseq, char **seq )
 		datafreq = AllocateDoubleVec( 20 );
 
 		if( ppenalty == NOTSPECIFIED ) ppenalty = DEFAULTGOP_J;
+		if( ppenalty_dist == NOTSPECIFIED ) ppenalty_dist = ppenalty;
 		if( ppenalty_OP == NOTSPECIFIED ) ppenalty_OP = DEFAULTGOP_J;
 		if( ppenalty_ex == NOTSPECIFIED ) ppenalty_ex = DEFAULTGEP_J;
 		if( ppenalty_EX == NOTSPECIFIED ) ppenalty_EX = DEFAULTGEP_J;
@@ -652,7 +1042,24 @@ void constants( int nseq, char **seq )
 		if( pamN == NOTSPECIFIED )    pamN    = DEFAULTPAMN;
 		if( kimuraR == NOTSPECIFIED ) kimuraR = 1;
 
+		if( pamN == 0 )
+		{
+			reporterr( "pamN=%d??\n", pamN );
+			exit( 1 );
+		}
+		if( pamN < 0 )
+		{
+			pamN *= -1;
+			makeaverage0 = 0;
+		}
+		else
+		{
+			makeaverage0 = 1;
+		}
+
 		penalty = (int)( 600.0 / 1000.0 * ppenalty + 0.5 );
+		penalty_dist = (int)( 600.0 / 1000.0 * ppenalty_dist + 0.5 );
+		penalty_shift = (int)( penalty_shift_factor * penalty );
 		penalty_OP = (int)( 600.0 / 1000.0 * ppenalty_OP + 0.5 );
 		penalty_ex = (int)( 600.0 / 1000.0 * ppenalty_ex + 0.5 );
 		penalty_EX = (int)( 600.0 / 1000.0 * ppenalty_EX + 0.5 );
@@ -662,10 +1069,24 @@ void constants( int nseq, char **seq )
 		penaltyLN = (int)( 600.0 / 1000.0 * -2000 + 0.5);
 		penalty_exLN = (int)( 600.0 / 1000.0 * -100 + 0.5);
 
-		sprintf( modelname, "%s %dPAM, %6.3f, %6.3f", (TMorJTT==TM)?"Transmembrane":"JTT", pamN, -(double)ppenalty/1000, -(double)poffset/1000 );
+		if( trywarp ) sprintf( shiftmodel, "%4.2f", -(double)penalty_shift/600 );
+		else sprintf( shiftmodel, "noshift" );
+
+		sprintf( modelname, "%s %dPAM, %4.2f, %4.2f, %s", (TMorJTT==TM)?"Transmembrane":"JTT", pamN, -(double)ppenalty/1000, -(double)poffset/1000, shiftmodel );
 
 		JTTmtx( rsr, freq, amino, amino_grp, (int)(TMorJTT==TM) );
 
+		for( i=0; i<0x80; i++ ) amino_n[i] = -1;
+		for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
+		if( fmodel == 1 )
+		{
+			calcfreq( nseq, seq, datafreq );
+			freq1 = datafreq;
+		}
+		else
+			freq1 = freq;
+
+
 #if TEST
 		fprintf( stdout, "rsr = \n" );
 		for( i=0; i<20; i++ )
@@ -678,26 +1099,16 @@ void constants( int nseq, char **seq )
 		}
 #endif
 
-		for( i=0; i<0x80; i++ ) amino_n[i] = -1;
-		for( i=0; i<26; i++ ) amino_n[(int)amino[i]] = i;
 
-		if( fmodel == 1 )
-		{
-			calcfreq( nseq, seq, datafreq );
-			freq1 = datafreq;
-		}
-		else
-			freq1 = freq;
-
-		fprintf( stderr, "generating %dPAM %s scoring matrix for amino acids ... ", pamN, (TMorJTT==TM)?"Transmembrane":"JTT" );
+		reporterr(       "generating %dPAM %s scoring matrix for amino acids ... ", pamN, (TMorJTT==TM)?"Transmembrane":"JTT" );
 
 		tmp = 0.0;
 		for( i=0; i<20; i++ )
 		{
 			mutab[i] = 0.0;
 			for( j=0; j<20; j++ )
-				mutab[i] += rsr[i][j] * freq[j];
-			tmp += mutab[i] * freq[i];
+				mutab[i] += rsr[i][j] * freq1[j];
+			tmp += mutab[i] * freq1[i];
 		}
 #if TEST
 		fprintf( stdout, "mutability = \n" );
@@ -712,7 +1123,7 @@ void constants( int nseq, char **seq )
 			for( j=0; j<20; j++ )
 			{
 				if( i != j )
-					pam1[i][j] = delta * rsr[i][j] * freq[i];
+					pam1[i][j] = delta * rsr[i][j] * freq1[j];
 				else
 					pam1[i][j] = 1.0 - delta * mutab[i];
 			}
@@ -735,13 +1146,13 @@ void constants( int nseq, char **seq )
 		for( x=0; x < pamN; x++ ) MtxmltDouble( pamx, pam1, 20 );
 
 		for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) 
-			pamx[i][j] /= freq[j];
+			pamx[i][j] /= freq1[j];
 
         for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
 		{
 			if( pamx[i][j] == 0.0 ) 
 			{
-				fprintf( stderr, "WARNING: pamx[%d][%d] = 0.0?\n", i, j );
+				reporterr(       "WARNING: pamx[%d][%d] = 0.0?\n", i, j );
 				pamx[i][j] = 0.00001; /* by J. Thompson */
 			}
             pamx[i][j] = log10( pamx[i][j] ) * 1000.0;
@@ -800,8 +1211,11 @@ void constants( int nseq, char **seq )
 		fprintf( stdout, "####### average2  = %f\n", average );
 #endif
 
-		for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) 
-			pamx[i][j] -= average;
+		if( makeaverage0 )
+		{
+			for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) 
+				pamx[i][j] -= average;
+		}
 #if TEST
 		fprintf( stdout, "average2 = %f\n", average );
 		fprintf( stdout, "after average substruction : \n" );
@@ -888,13 +1302,13 @@ void constants( int nseq, char **seq )
 			average = 0.0;
         	for( i=0; i<20; i++ ) for( j=0; j<20; j++ )
 				average += pamx[i][j] * freq1[i] * freq1[j];
-			fprintf( stdout, "average = %f\n", average );
+			fprintf( stdout, "\naverage = %f\n", average );
 
-			average = 0.0;
+			iaverage = 0.0;
         	for( i=0; i<20; i++ )
-				average += pamx[i][i] * freq1[i];
-			fprintf( stdout, "itch average = %f\n", average );
-			fprintf( stderr, "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
+				iaverage += pamx[i][i] * freq1[i];
+			fprintf( stdout, "itch average = %f, E=%f\n", average, average/iaverage );
+			reporterr(       "parameters: %d, %d, %d\n", penalty, penalty_ex, offset );
 
 			
   			exit( 1 );
@@ -903,7 +1317,7 @@ void constants( int nseq, char **seq )
 		for( i=0; i<26; i++ ) for( j=0; j<26; j++ ) n_dis[i][j] = 0;
 		for( i=0; i<20; i++ ) for( j=0; j<20; j++ ) n_dis[i][j] = (int)pamx[i][j];
 
-		fprintf( stderr, "done.\n" );
+		reporterr(       "done.\n" );
 		FreeDoubleMtx( rsr );
 		FreeDoubleMtx( pam1 );
 		FreeDoubleMtx( pamx );
@@ -911,83 +1325,97 @@ void constants( int nseq, char **seq )
 		FreeDoubleVec( mutab );
 		FreeDoubleVec( datafreq );
 	}
-	fprintf( stderr, "scoremtx = %d\n", scoremtx );
 
 #if DEBUG
-	fprintf( stderr, "scoremtx = %d\n", scoremtx );
-	fprintf( stderr, "amino[] = %s\n", amino );
+	reporterr(       "scoremtx = %d\n", scoremtx );
+	reporterr(       "amino[] = %s\n", amino );
 #endif
 
-	for( i=0; i<0x80; i++ )amino_n[i] = -1;
-	for( i=0; i<26; i++) amino_n[(int)amino[i]] = i;
-    for( i=0; i<0x80; i++ ) for( j=0; j<0x80; j++ ) amino_dis[i][j] = 0;
-    for( i=0; i<0x80; i++ ) for( j=0; j<0x80; j++ ) amino_disLN[i][j] = 0;
-    for( i=0; i<0x80; i++ ) for( j=0; j<0x80; j++ ) amino_dis_consweight_multi[i][j] = 0.0;
-    for( i=0; i<26; i++) for( j=0; j<26; j++ )
+	amino_dis = AllocateIntMtx( charsize, charsize );
+	amino_dis_consweight_multi = AllocateDoubleMtx( charsize, charsize );
+
+//	reporterr( "charsize=%d\n", charsize );
+
+	for( i=0; i<charsize; i++ )amino_n[i] = -1;
+	for( i=0; i<nalphabets; i++) amino_n[(int)amino[i]] = i;
+    for( i=0; i<charsize; i++ ) for( j=0; j<charsize; j++ ) amino_dis[i][j] = 0;
+    for( i=0; i<nalphabets; i++ ) for( j=0; j<nalphabets; j++ ) n_disLN[i][j] = 0;
+    for( i=0; i<charsize; i++ ) for( j=0; j<charsize; j++ ) amino_dis_consweight_multi[i][j] = 0.0;
+
+	n_dis_consweight_multi = AllocateDoubleMtx( nalphabets, nalphabets );
+	n_disFFT = AllocateIntMtx( nalphabets, nalphabets );
+    for( i=0; i<nalphabets; i++) for( j=0; j<nalphabets; j++ )
 	{
         amino_dis[(int)amino[i]][(int)amino[j]] = n_dis[i][j];
-		n_dis_consweight_multi[i][j] = (float)n_dis[i][j] * consweight_multi;
+		n_dis_consweight_multi[i][j] = (double)n_dis[i][j] * consweight_multi;
 		amino_dis_consweight_multi[(int)amino[i]][(int)amino[j]] = (double)n_dis[i][j] * consweight_multi;
 	}
 
 	if( dorp == 'd' )  /* DNA */
 	{
+#if 0 // ???
 	    for( i=0; i<5; i++) for( j=0; j<5; j++ )
-        	amino_disLN[(int)amino[i]][(int)amino[j]] = n_dis[i][j] + offset - offsetLN;
+        	n_disLN[i][j] = (double)n_dis[i][j] + offset - offsetLN;
 	    for( i=5; i<10; i++) for( j=5; j<10; j++ )
-        	amino_disLN[(int)amino[i]][(int)amino[j]] = n_dis[i][j] + offset - offsetLN;
+        	n_disLN[i][j] = (double)n_dis[i][j] + offset - offsetLN;
 	    for( i=0; i<5; i++) for( j=0; j<5; j++ )
         	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
 	    for( i=5; i<10; i++) for( j=5; j<10; j++ )
         	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
+#else
+	    for( i=0; i<10; i++) for( j=0; j<10; j++ )
+        	n_disLN[i][j] = (double)n_dis[i][j] + offset - offsetLN;
+	    for( i=0; i<10; i++) for( j=0; j<10; j++ )
+        	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
+#endif
 	}
 	else // protein
 	{
 	    for( i=0; i<20; i++) for( j=0; j<20; j++ )
-        	amino_disLN[(int)amino[i]][(int)amino[j]] = n_dis[i][j] + offset - offsetLN;
+        	n_disLN[i][j] = (double)n_dis[i][j] + offset - offsetLN;
 	    for( i=0; i<20; i++) for( j=0; j<20; j++ )
         	n_disFFT[i][j] = n_dis[i][j] + offset - offsetFFT;
 	}
 
 #if 0
-		fprintf( stderr, "amino_dis (offset = %d): \n", offset );
+		reporterr(       "amino_dis (offset = %d): \n", offset );
 		for( i=0; i<20; i++ )
 		{
 			for( j=0; j<20; j++ ) 
 			{
-				fprintf( stderr, "%5d", amino_dis[(int)amino[i]][(int)amino[j]] );
+				reporterr(       "%5d", amino_dis[(int)amino[i]][(int)amino[j]] );
 			}
-			fprintf( stderr, "\n" );
+			reporterr(       "\n" );
 		}
 
-		fprintf( stderr, "amino_disLN (offsetLN = %d): \n", offsetLN );
+		reporterr(       "amino_disLN (offsetLN = %d): \n", offsetLN );
 		for( i=0; i<20; i++ )
 		{
 			for( j=0; j<20; j++ ) 
 			{
-				fprintf( stderr, "%5d", amino_disLN[(int)amino[i]][(int)amino[j]] );
+				reporterr(       "%5d", amino_disLN[(int)amino[i]][(int)amino[j]] );
 			}
-			fprintf( stderr, "\n" );
+			reporterr(       "\n" );
 		}
 
-		fprintf( stderr, "n_dis (offset = %d): \n", offset );
+		reporterr(       "n_dis (offset = %d): \n", offset );
 		for( i=0; i<26; i++ )
 		{
 			for( j=0; j<26; j++ ) 
 			{
-				fprintf( stderr, "%5d", n_dis[i][j] );
+				reporterr(       "%5d", n_dis[i][j] );
 			}
-			fprintf( stderr, "\n" );
+			reporterr(       "\n" );
 		}
 
-		fprintf( stderr, "n_disFFT (offsetFFT = %d): \n", offsetFFT );
+		reporterr(       "n_disFFT (offsetFFT = %d): \n", offsetFFT );
 		for( i=0; i<26; i++ )
 		{
 			for( j=0; j<26; j++ ) 
 			{
-				fprintf( stderr, "%5d", n_disFFT[i][j] );
+				reporterr(       "%5d", n_disFFT[i][j] );
 			}
-			fprintf( stderr, "\n" );
+			reporterr(       "\n" );
 		}
 exit( 1 );
 #endif
@@ -1035,3 +1463,13 @@ exit( 1 );
 #endif
 	}
 }
+
+void freeconstants()
+{
+	if( n_disLN ) FreeDoubleMtx( n_disLN ); n_disLN = NULL;
+	if( n_dis ) FreeIntMtx( n_dis ); n_dis = NULL;
+	if( n_disFFT ) FreeIntMtx( n_disFFT ); n_disFFT = NULL;
+	if( n_dis_consweight_multi ) FreeDoubleMtx( n_dis_consweight_multi ); n_dis_consweight_multi = NULL;
+	if( amino_dis ) FreeIntMtx( amino_dis ); amino_dis = NULL;
+	if( amino_dis_consweight_multi ) FreeDoubleMtx( amino_dis_consweight_multi ); amino_dis_consweight_multi = NULL;
+}