new mafft v 6.857 with extensions
[jabaws.git] / binaries / src / mafft / core / Salignmm.c
index 1170040..06e402e 100644 (file)
@@ -2,15 +2,15 @@
 #include "dp.h"
 
 #define MACHIGAI 0
-#define OUTGAP0TRY 1
+#define OUTGAP0TRY 0
 #define DEBUG 0
 #define XXXXXXX    0
 #define USE_PENALTY_EX  0
 #define FASTMATCHCALC 1
 
 
-static float **impmtx = NULL;
-static int impalloclen = 0;
+static TLS float **impmtx = NULL;
+static TLS int impalloclen = 0;
 float imp_match_out_sc( int i1, int j1 )
 {
 //     fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
@@ -79,10 +79,22 @@ void imp_match_init_strict( float *imp, int clus1, int clus2, int lgth1, int lgt
        float effij_kozo;
        double effijx;
        char *pt, *pt1, *pt2;
-       static char *nocount1 = NULL;
-       static char *nocount2 = NULL;
+       static TLS char *nocount1 = NULL;
+       static TLS char *nocount2 = NULL;
        LocalHom *tmpptr;
 
+       if( seq1 == NULL )
+       {
+               if( impmtx ) FreeFloatMtx( impmtx );
+               impmtx = NULL;
+               if( nocount1 ) free( nocount1 );
+               nocount1 = NULL;
+               if( nocount2 ) free( nocount2 );
+               nocount2 = NULL;
+               
+               return;
+       }
+
        if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
        {
                if( impmtx ) FreeFloatMtx( impmtx );
@@ -295,11 +307,11 @@ for( i = 0; i<impalloclen; i++ )
 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
 {
        int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
-       static int impalloclen = 0;
+       static TLS int impalloclen = 0;
        char *pt;
        int allgap;
-       static char *nocount1 = NULL;
-       static char *nocount2 = NULL;
+       static TLS char *nocount1 = NULL;
+       static TLS char *nocount2 = NULL;
 
        if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
        {
@@ -761,7 +773,8 @@ static float Atracking( float *lasthorizontalw, float *lastverticalw,
                                                char **seq1, char **seq2, 
                         char **mseq1, char **mseq2, 
                         float **cpmx1, float **cpmx2, 
-                        int **ijp, int icyc, int jcyc )
+                        int **ijp, int icyc, int jcyc,
+                                               int tailgp )
 {
        int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
        float wm;
@@ -780,7 +793,7 @@ static float Atracking( float *lasthorizontalw, float *lastverticalw,
        }
 #endif
  
-       if( outgap == 1 )
+       if( tailgp == 1 )
                ;
        else
        {
@@ -864,7 +877,7 @@ static float Atracking( float *lasthorizontalw, float *lastverticalw,
        return( 0.0 );
 }
 
-float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2 )
+float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp )
 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 {
 //     int k;
@@ -885,25 +898,25 @@ float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc,
        float *mjpt, *prept, *curpt;
        int *mpjpt;
 #endif
-       static float mi, *m;
-       static int **ijp;
-       static int mpi, *mp;
-       static float *w1, *w2;
-       static float *match;
-       static float *initverticalw;    /* kufuu sureba iranai */
-       static float *lastverticalw;    /* kufuu sureba iranai */
-       static char **mseq1;
-       static char **mseq2;
-       static char **mseq;
-       static float *ogcp1;
-       static float *ogcp2;
-       static float *fgcp1;
-       static float *fgcp2;
-       static float **cpmx1;
-       static float **cpmx2;
-       static int **intwork;
-       static float **floatwork;
-       static int orlgth1 = 0, orlgth2 = 0;
+       static TLS float mi, *m;
+       static TLS int **ijp;
+       static TLS int mpi, *mp;
+       static TLS float *w1, *w2;
+       static TLS float *match;
+       static TLS float *initverticalw;    /* kufuu sureba iranai */
+       static TLS float *lastverticalw;    /* kufuu sureba iranai */
+       static TLS char **mseq1;
+       static TLS char **mseq2;
+       static TLS char **mseq;
+       static TLS float *ogcp1;
+       static TLS float *ogcp2;
+       static TLS float *fgcp1;
+       static TLS float *fgcp2;
+       static TLS float **cpmx1;
+       static TLS float **cpmx2;
+       static TLS int **intwork;
+       static TLS float **floatwork;
+       static TLS int orlgth1 = 0, orlgth2 = 0;
        float fpenalty = (float)penalty;
        float *fgcp2pt;
        float *ogcp2pt;
@@ -911,12 +924,93 @@ float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc,
        float ogcp1va;
 
 
+       if( seq1 == NULL )
+       {
+               if( orlgth1 )
+               {
+//                     fprintf( stderr, "## Freeing local arrays in A__align\n" );
+                       orlgth1 = 0;
+                       orlgth2 = 0;
+
+                       imp_match_init_strict( NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0 );
+
+                       free( mseq1 );
+                       free( mseq2 );
+                       FreeFloatVec( w1 );
+                       FreeFloatVec( w2 );
+                       FreeFloatVec( match );
+                       FreeFloatVec( initverticalw );
+                       FreeFloatVec( lastverticalw );
+
+                       FreeFloatVec( m );
+                       FreeIntVec( mp );
+
+                       FreeCharMtx( mseq );
+
+                       FreeFloatVec( ogcp1 );
+                       FreeFloatVec( ogcp2 );
+                       FreeFloatVec( fgcp1 );
+                       FreeFloatVec( fgcp2 );
+
+
+                       FreeFloatMtx( cpmx1 );
+                       FreeFloatMtx( cpmx2 );
+
+                       FreeFloatMtx( floatwork );
+                       FreeIntMtx( intwork );
+
+               }
+               else
+               {
+//                     fprintf( stderr, "## Not allocated\n" );
+               }
+               return( 0.0 );
+       }
+
+       lgth1 = strlen( seq1[0] );
+       lgth2 = strlen( seq2[0] );
+#if 1
+       if( lgth1 == 0 || lgth2 == 0 )
+       {
+               fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
+       }
+       if( lgth1 == 0 && lgth2 == 0 )
+               return( 0.0 );
+
+       if( lgth1 == 0 )
+       {
+               for( i=0; i<icyc; i++ )
+               {
+                       j = lgth2;
+                       seq1[i][j] = 0;
+                       while( j ) seq1[i][--j] = '-';
+//                     fprintf( stderr, "seq1[i] = %s\n", seq1[i] );
+               }
+               return( 0.0 );
+       }
+
+       if( lgth2 == 0 )
+       {
+               for( i=0; i<jcyc; i++ )
+               {
+                       j = lgth1;
+                       seq2[i][j] = 0;
+                       while( j ) seq2[i][--j] = '-';
+//                     fprintf( stderr, "seq2[i] = %s\n", seq2[i] );
+               }
+               return( 0.0 );
+       }
+#endif
+
 
 #if 0
        fprintf( stderr, "####  eff in SA+++align\n" );
        fprintf( stderr, "####  seq1[0] = %s\n", seq1[0] );
        fprintf( stderr, "####  strlen( seq1[0] ) = %d\n", strlen( seq1[0] ) );
        for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
+       fprintf( stderr, "####  seq2[0] = %s\n", seq2[0] );
+       fprintf( stderr, "####  strlen( seq2[0] ) = %d\n", strlen( seq2[0] ) );
+       for( i=0; i<jcyc; i++ ) fprintf( stderr, "eff2[%d] = %f\n", i, eff2[i] );
 #endif
        if( orlgth1 == 0 )
        {
@@ -925,14 +1019,6 @@ float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc,
        }
 
 
-       lgth1 = strlen( seq1[0] );
-       lgth2 = strlen( seq2[0] );
-#if 0
-       if( lgth1 == 0 || lgth2 == 0 )
-       {
-               fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
-       }
-#endif
 
        if( lgth1 > orlgth1 || lgth2 > orlgth2 )
        {
@@ -1105,7 +1191,7 @@ float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc,
 
 #endif
 
-       if( outgap == 1 )
+       if( headgp == 1 )
        {
                for( i=1; i<lgth1+1; i++ )
                {
@@ -1135,7 +1221,7 @@ float A__align( char **seq1, char **seq2, double *eff1, double *eff2, int icyc,
        else
                lastverticalw[0] = currentw[lgth2-1];
 
-       if( outgap ) lasti = lgth1+1; else lasti = lgth1;
+       if( tailgp ) lasti = lgth1+1; else lasti = lgth1;
 
 #if XXXXXXX
 fprintf( stderr, "currentw = \n" );
@@ -1159,6 +1245,16 @@ for( i=0; i<lgth2; i++ )
 
        for( i=1; i<lasti; i++ )
        {
+
+#ifdef enablemultithread
+//             fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
+               if( chudanpt && *chudanpt != chudanref ) 
+               {
+//                     fprintf( stderr, "\n\n## CHUUDAN!!! S\n" );
+                       *chudanres = 1;
+                       return( -1.0 );
+               }
+#endif
                wtmp = previousw; 
                previousw = currentw;
                currentw = wtmp;
@@ -1211,6 +1307,15 @@ fprintf( stderr, "\n" );
                lastj = lgth2+1;
                for( j=1; j<lastj; j++ )
                {
+#ifdef xxxenablemultithread
+//                     fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref );
+                       if( chudanpt && *chudanpt != chudanref ) 
+                       {
+//                             fprintf( stderr, "\n\n## CHUUDAN!!! S\n" );
+                               *chudanres = 1;
+                               return( -1.0 );
+                       }
+#endif
                        wm = *prept;
                        *ijppt = 0;
 
@@ -1291,7 +1396,7 @@ fprintf( stderr, "\n" );
                Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
        }
        else
-               Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
+               Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, tailgp );
 
 //     fprintf( stderr, "### impmatch = %f\n", *impmatch );
 
@@ -1305,12 +1410,12 @@ fprintf( stderr, "\n" );
 
        for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
        for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
-       /*
+#if 0
        fprintf( stderr, "\n" );
        for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
        fprintf( stderr, "#####\n" );
        for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
-       */
+#endif
 
 //     fprintf( stderr, "wm = %f\n", wm );
 
@@ -1338,25 +1443,25 @@ float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int
        float *mjpt, *prept, *curpt;
        int *mpjpt;
 #endif
-       static float mi, *m;
-       static int **ijp;
-       static int mpi, *mp;
-       static float *w1, *w2;
-       static float *match;
-       static float *initverticalw;    /* kufuu sureba iranai */
-       static float *lastverticalw;    /* kufuu sureba iranai */
-       static char **mseq1;
-       static char **mseq2;
-       static char **mseq;
-       static float *ogcp1;
-       static float *ogcp2;
-       static float *fgcp1;
-       static float *fgcp2;
-       static float **cpmx1;
-       static float **cpmx2;
-       static int **intwork;
-       static float **floatwork;
-       static int orlgth1 = 0, orlgth2 = 0;
+       static TLS float mi, *m;
+       static TLS int **ijp;
+       static TLS int mpi, *mp;
+       static TLS float *w1, *w2;
+       static TLS float *match;
+       static TLS float *initverticalw;    /* kufuu sureba iranai */
+       static TLS float *lastverticalw;    /* kufuu sureba iranai */
+       static TLS char **mseq1;
+       static TLS char **mseq2;
+       static TLS char **mseq;
+       static TLS float *ogcp1;
+       static TLS float *ogcp2;
+       static TLS float *fgcp1;
+       static TLS float *fgcp2;
+       static TLS float **cpmx1;
+       static TLS float **cpmx2;
+       static TLS int **intwork;
+       static TLS float **floatwork;
+       static TLS int orlgth1 = 0, orlgth2 = 0;
        float *fgcp2pt;
        float *ogcp2pt;
        float fgcp1va;
@@ -1537,7 +1642,7 @@ float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int
 
 #endif
 
-       if( outgap == 1 )
+       if( 1 ) // tsuneni outgap=1
        {
                for( i=1; i<lgth1+1; i++ )
                {
@@ -1568,7 +1673,7 @@ float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int
        else
                lastverticalw[0] = currentw[lgth2-1];
 
-       if( outgap ) lasti = lgth1+1; else lasti = lgth1;
+       if( 1 ) lasti = lgth1+1; else lasti = lgth1; // tsuneni outgap=1
 
 #if XXXXXXX
 fprintf( stderr, "currentw = \n" );
@@ -1727,7 +1832,7 @@ fprintf( stderr, "\n" );
                Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
        }
        else
-               Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
+               Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, 1 );
 
 //     fprintf( stderr, "### impmatch = %f\n", *impmatch );
 
@@ -1753,40 +1858,3 @@ fprintf( stderr, "\n" );
        return( wm );
 }
 
-float translate_and_Calign( char **mseq1, char **mseq2, double *effarr1, double *effarr2, int clus1, int clus2, int alloclen )
-{
-    int i;
-    float wm;
-    char **result;
-    char *seq = NULL, **aseq = NULL;
-    double *effarr = NULL;
-    int nseq = 0;
-       int resultlen;
-
-    if     ( clus1 == 1 && clus2 != 1 ) 
-    {
-        seq = mseq1[0]; aseq = mseq2; effarr = effarr2; nseq = clus2+1;
-#if 0
-               printf( "effarr in transl... = \n" );
-               for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
-#endif
-    }
-    else if( clus1 != 1 && clus2 == 1 ) 
-    {
-        seq = mseq2[0]; aseq = mseq1; effarr = effarr1; nseq = clus1+1;
-    }
-    else ErrorExit( "ERROR in translate_and_Calign" );
-
-    result = Calignm1( &wm, aseq, seq, effarr, nseq-2, 0 );
-
-       resultlen = strlen( result[0] );
-       if( alloclen < resultlen || resultlen > N )
-       {
-               fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
-               ErrorExit( "LENGTH OVER!\n" );
-       }
-    for( i=0; i<nseq-1; i++ ) strcpy( aseq[i], result[i] );
-    strcpy( seq, result[nseq-1] );
-
-    return( 0.0 );
-}