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