binaries/src/mafft/core/Galign11.c.org

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define DEBUG 0
   5 #define XXXXXXX    0
   6 #define USE_PENALTY_EX  0
   7
   8
   9 #if 0
  10 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
  11 {
  12         int j;
  13         char tmpc = s1[0][i1];
  14         char *seq2 = s2[0];
  15
  16         while( lgth2-- )
  17                 *match++ = amino_dis[tmpc][*seq2++];
  18 }
  19 #else
  20 static void match_calc( float *match, char **s1, char **s2, int i1, int lgth2 )
  21 {
  22         int j;
  23
  24         for( j=0; j<lgth2; j++ )
  25                 match[j] = amino_dis[(*s1)[i1]][(*s2)[j]];
  26 }
  27 #endif
  28
  29 static float Atracking( float *lasthorizontalw, float *lastverticalw,
  30                                                 char **seq1, char **seq2,
  31                         char **mseq1, char **mseq2,
  32                         float **cpmx1, float **cpmx2,
  33                         short **ijp )
  34 {
  35         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk;
  36         char gap[] = "-";
  37         float wm;
  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         ijp[i][0] = i + 1;
  52     }
  53     for( j=0; j<lgth2+1; j++ )
  54     {
  55         ijp[0][j] = -( 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                 if( ijp[iin][jin] < 0 )
  68                 {
  69                         ifi = iin-1; jfi = jin+ijp[iin][jin];
  70                 }
  71                 else if( ijp[iin][jin] > 0 )
  72                 {
  73                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
  74                 }
  75                 else
  76                 {
  77                         ifi = iin-1; jfi = jin-1;
  78                 }
  79                 l = iin - ifi;
  80                 while( --l )
  81                 {
  82                         *--mseq1[0] = seq1[0][ifi+l];
  83                         *--mseq2[0] = *gap;
  84                         k++;
  85                 }
  86                 l= jin - jfi;
  87                 while( --l )
  88                 {
  89                         *--mseq1[0] = *gap;
  90                         *--mseq2[0] = seq2[0][jfi+l];
  91                         k++;
  92                 }
  93                 if( iin <= 0 || jin <= 0 ) break;
  94                 *--mseq1[0] = seq1[0][ifi];
  95                 *--mseq2[0] = seq2[0][jfi];
  96                 k++;
  97                 iin = ifi; jin = jfi;
  98         }
  99         return( 0.0 );
 100 }
 101
 102
 103 float G__align11( char **seq1, char **seq2, int alloclen )
 104 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 105 {
 106 //      int k;
 107         register int i, j;
 108         int lasti;                      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 109         int lgth1, lgth2;
 110         int resultlen;
 111         float wm;   /* int ?????? */
 112         float g;
 113         float *currentw, *previousw;
 114         float fpenalty = (float)penalty;
 115         float fpenalty_ex = (float)penalty_ex;
 116 #if 1
 117         float *wtmp;
 118         short *ijppt;
 119         float *mjpt, *prept, *curpt;
 120         int *mpjpt;
 121 #endif
 122         static float mi, *m;
 123         static short **ijp;
 124         static int mpi, *mp;
 125         static float *w1, *w2;
 126         static float *match;
 127         static float *initverticalw;    /* kufuu sureba iranai */
 128         static float *lastverticalw;    /* kufuu sureba iranai */
 129         static char **mseq1;
 130         static char **mseq2;
 131         static char **mseq;
 132         static float **cpmx1;
 133         static float **cpmx2;
 134         static int **intwork;
 135         static float **floatwork;
 136         static int orlgth1 = 0, orlgth2 = 0;
 137
 138         if( orlgth1 == 0 )
 139         {
 140                 mseq1 = AllocateCharMtx( njob, 0 );
 141                 mseq2 = AllocateCharMtx( njob, 0 );
 142         }
 143
 144
 145         lgth1 = strlen( seq1[0] );
 146         lgth2 = strlen( seq2[0] );
 147
 148         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 149         {
 150                 int ll1, ll2;
 151
 152                 if( orlgth1 > 0 && orlgth2 > 0 )
 153                 {
 154                         FreeFloatVec( w1 );
 155                         FreeFloatVec( w2 );
 156                         FreeFloatVec( match );
 157                         FreeFloatVec( initverticalw );
 158                         FreeFloatVec( lastverticalw );
 159
 160                         FreeFloatVec( m );
 161                         FreeIntVec( mp );
 162
 163                         FreeCharMtx( mseq );
 164
 165
 166                         FreeFloatMtx( cpmx1 );
 167                         FreeFloatMtx( cpmx2 );
 168
 169                         FreeFloatMtx( floatwork );
 170                         FreeIntMtx( intwork );
 171                 }
 172
 173                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 174                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 175
 176 #if DEBUG
 177                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 178 #endif
 179
 180                 w1 = AllocateFloatVec( ll2+2 );
 181                 w2 = AllocateFloatVec( ll2+2 );
 182                 match = AllocateFloatVec( ll2+2 );
 183
 184                 initverticalw = AllocateFloatVec( ll1+2 );
 185                 lastverticalw = AllocateFloatVec( ll1+2 );
 186
 187                 m = AllocateFloatVec( ll2+2 );
 188                 mp = AllocateIntVec( ll2+2 );
 189
 190                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 191
 192                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 193                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 194
 195                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 196                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
 197
 198 #if DEBUG
 199                 fprintf( stderr, "succeeded\n" );
 200 #endif
 201
 202                 orlgth1 = ll1 - 100;
 203                 orlgth2 = ll2 - 100;
 204         }
 205
 206
 207         mseq1[0] = mseq[0];
 208         mseq2[0] = mseq[1];
 209
 210
 211         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
 212         {
 213                 int ll1, ll2;
 214
 215                 if( commonAlloc1 && commonAlloc2 )
 216                 {
 217                         FreeShortMtx( commonIP );
 218                 }
 219
 220                 ll1 = MAX( orlgth1, commonAlloc1 );
 221                 ll2 = MAX( orlgth2, commonAlloc2 );
 222
 223 #if DEBUG
 224                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
 225 #endif
 226
 227                 commonIP = AllocateShortMtx( ll1+10, ll2+10 );
 228
 229 #if DEBUG
 230                 fprintf( stderr, "succeeded\n\n" );
 231 #endif
 232
 233                 commonAlloc1 = ll1;
 234                 commonAlloc2 = ll2;
 235         }
 236         ijp = commonIP;
 237
 238
 239 #if 0
 240         for( i=0; i<lgth1; i++ )
 241                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
 242 #endif
 243
 244         currentw = w1;
 245         previousw = w2;
 246
 247
 248         match_calc( initverticalw, seq2, seq1, 0, lgth1 );
 249
 250
 251         match_calc( currentw, seq1, seq2, 0, lgth2 );
 252
 253         if( outgap == 1 )
 254         {
 255                 for( i=1; i<lgth1+1; i++ )
 256                 {
 257                         initverticalw[i] += fpenalty;
 258                 }
 259                 for( j=1; j<lgth2+1; j++ )
 260                 {
 261                         currentw[j] += fpenalty;
 262                 }
 263         }
 264
 265         for( j=1; j<lgth2+1; ++j )
 266         {
 267                 m[j] = currentw[j-1]; mp[j] = 0;
 268         }
 269
 270         lastverticalw[0] = currentw[lgth2-1];
 271
 272         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
 273
 274 #if XXXXXXX
 275 fprintf( stderr, "currentw = \n" );
 276 for( i=0; i<lgth1+1; i++ )
 277 {
 278         fprintf( stderr, "%5.2f ", currentw[i] );
 279 }
 280 fprintf( stderr, "\n" );
 281 fprintf( stderr, "initverticalw = \n" );
 282 for( i=0; i<lgth2+1; i++ )
 283 {
 284         fprintf( stderr, "%5.2f ", initverticalw[i] );
 285 }
 286 fprintf( stderr, "\n" );
 287 #endif
 288
 289         for( i=1; i<lasti; i++ )
 290         {
 291                 wtmp = previousw;
 292                 previousw = currentw;
 293                 currentw = wtmp;
 294
 295                 previousw[0] = initverticalw[i-1];
 296
 297                 match_calc( currentw, seq1, seq2, i, lgth2 );
 298 #if XXXXXXX
 299 fprintf( stderr, "\n" );
 300 fprintf( stderr, "i=%d\n", i );
 301 fprintf( stderr, "currentw = \n" );
 302 for( j=0; j<lgth2; j++ )
 303 {
 304         fprintf( stderr, "%5.2f ", currentw[j] );
 305 }
 306 fprintf( stderr, "\n" );
 307 #endif
 308 #if XXXXXXX
 309 fprintf( stderr, "\n" );
 310 fprintf( stderr, "i=%d\n", i );
 311 fprintf( stderr, "currentw = \n" );
 312 for( j=0; j<lgth2; j++ )
 313 {
 314         fprintf( stderr, "%5.2f ", currentw[j] );
 315 }
 316 fprintf( stderr, "\n" );
 317 #endif
 318                 currentw[0] = initverticalw[i];
 319
 320                 mi = previousw[0]; mpi = 0;
 321
 322                 ijppt = ijp[i] + 1;
 323                 mjpt = m + 1;
 324                 prept = previousw;
 325                 curpt = currentw + 1;
 326                 mpjpt = mp + 1;
 327                 for( j=1; j<lgth2+1; j++ )
 328                 {
 329                         wm = *prept;
 330                         *ijppt = 0;
 331
 332 #if 0
 333                         fprintf( stderr, "%5.0f->", wm );
 334 #endif
 335                         g = mi + fpenalty;
 336 #if 0
 337                         fprintf( stderr, "%5.0f?", g );
 338 #endif
 339                         if( g > wm )
 340                         {
 341                                 wm = g;
 342                                 *ijppt = -( j - mpi );
 343                         }
 344                         g = *prept;
 345                         if( g >= mi )
 346                         {
 347                                 mi = g;
 348                                 mpi = j-1;
 349                         }
 350 #if USE_PENALTY_EX
 351                         mi += fpenalty_ex;
 352 #endif
 353
 354                         g = *mjpt + fpenalty;
 355 #if 0
 356                         fprintf( stderr, "%5.0f?", g );
 357 #endif
 358                         if( g > wm )
 359                         {
 360                                 wm = g;
 361                                 *ijppt = +( i - *mpjpt );
 362                         }
 363                         g = *prept;
 364                         if( g >= *mjpt )
 365                         {
 366                                 *mjpt = g;
 367                                 *mpjpt = i-1;
 368                         }
 369 #if USE_PENALTY_EX
 370                         m[j] += fpenalty_ex;
 371 #endif
 372
 373 #if 0
 374                         fprintf( stderr, "%5.0f ", wm );
 375 #endif
 376                         *curpt += wm;
 377                         ijppt++;
 378                         mjpt++;
 379                         prept++;
 380                         mpjpt++;
 381                         curpt++;
 382                 }
 383                 lastverticalw[i] = currentw[lgth2-1];
 384         }
 385
 386         Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp );
 387
 388
 389         resultlen = strlen( mseq1[0] );
 390         if( alloclen < resultlen || resultlen > N )
 391         {
 392                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
 393                 ErrorExit( "LENGTH OVER!\n" );
 394         }
 395
 396
 397         strcpy( seq1[0], mseq1[0] );
 398         strcpy( seq2[0], mseq2[0] );
 399 #if 0
 400         fprintf( stdout, "\n" );
 401         fprintf( stdout, ">\n%s\n", mseq1[0] );
 402         fprintf( stdout, ">\n%s\n", mseq2[0] );
 403 #endif
 404
 405         fprintf( stderr, "wm = %f\n", wm );
 406
 407
 408         return( wm );
 409 }
 410