binaries/src/mafft/core/genGalign11.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define DEBUG 0
   5 #define XXXXXXX    0
   6 #define USE_PENALTY_EX  1
   7
   8
   9 #if 1
  10 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
  11 {
  12         char *seq2 = s2[0];
  13         int *intptr = amino_dis[(int)s1[0][i1]];
  14
  15         while( lgth2-- )
  16                 *match++ = intptr[(int)*seq2++];
  17 }
  18 #else
  19 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
  20 {
  21         int j;
  22
  23         for( j=0; j<lgth2; j++ )
  24                 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
  25 }
  26 #endif
  27
  28 static float genGtracking( float *lasthorizontalw, float *lastverticalw,
  29                                                 char **seq1, char **seq2,
  30                         char **mseq1, char **mseq2,
  31                         float **cpmx1, float **cpmx2,
  32                         int **ijpi, int **ijpj )
  33 {
  34         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
  35 //      char gap[] = "-";
  36         char *gap;
  37         gap = newgapstr;
  38         lgth1 = strlen( seq1[0] );
  39         lgth2 = strlen( seq2[0] );
  40
  41
  42 #if 0
  43         for( i=0; i<lgth1; i++ )
  44         {
  45                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
  46         }
  47 #endif
  48
  49     for( i=0; i<lgth1+1; i++ )
  50     {
  51         ijpi[i][0] = -1;
  52                 ijpj[i][0] = -1; // ???
  53     }
  54     for( j=0; j<lgth2+1; j++ )
  55     {
  56                 ijpi[0][j] = -1; // ???
  57         ijpj[0][j] = -1;
  58     }
  59
  60
  61         mseq1[0] += lgth1+lgth2;
  62         *mseq1[0] = 0;
  63         mseq2[0] += lgth1+lgth2;
  64         *mseq2[0] = 0;
  65         iin = lgth1; jin = lgth2;
  66         limk = lgth1+lgth2 + 1;
  67         for( k=0; k<limk; k++ )
  68         {
  69                 ifi = ( ijpi[iin][jin] );
  70                 jfi = ( ijpj[iin][jin] );
  71                 l = iin - ifi;
  72                 while( --l )
  73                 {
  74                         *--mseq1[0] = seq1[0][ifi+l];
  75                         *--mseq2[0] = *gap;
  76                         k++;
  77                 }
  78                 l= jin - jfi;
  79                 while( --l )
  80                 {
  81                         *--mseq1[0] = *gap;
  82                         *--mseq2[0] = seq2[0][jfi+l];
  83                         k++;
  84                 }
  85                 if( iin <= 0 || jin <= 0 ) break;
  86                 *--mseq1[0] = seq1[0][ifi];
  87                 *--mseq2[0] = seq2[0][jfi];
  88
  89 //              fprintf( stdout, "mseq1 = %s\n", mseq1[0] );
  90 //              fprintf( stdout, "mseq2 = %s\n", mseq2[0] );
  91
  92                 k++;
  93                 iin = ifi; jin = jfi;
  94         }
  95         return( 0.0 );
  96 }
  97
  98
  99 float genG__align11( char **seq1, char **seq2, int alloclen )
 100 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 101 {
 102 //      int k;
 103         register int i, j;
 104         int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 105         int lgth1, lgth2;
 106         int resultlen;
 107         float wm;   /* int ?????? */
 108         float g;
 109         float *currentw, *previousw;
 110         float fpenalty = (float)penalty;
 111         float fpenalty_OP = (float)penalty_OP;
 112 #if USE_PENALTY_EX
 113         float fpenalty_ex = (float)penalty_ex;
 114 #endif
 115 #if 1
 116         float *wtmp;
 117         int *ijpipt;
 118         int *ijpjpt;
 119         float *mjpt, *Mjpt, *prept, *curpt;
 120         int *mpjpt, *Mpjpt;
 121 #endif
 122         static float mi, *m;
 123         static float Mi, *largeM;
 124         static int **ijpi;
 125         static int **ijpj;
 126         static int mpi, *mp;
 127         static int Mpi, *Mp;
 128         static float *w1, *w2;
 129         static float *match;
 130         static float *initverticalw;    /* kufuu sureba iranai */
 131         static float *lastverticalw;    /* kufuu sureba iranai */
 132         static char **mseq1;
 133         static char **mseq2;
 134         static char **mseq;
 135         static float **cpmx1;
 136         static float **cpmx2;
 137         static int **intwork;
 138         static float **floatwork;
 139         static int orlgth1 = 0, orlgth2 = 0;
 140         float tbk;
 141         int tbki, tbkj;
 142
 143         wm = 0.0;
 144
 145         if( orlgth1 == 0 )
 146         {
 147                 mseq1 = AllocateCharMtx( njob, 0 );
 148                 mseq2 = AllocateCharMtx( njob, 0 );
 149         }
 150
 151
 152         lgth1 = strlen( seq1[0] );
 153         lgth2 = strlen( seq2[0] );
 154
 155
 156
 157         if( lgth1 <= 0 || lgth2 <= 0 )
 158         {
 159                 fprintf( stderr, "WARNING (g11): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 160         }
 161
 162         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 163         {
 164                 int ll1, ll2;
 165
 166                 if( orlgth1 > 0 && orlgth2 > 0 )
 167                 {
 168                         FreeFloatVec( w1 );
 169                         FreeFloatVec( w2 );
 170                         FreeFloatVec( match );
 171                         FreeFloatVec( initverticalw );
 172                         FreeFloatVec( lastverticalw );
 173
 174                         FreeFloatVec( m );
 175                         FreeIntVec( mp );
 176                         FreeFloatVec( largeM );
 177                         FreeIntVec( Mp );
 178
 179                         FreeCharMtx( mseq );
 180
 181
 182                         FreeFloatMtx( cpmx1 );
 183                         FreeFloatMtx( cpmx2 );
 184
 185                         FreeFloatMtx( floatwork );
 186                         FreeIntMtx( intwork );
 187                 }
 188
 189                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 190                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 191
 192 #if DEBUG
 193                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 194 #endif
 195
 196                 w1 = AllocateFloatVec( ll2+2 );
 197                 w2 = AllocateFloatVec( ll2+2 );
 198                 match = AllocateFloatVec( ll2+2 );
 199
 200                 initverticalw = AllocateFloatVec( ll1+2 );
 201                 lastverticalw = AllocateFloatVec( ll1+2 );
 202
 203                 m = AllocateFloatVec( ll2+2 );
 204                 mp = AllocateIntVec( ll2+2 );
 205                 largeM = AllocateFloatVec( ll2+2 );
 206                 Mp = AllocateIntVec( ll2+2 );
 207
 208                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 209
 210                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 211                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 212
 213                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 214                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 215
 216 #if DEBUG
 217                 fprintf( stderr, "succeeded\n" );
 218 #endif
 219
 220                 orlgth1 = ll1 - 100;
 221                 orlgth2 = ll2 - 100;
 222         }
 223
 224
 225         mseq1[0] = mseq[0];
 226         mseq2[0] = mseq[1];
 227
 228
 229         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 230         {
 231                 int ll1, ll2;
 232
 233                 if( commonAlloc1 && commonAlloc2 )
 234                 {
 235                         FreeIntMtx( commonIP );
 236                         FreeIntMtx( commonJP );
 237                 }
 238
 239                 ll1 = MAX( orlgth1, commonAlloc1 );
 240                 ll2 = MAX( orlgth2, commonAlloc2 );
 241
 242 #if DEBUG
 243                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 244 #endif
 245
 246                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
 247                 commonJP = AllocateIntMtx( ll1+10, ll2+10 );
 248
 249 #if DEBUG
 250                 fprintf( stderr, "succeeded\n\n" );
 251 #endif
 252
 253                 commonAlloc1 = ll1;
 254                 commonAlloc2 = ll2;
 255         }
 256         ijpi = commonIP;
 257         ijpj = commonJP;
 258
 259
 260 #if 0
 261         for( i=0; i<lgth1; i++ )
 262                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
 263 #endif
 264
 265         currentw = w1;
 266         previousw = w2;
 267
 268
 269         match_calc( initverticalw, seq2, seq1, 0, lgth1 );
 270
 271
 272         match_calc( currentw, seq1, seq2, 0, lgth2 );
 273
 274         if( outgap == 1 )
 275         {
 276                 for( i=1; i<lgth1+1; i++ )
 277                 {
 278                         initverticalw[i] += fpenalty;
 279                 }
 280                 for( j=1; j<lgth2+1; j++ )
 281                 {
 282                         currentw[j] += fpenalty;
 283                 }
 284         }
 285
 286         for( j=1; j<lgth2+1; ++j )
 287         {
 288                 m[j] = currentw[j-1]; mp[j] = 0;
 289                 largeM[j] = currentw[j-1]; Mp[j] = 0;
 290         }
 291
 292         if( lgth2 == 0 )
 293                 lastverticalw[0] = 0.0;               // lgth2==0 no toki error
 294         else
 295                 lastverticalw[0] = currentw[lgth2-1]; // lgth2==0 no toki error
 296
 297         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 298
 299 #if XXXXXXX
 300 fprintf( stderr, "currentw = \n" );
 301 for( i=0; i<lgth1+1; i++ )
 302 {
 303         fprintf( stderr, "%5.2f ", currentw[i] );
 304 }
 305 fprintf( stderr, "\n" );
 306 fprintf( stderr, "initverticalw = \n" );
 307 for( i=0; i<lgth2+1; i++ )
 308 {
 309         fprintf( stderr, "%5.2f ", initverticalw[i] );
 310 }
 311 fprintf( stderr, "\n" );
 312 #endif
 313
 314         for( i=1; i<lasti; i++ )
 315         {
 316                 wtmp = previousw;
 317                 previousw = currentw;
 318                 currentw = wtmp;
 319
 320                 previousw[0] = initverticalw[i-1];
 321
 322                 match_calc( currentw, seq1, seq2, i, lgth2 );
 323 #if XXXXXXX
 324 fprintf( stderr, "\n" );
 325 fprintf( stderr, "i=%d\n", i );
 326 fprintf( stderr, "currentw = \n" );
 327 for( j=0; j<lgth2; j++ )
 328 {
 329         fprintf( stderr, "%5.2f ", currentw[j] );
 330 }
 331 fprintf( stderr, "\n" );
 332 #endif
 333 #if XXXXXXX
 334 fprintf( stderr, "\n" );
 335 fprintf( stderr, "i=%d\n", i );
 336 fprintf( stderr, "currentw = \n" );
 337 for( j=0; j<lgth2; j++ )
 338 {
 339         fprintf( stderr, "%5.2f ", currentw[j] );
 340 }
 341 fprintf( stderr, "\n" );
 342 #endif
 343                 currentw[0] = initverticalw[i];
 344
 345                 mi = previousw[0]; mpi = 0;
 346                 Mi = previousw[0]; Mpi = 0;
 347
 348                 ijpipt = ijpi[i] + 1;
 349                 ijpjpt = ijpj[i] + 1;
 350                 mjpt = m + 1;
 351                 Mjpt = largeM + 1;
 352                 prept = previousw;
 353                 curpt = currentw + 1;
 354                 mpjpt = mp + 1;
 355                 Mpjpt = Mp + 1;
 356                 tbk = -9999999.9;
 357                 tbki = 0;
 358                 tbkj = 0;
 359                 for( j=1; j<lgth2+1; j++ )
 360                 {
 361                         wm = *prept;
 362                         *ijpipt = i-1;
 363                         *ijpjpt = j-1;
 364
 365 #if 0
 366                         fprintf( stderr, "%5.0f->", wm );
 367 #endif
 368 #if 0
 369                         fprintf( stderr, "%5.0f?", g );
 370 #endif
 371                         if( (g=mi+fpenalty) > wm )
 372                         {
 373                                 wm = g;
 374 //                              *ijpipt = i - 1; // iranai
 375                                 *ijpjpt = mpi;
 376                         }
 377                         if( (g=*prept) >= mi )
 378                         {
 379                                 mi = g;
 380                                 mpi = j-1;
 381                         }
 382 #if USE_PENALTY_EX
 383                         mi += fpenalty_ex;
 384 #endif
 385
 386 #if 0
 387                         fprintf( stderr, "%5.0f?", g );
 388 #endif
 389                         if( (g=*mjpt + fpenalty) > wm )
 390                         {
 391                                 wm = g;
 392                                 *ijpipt = *mpjpt;
 393                                 *ijpjpt = j - 1; //IRU!
 394                         }
 395                         if( (g=*prept) >= *mjpt )
 396                         {
 397                                 *mjpt = g;
 398                                 *mpjpt = i-1;
 399                         }
 400 #if USE_PENALTY_EX
 401                         m[j] += fpenalty_ex;
 402 #endif
 403
 404 #if 1
 405                         g =  tbk + fpenalty_OP;
 406                         if( g > wm )
 407                         {
 408                                 wm = g;
 409                                 *ijpipt = tbki;
 410                                 *ijpjpt = tbkj;
 411 //                              fprintf( stderr, "hit! i%d, j%d, ijpi = %d, ijpj = %d\n", i, j, *ijpipt, *ijpjpt );
 412                         }
 413                         if( Mi > tbk )
 414                         {
 415                                 tbk = Mi; //error desu.
 416                                 tbki = i-1;
 417                                 tbkj = Mpi;
 418                         }
 419                         if( *Mjpt > tbk )
 420                         {
 421                                 tbk = *Mjpt;
 422                                 tbki = *Mpjpt;
 423                                 tbkj = j-1;
 424                         }
 425
 426                         if( *prept > *Mjpt )
 427                         {
 428                                 *Mjpt = *prept;
 429                                 *Mpjpt = i-1;
 430                         }
 431                         if( *prept > Mi )
 432                         {
 433                                 Mi = *prept;
 434                                 Mpi = j-1;
 435                         }
 436
 437 #endif
 438
 439
 440 #if 0
 441                         fprintf( stderr, "%5.0f ", wm );
 442 #endif
 443                         *curpt++ += wm;
 444                         ijpipt++;
 445                         ijpjpt++;
 446                         mjpt++;
 447                         Mjpt++;
 448                         prept++;
 449                         mpjpt++;
 450                         Mpjpt++;
 451                 }
 452                 lastverticalw[i] = currentw[lgth2-1]; // lgth2==0 no toki error
 453         }
 454 #if 0
 455         for( i=0; i<lgth1; i++ )
 456         {
 457                 for( j=0; j<lgth2; j++ )
 458                 {
 459                         fprintf( stdout, "i,j=%d,%d - ijpi,ijpj=%d,%d\n", i, j, ijpi[i][j], ijpj[i][j] );
 460                 }
 461         }
 462         fflush( stdout );
 463 #endif
 464
 465         genGtracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijpi, ijpj );
 466
 467         resultlen = strlen( mseq1[0] );
 468         if( alloclen < resultlen || resultlen > N )
 469         {
 470                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
 471                 ErrorExit( "LENGTH OVER!\n" );
 472         }
 473
 474
 475         strcpy( seq1[0], mseq1[0] );
 476         strcpy( seq2[0], mseq2[0] );
 477 #if 0
 478         fprintf( stderr, "\n" );
 479         fprintf( stderr, ">\n%s\n", mseq1[0] );
 480         fprintf( stderr, ">\n%s\n", mseq2[0] );
 481         fprintf( stderr, "wm = %f\n", wm );
 482 #endif
 483
 484         return( wm );
 485 }
 486