JWS-112 Bumping version of Mafft to version 7.310.
[jabaws.git] / binaries / src / mafft / core / Lalignmm.c
index 19f26ec..4346bf1 100644 (file)
 
 static int reccycle = 0;
 
-static float localthr;
+static double localthr;
 
-static void match_ribosum( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
+static void match_ribosum( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
 {
        int j, k, l;
-       float scarr[38];
-       float **cpmxpd = floatwork;
+       double scarr[38];
+       double **cpmxpd = doublework;
        int **cpmxpdn = intwork;
        int count = 0;
-       float *matchpt;
-       float **cpmxpdpt;
+       double *matchpt;
+       double **cpmxpdpt;
        int **cpmxpdnpt;
        int cpkd;
 
@@ -53,9 +53,9 @@ static void match_ribosum( float *match, float **cpmx1, float **cpmx2, int i1, i
                        scarr[l] += ribosumdis[k][l] * cpmx1[i1][k];
                }
        }
-#if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
+#if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïdoublework¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
        {
-               float *fpt, **fptpt, *fpt2;
+               double *fpt, **fptpt, *fpt2;
                int *ipt, **iptpt;
                fpt2 = match;
                iptpt = cpmxpdn;
@@ -91,24 +91,26 @@ static void match_ribosum( float *match, float **cpmx1, float **cpmx2, int i1, i
 #endif
 }
 
-static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
+static void match_calc( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
 {
        int j, k, l;
-       float scarr[26];
-       float **cpmxpd = floatwork;
+//     double scarr[26];
+       double **cpmxpd = doublework;
        int **cpmxpdn = intwork;
        int count = 0;
-       float *matchpt;
-       float **cpmxpdpt;
+       double *matchpt;
+       double **cpmxpdpt;
        int **cpmxpdnpt;
        int cpkd;
+       double *scarr;
+       scarr = calloc( nalphabets, sizeof( double ) );
 
        if( initialize )
        {
                for( j=0; j<lgth2; j++ )
                {
                        count = 0;
-                       for( l=0; l<26; l++ )
+                       for( l=0; l<nalphabets; l++ )
                        {
                                if( cpmx2[j][l] )
                                {
@@ -121,17 +123,17 @@ static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int
                }
        }
 
-       for( l=0; l<26; l++ )
+       for( l=0; l<nalphabets; l++ )
        {
                scarr[l] = 0.0;
-               for( k=0; k<26; k++ )
+               for( k=0; k<nalphabets; k++ )
                {
                        scarr[l] += (n_dis[k][l]-RNAthr) * cpmx1[i1][k];
                }
        }
-#if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
+#if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïdoublework¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
        {
-               float *fpt, **fptpt, *fpt2;
+               double *fpt, **fptpt, *fpt2;
                int *ipt, **iptpt;
                fpt2 = match;
                iptpt = cpmxpdn;
@@ -165,18 +167,19 @@ static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int
                cpmxpdpt++;
        }
 #endif
+       free( scarr );
 }
 
 #if 0
-static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
+static void match_add( double *match, double **cpmx1, double **cpmx2, int i1, int lgth2, double **doublework, int **intwork, int initialize )
 {
        int j, k, l;
-       float scarr[26];
-       float **cpmxpd = floatwork;
+       double scarr[nalphabets];
+       double **cpmxpd = doublework;
        int **cpmxpdn = intwork;
        int count = 0;
-       float *matchpt;
-       float **cpmxpdpt;
+       double *matchpt;
+       double **cpmxpdpt;
        int **cpmxpdnpt;
        int cpkd;
 
@@ -186,7 +189,7 @@ static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int l
                for( j=0; j<lgth2; j++ )
                {
                        count = 0;
-                       for( l=0; l<26; l++ )
+                       for( l=0; l<nalphabets; l++ )
                        {
                                if( cpmx2[j][l] )
                                {
@@ -199,17 +202,17 @@ static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int l
                }
        }
 
-       for( l=0; l<26; l++ )
+       for( l=0; l<nalphabets; l++ )
        {
                scarr[l] = 0.0;
-               for( k=0; k<26; k++ )
+               for( k=0; k<nalphabets; k++ )
                {
                        scarr[l] += n_dis[k][l] * cpmx1[i1][k];
                }
        }
-#if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïfloatwork¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
+#if 0 /* ¤³¤ì¤ò»È¤¦¤È¤\ad¤Ïdoublework¤Î¥¢¥í¥±¡¼¥È¤òµÕ¤Ë¤¹¤ë */
        {
-               float *fpt, **fptpt, *fpt2;
+               double *fpt, **fptpt, *fpt2;
                int *ipt, **iptpt;
                fpt2 = match;
                iptpt = cpmxpdn;
@@ -247,7 +250,7 @@ static void match_add( float *match, float **cpmx1, float **cpmx2, int i1, int l
 #endif
 
 #if 0
-static float Atracking( 
+static double Atracking( 
                                                char **seq1, char **seq2, 
                         char **mseq1, char **mseq2, 
                         int **ijp, int icyc, int jcyc,
@@ -340,46 +343,46 @@ static float Atracking(
 #endif
 
 
-static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, float **gapinfo, float **map )
+static double MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, double **cpmx1, double **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, double **gapinfo, double **map )
 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 {
-       float value = 0.0;
+       double value = 0.0;
        register int i, j;
        char **aseq1, **aseq2;
        int ll1, ll2;
        int lasti, lastj, imid, jmid = 0;
-       float wm = 0.0;   /* int ?????? */
-       float g;
-       float *currentw, *previousw;
+       double wm = 0.0;   /* int ?????? */
+       double g;
+       double *currentw, *previousw;
 #if USE_PENALTY_EX
-       float fpenalty_ex = (float)penalty_ex;
+       double fpenalty_ex = (double)penalty_ex;
 #endif
-//     float fpenalty = (float)penalty;
-       float *wtmp;
+//     double fpenalty = (double)penalty;
+       double *wtmp;
 //     short *ijppt;
        int *mpjpt;
 //     short **ijp;
        int *mp;
        int mpi;
-       float *mjpt, *prept, *curpt;
-       float mi;
-       float *m;
-       float *w1, *w2;
-//     float *match;
-       float *initverticalw;    /* kufuu sureba iranai */
-       float *lastverticalw;    /* kufuu sureba iranai */
+       double *mjpt, *prept, *curpt;
+       double mi;
+       double *m;
+       double *w1, *w2;
+//     double *match;
+       double *initverticalw;    /* kufuu sureba iranai */
+       double *lastverticalw;    /* kufuu sureba iranai */
        int **intwork;
-       float **floatwork;
+       double **doublework;
 //     short **shortmtx;
 #if STOREWM
-       float **WMMTX;
-       float **WMMTX2;
+       double **WMMTX;
+       double **WMMTX2;
 #endif
-       float *midw;
-       float *midm;
-       float *midn;
+       double *midw;
+       double *midm;
+       double *midn;
        int lgth1, lgth2;
-       float maxwm = 0.0;
+       double maxwm = 0.0;
        int *jumpforwi;
        int *jumpforwj;
        int *jumpbacki;
@@ -389,11 +392,11 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
        int jumpi, jumpj = 0;
        char *gaps;
        int ijpi, ijpj;
-       float *ogcp1;
-       float *fgcp1;
-       float *ogcp2;
-       float *fgcp2;
-       float firstm;
+       double *ogcp1;
+       double *fgcp1;
+       double *ogcp2;
+       double *fgcp2;
+       double firstm;
        int firstmp;
 #if 0
        static char ttt1[50000];
@@ -511,8 +514,8 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
        mp = AllocateIntVec( ll2+2 );
        gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
 
-       floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
-       intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); 
+       doublework = AllocateFloatMtx( MAX( ll1, ll2 )+2, nalphabets ); 
+       intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, nalphabets ); 
 
 #if DEBUG
        fprintf( stderr, "succeeded\n" );
@@ -535,9 +538,9 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
        currentw = w1;
        previousw = w2;
 
-       match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
+       match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, doublework, intwork, 1 );
 
-       match_ribosum( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
+       match_ribosum( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, doublework, intwork, 1 );
 
        for( i=1; i<lgth1+1; i++ )
        {
@@ -586,7 +589,7 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
 
                previousw[0] = initverticalw[i-1];
 
-               match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
+               match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
                currentw[0] = initverticalw[i];
 
                m[0] = ogcp1[i];
@@ -737,8 +740,8 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
 
 // gyakudp
 
-       match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
-       match_ribosum( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
+       match_ribosum( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, doublework, intwork, 1 );
+       match_ribosum( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, doublework, intwork, 1 );
 
        for( i=0; i<lgth1-1; i++ )
        {
@@ -790,7 +793,7 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
                currentw = wtmp;
                previousw[lgth2-1] = initverticalw[i+1];
 //             match_calc( currentw, seq1, seq2, i, lgth2 );
-               match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
+               match_ribosum( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
 
                currentw[lgth2-1] = initverticalw[i];
 
@@ -1054,8 +1057,8 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
 #if 0
        for( i=0; i<lgth1; i++ )
        {
-               float maxpairscore = -9999.9;
-               float tmpscore;
+               double maxpairscore = -9999.9;
+               double tmpscore;
 
                for( j=0; j<lgth2; j++ )
                {
@@ -1082,8 +1085,8 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
        }
        for( j=0; j<lgth2; j++ )
        {
-               float maxpairscore = -9999.9;
-               float tmpscore;
+               double maxpairscore = -9999.9;
+               double tmpscore;
 
                for( i=0; i<lgth1; i++ )
                {
@@ -1229,7 +1232,7 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
        FreeFloatVec( m );
        FreeIntVec( mp );
 
-       FreeFloatMtx( floatwork );
+       FreeFloatMtx( doublework );
        FreeIntMtx( intwork );
 
 #if STOREWM
@@ -1240,47 +1243,46 @@ static float MSalign2m2m_rec( int icyc, int jcyc, double *eff1, double *eff2, ch
        return( value );
 
 }
-static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, float **cpmx1, float **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, float **gapinfo, float **map )
+static double MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char **seq1, char **seq2, double **cpmx1, double **cpmx2, int ist, int ien, int jst, int jen, int alloclen, char **mseq1, char **mseq2, int depth, double **gapinfo, double **map )
 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 {
-       int alnlen;
-       float value = 0.0;
+       double value = 0.0;
        register int i, j;
        char **aseq1, **aseq2;
-       int ll1, ll2, l, len;
+       int ll1, ll2;
        int lasti, lastj, imid, jmid=0;
-       float wm = 0.0;   /* int ?????? */
-       float g;
-       float *currentw, *previousw;
+       double wm = 0.0;   /* int ?????? */
+       double g;
+       double *currentw, *previousw;
 #if USE_PENALTY_EX
-       float fpenalty_ex = (float)RNApenalty_ex;
+       double fpenalty_ex = (double)RNApenalty_ex;
 #endif
-//     float fpenalty = (float)penalty;
-       float *wtmp;
+//     double fpenalty = (double)penalty;
+       double *wtmp;
 //     short *ijppt;
        int *mpjpt;
 //     short **ijp;
        int *mp;
        int mpi;
-       float *mjpt, *prept, *curpt;
-       float mi;
-       float *m;
-       float *w1, *w2;
-//     float *match;
-       float *initverticalw;    /* kufuu sureba iranai */
-       float *lastverticalw;    /* kufuu sureba iranai */
+       double *mjpt, *prept, *curpt;
+       double mi;
+       double *m;
+       double *w1, *w2;
+//     double *match;
+       double *initverticalw;    /* kufuu sureba iranai */
+       double *lastverticalw;    /* kufuu sureba iranai */
        int **intwork;
-       float **floatwork;
+       double **doublework;
 //     short **shortmtx;
 #if STOREWM
-       float **WMMTX;
-       float **WMMTX2;
+       double **WMMTX;
+       double **WMMTX2;
 #endif
-       float *midw;
-       float *midm;
-       float *midn;
+       double *midw;
+       double *midm;
+       double *midn;
        int lgth1, lgth2;
-       float maxwm = 0.0;
+       double maxwm = 0.0;
        int *jumpforwi;
        int *jumpforwj;
        int *jumpbacki;
@@ -1290,11 +1292,11 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
        int jumpi, jumpj = 0;
        char *gaps;
        int ijpi, ijpj;
-       float *ogcp1;
-       float *fgcp1;
-       float *ogcp2;
-       float *fgcp2;
-       float firstm;
+       double *ogcp1;
+       double *fgcp1;
+       double *ogcp2;
+       double *fgcp2;
+       double firstm;
        int firstmp;
 #if 0
        static char ttt1[50000];
@@ -1412,8 +1414,8 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
        mp = AllocateIntVec( ll2+2 );
        gaps = AllocateCharVec( MAX( ll1, ll2 ) + 2 );
 
-       floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 ); 
-       intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 ); 
+       doublework = AllocateFloatMtx( MAX( ll1, ll2 )+2, nalphabets ); 
+       intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, nalphabets ); 
 
 #if DEBUG
        fprintf( stderr, "succeeded\n" );
@@ -1436,9 +1438,9 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
        currentw = w1;
        previousw = w2;
 
-       match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, floatwork, intwork, 1 );
+       match_calc( initverticalw, cpmx2+jst, cpmx1+ist, 0, lgth1, doublework, intwork, 1 );
 
-       match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, floatwork, intwork, 1 );
+       match_calc( currentw, cpmx1+ist, cpmx2+jst, 0, lgth2, doublework, intwork, 1 );
 
        for( i=1; i<lgth1+1; i++ )
        {
@@ -1487,7 +1489,7 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
 
                previousw[0] = initverticalw[i-1];
 
-               match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
+               match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
                currentw[0] = initverticalw[i];
 
                m[0] = ogcp1[i];
@@ -1638,8 +1640,8 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
 
 // gyakudp
 
-       match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, floatwork, intwork, 1 );
-       match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, floatwork, intwork, 1 );
+       match_calc( initverticalw, cpmx2+jst, cpmx1+ist, lgth2-1, lgth1, doublework, intwork, 1 );
+       match_calc( currentw, cpmx1+ist, cpmx2+jst, lgth1-1, lgth2, doublework, intwork, 1 );
 
        for( i=0; i<lgth1-1; i++ )
        {
@@ -1691,7 +1693,7 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
                currentw = wtmp;
                previousw[lgth2-1] = initverticalw[i+1];
 //             match_calc( currentw, seq1, seq2, i, lgth2 );
-               match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, floatwork, intwork, 0 );
+               match_calc( currentw, cpmx1+ist, cpmx2+jst, i, lgth2, doublework, intwork, 0 );
 
                currentw[lgth2-1] = initverticalw[i];
 
@@ -1955,8 +1957,8 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
 #if 0
        for( i=0; i<lgth1; i++ )
        {
-               float maxpairscore = -9999.9;
-               float tmpscore;
+               double maxpairscore = -9999.9;
+               double tmpscore;
 
                for( j=0; j<lgth2; j++ )
                {
@@ -1983,8 +1985,8 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
        }
        for( j=0; j<lgth2; j++ )
        {
-               float maxpairscore = -9999.9;
-               float tmpscore;
+               double maxpairscore = -9999.9;
+               double tmpscore;
 
                for( i=0; i<lgth1; i++ )
                {
@@ -2132,7 +2134,7 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
        FreeFloatVec( m );
        FreeIntVec( mp );
 
-       FreeFloatMtx( floatwork );
+       FreeFloatMtx( doublework );
        FreeIntMtx( intwork );
 
 #if STOREWM
@@ -2140,130 +2142,6 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
        FreeFloatMtx( WMMTX2 );
 #endif
 
-       return( value );
-
-//     fprintf( stderr, "==== calling myself (first)\n" );
-
-#if 0
-               fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
-               fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
-#endif
-       value = MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist, ist+jumpi, jst, jst+jumpj, alloclen, aseq1, aseq2, depth, gapinfo, map ); 
-#if 0
-               fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
-               fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
-#endif
-#if MEMSAVE
-#else
-       for( i=0; i<icyc; i++ ) strcpy( mseq1[i], aseq1[i] );
-       for( i=0; i<jcyc; i++ ) strcpy( mseq2[i], aseq2[i] );
-#endif
-
-//     fprintf( stderr, "====(f) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
-//     fprintf( stderr, "====(f) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
-
-       len = strlen( mseq1[0] );
-//     fprintf( stderr, "len = %d\n", len );
-       l = jmid - jumpj - 1;
-//     fprintf( stderr, "l=%d\n", l );
-       if( l > 0 )
-       {
-               for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
-               for( i=0; i<icyc; i++ ) 
-               {
-                       strcat( mseq1[i], gaps );
-                       mseq1[i][len+l] = 0;
-               }
-               for( j=0; j<jcyc; j++ )
-               {
-                       strncat( mseq2[j], seq2[j]+jst+jumpj+1, l );
-                       mseq2[j][len+l] = 0;
-               }
-//             fprintf( stderr, "penalizing (2) .. %f(%d), %f(%d)\n", ogcp2[jumpj+1], jumpj+1, fgcp2[jmid-1], jmid-1 );
-               value +=  ( ogcp2[jumpj+1] + fgcp2[jmid-1] );
-//             value += fpenalty;
-       }
-       len = strlen( mseq1[0] );
-       l = imid - jumpi - 1;
-//     fprintf( stderr, "l=%d\n", l );
-       if( l > 0 )
-       {
-               for( i=0; i<l; i++ ) gaps[i] = '-'; gaps[i] = 0;
-               for( i=0; i<icyc; i++ )
-               {
-                       strncat( mseq1[i], seq1[i]+ist+jumpi+1, l );
-                       mseq1[i][len+l] = 0;
-               }
-               for( j=0; j<jcyc; j++ ) 
-               {
-                       strcat( mseq2[j], gaps );
-                       mseq2[j][len+l] = 0;
-               }
-
-//             for( i=0; i<lgth1; i++ ) fprintf( stderr, "ogcp1[%d] = %f\n", i, ogcp1[i] );
-//             for( i=0; i<lgth1; i++ ) fprintf( stderr, "fgcp1[%d] = %f\n", i, fgcp1[i] );
-
-
-//             fprintf( stderr, "penalizing (1) .. ogcp1[%d] = %f, fgcp1[%d] = %f\n", jumpi+1, ogcp1[jumpi+1], imid-1, fgcp1[imid-1] );
-               value += ( ogcp1[jumpi+1] + fgcp1[imid-1] );
-//             value += fpenalty;
-       }
-#if 0
-       for( i=0; i<icyc; i++ ) fprintf( stderr, "after gapfill mseq1[%d]=%s\n", i, mseq1[i] );
-       for( i=0; i<jcyc; i++ ) fprintf( stderr, "after gapfill mseq2[%d]=%s\n", i, mseq2[i] );
-#endif
-
-//     fprintf( stderr, "==== calling myself (second)\n" );
-
-#if MEMSAVE
-       alnlen = strlen( aseq1[0] );
-       for( i=0; i<icyc; i++ ) aseq1[i] += alnlen;
-       for( i=0; i<jcyc; i++ ) aseq2[i] += alnlen;
-#endif
-
-#if 0
-               fprintf( stderr, "seq1[0] = %s\n", seq1[0] );
-               fprintf( stderr, "seq2[0] = %s\n", seq2[0] );
-#endif
-       value += MSalignmm_rec( icyc, jcyc, eff1, eff2, seq1, seq2, cpmx1, cpmx2, ist+imid, ien, jst+jmid, jen, alloclen, aseq1, aseq2, depth, gapinfo, map );  
-#if 0
-               fprintf( stderr, "aseq1[0] = %s\n", aseq1[0] );
-               fprintf( stderr, "aseq2[0] = %s\n", aseq2[0] );
-#endif
-
-
-
-#if DEBUG
-       if( value - maxwm > 1 || maxwm - value > 1 )
-       {
-               fprintf( stderr, "WARNING value  = %f, but maxwm = %f\n", value, maxwm );
-               for( i=0; i<icyc; i++ )
-               {
-                       fprintf( stderr, ">1-%d\n%s\n", i, mseq1[i] );
-                       fprintf( stderr, "%s\n", aseq1[i] );
-               }
-               for( i=0; i<jcyc; i++ )
-               {
-                       fprintf( stderr, ">2-%d\n%s\n", i, mseq2[i] );
-                       fprintf( stderr, "%s\n", aseq2[i] );
-               }
-
-//             exit( 1 );
-       }
-       else
-       {
-               fprintf( stderr, "value = %.0f, maxwm = %.0f -> ok\n", value, maxwm );
-       }
-#endif
-
-#if MEMSAVE
-#else
-       for( i=0; i<icyc; i++ ) strcat( mseq1[i], aseq1[i] );
-       for( i=0; i<jcyc; i++ ) strcat( mseq2[i], aseq2[i] );
-#endif
-
-//     fprintf( stderr, "====(s) aseq1[0] (%d) = %s (%d-%d)\n", depth, aseq1[0], ist, ien );
-//     fprintf( stderr, "====(s) aseq2[0] (%d) = %s (%d-%d)\n", depth, aseq2[0], jst, jen );
 
        free( gaps );
 #if MEMSAVE
@@ -2279,26 +2157,26 @@ static float MSalignmm_rec( int icyc, int jcyc, double *eff1, double *eff2, char
 
 
 
-float Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, float **map )
+double Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, double **map )
 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 {
 //     int k;
        int i, j;
        int ll1, ll2;
        int lgth1, lgth2;
-       float wm = 0.0;   /* int ?????? */
+       double wm = 0.0;   /* int ?????? */
        char **mseq1;
        char **mseq2;
 //     char **mseq;
-       float *ogcp1;
-       float *ogcp2;
-       float *fgcp1;
-       float *fgcp2;
-       float **cpmx1;
-       float **cpmx2;
-       float **gapinfo;
-//     float fpenalty;
-       float fpenalty = (float)RNApenalty;
+       double *ogcp1;
+       double *ogcp2;
+       double *fgcp1;
+       double *fgcp2;
+       double **cpmx1;
+       double **cpmx2;
+       double **gapinfo;
+//     double fpenalty;
+       double fpenalty = (double)RNApenalty;
        int nglen1, nglen2;
 
 
@@ -2329,8 +2207,8 @@ float Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int
        fgcp2 = AllocateFloatVec( ll2+2 );
 
 
-       cpmx1 = AllocateFloatMtx( ll1+2, 27 );
-       cpmx2 = AllocateFloatMtx( ll2+2, 27 );
+       cpmx1 = AllocateFloatMtx( ll1+2, nalphabets+1 );
+       cpmx2 = AllocateFloatMtx( ll2+2, nalphabets+1 );
 
        for( i=0; i<icyc; i++ ) 
        {
@@ -2482,24 +2360,24 @@ float Lalignmm_hmout( char **seq1, char **seq2, double *eff1, double *eff2, int
        return( wm );
 }
 
-float Lalign2m2m_hmout( char **seq1, char **seq2, char **seq1r, char **seq2r, char *dir1, char *dir2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, float **map )
+double Lalign2m2m_hmout( char **seq1, char **seq2, char **seq1r, char **seq2r, char *dir1, char *dir2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, char *sgap1, char *sgap2, char *egap1, char *egap2, double **map )
 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 {
 //     int k;
        int i, j;
        int ll1, ll2;
        int lgth1, lgth2;
-       float wm = 0.0;   /* int ?????? */
+       double wm = 0.0;   /* int ?????? */
        char **mseq1;
        char **mseq2;
-       float *ogcp1;
-       float *ogcp2;
-       float *fgcp1;
-       float *fgcp2;
-       float **cpmx1;
-       float **cpmx2;
-       float **gapinfo;
-       float fpenalty = (float)penalty;
+       double *ogcp1;
+       double *ogcp2;
+       double *fgcp1;
+       double *fgcp2;
+       double **cpmx1;
+       double **cpmx2;
+       double **gapinfo;
+       double fpenalty = (double)penalty;
        int nglen1, nglen2;
 
 #if 0