binaries/src/mafft/core/Salignmm.c

   1 #include "mltaln.h"
   2 #include "dp.h"
   3
   4 #define MACHIGAI 0
   5 #define OUTGAP0TRY 1
   6 #define DEBUG 0
   7 #define XXXXXXX    0
   8 #define USE_PENALTY_EX  0
   9 #define FASTMATCHCALC 1
  10
  11
  12 static float **impmtx = NULL;
  13 static int impalloclen = 0;
  14 float imp_match_out_sc( int i1, int j1 )
  15 {
  16 //      fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold );
  17 //      fprintf( stderr, "val = %f\n", impmtx[i1][j1] );
  18         return( impmtx[i1][j1] );
  19 }
  20 static void imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int *gapmap2 )
  21 {
  22 #if FASTMATCHCALC
  23         float *pt = impmtx[i1];
  24         int *gapmappt = gapmap2;
  25         while( lgth2-- )
  26                 *imp++ += pt[*gapmappt++];
  27 #else
  28         int j;
  29         float *pt = impmtx[i1];
  30         for( j=0; j<lgth2; j++ )
  31                 *imp++ += pt[gapmap2[j]];
  32 #endif
  33 }
  34
  35
  36 static void imp_match_out_vead_tate_gapmap( float *imp, int j1, int lgth1, int *gapmap1 )
  37 {
  38 #if FASTMATCHCALC
  39         int *gapmappt = gapmap1;
  40         while( lgth1-- )
  41                 *imp++ += impmtx[*gapmappt++][j1];
  42 #else
  43         int i;
  44         for( i=0; i<lgth1; i++ )
  45                 *imp++ += impmtx[gapmap1[i]][j1];
  46 #endif
  47 }
  48
  49 static void imp_match_out_vead( float *imp, int i1, int lgth2 )
  50 {
  51 #if FASTMATCHCALC
  52         float *pt = impmtx[i1];
  53         while( lgth2-- )
  54                 *imp++ += *pt++;
  55 #else
  56         int j;
  57         float *pt = impmtx[i1];
  58         for( j=0; j<lgth2; j++ )
  59                 *imp++ += pt[j];
  60 #endif
  61 }
  62 static void imp_match_out_vead_tate( float *imp, int j1, int lgth1 )
  63 {
  64         int i;
  65         for( i=0; i<lgth1; i++ )
  66                 *imp++ += impmtx[i][j1];
  67 }
  68
  69 void imp_rna( int nseq1, int nseq2, char **seq1, char **seq2, double *eff1, double *eff2, RNApair ***grouprna1, RNApair ***grouprna2, int *gapmap1, int *gapmap2, RNApair *pair )
  70 {
  71         foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, pair );
  72 }
  73
  74
  75 void imp_match_init_strict( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, double *eff1_kozo, double *eff2_kozo, LocalHom ***localhom, int forscore )
  76 {
  77         int i, j, k1, k2, tmpint, start1, start2, end1, end2;
  78         float effij;
  79         float effij_kozo;
  80         double effijx;
  81         char *pt, *pt1, *pt2;
  82         static char *nocount1 = NULL;
  83         static char *nocount2 = NULL;
  84         LocalHom *tmpptr;
  85
  86         if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
  87         {
  88                 if( impmtx ) FreeFloatMtx( impmtx );
  89                 if( nocount1 ) free( nocount1 );
  90                 if( nocount2 ) free( nocount2 );
  91                 impalloclen = MAX( lgth1, lgth2 ) + 2;
  92                 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
  93                 nocount1 = AllocateCharVec( impalloclen );
  94                 nocount2 = AllocateCharVec( impalloclen );
  95         }
  96
  97         for( i=0; i<lgth1; i++ )
  98         {
  99                 for( j=0; j<clus1; j++ )
 100                         if( seq1[j][i] == '-' ) break;
 101                 if( j != clus1 ) nocount1[i] = 1;
 102                 else                     nocount1[i] = 0;
 103         }
 104         for( i=0; i<lgth2; i++ )
 105         {
 106                 for( j=0; j<clus2; j++ )
 107                         if( seq2[j][i] == '-' ) break;
 108                 if( j != clus2 ) nocount2[i] = 1;
 109                 else                     nocount2[i] = 0;
 110         }
 111
 112 #if 0
 113 fprintf( stderr, "nocount2 =\n" );
 114 for( i = 0; i<impalloclen; i++ )
 115 {
 116         fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
 117 }
 118 #endif
 119
 120
 121
 122 #if 0
 123         fprintf( stderr, "eff1 in _init_strict = \n" );
 124         for( i=0; i<clus1; i++ )
 125                 fprintf( stderr, "eff1[] = %f\n", eff1[i] );
 126         for( i=0; i<clus2; i++ )
 127                 fprintf( stderr, "eff2[] = %f\n", eff2[i] );
 128 #endif
 129
 130         for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
 131                 impmtx[i][j] = 0.0;
 132         effijx =  fastathreshold;
 133         for( i=0; i<clus1; i++ )
 134         {
 135                 for( j=0; j<clus2; j++ )
 136                 {
 137                         effij = (float)( eff1[i] * eff2[j] * effijx );
 138                         effij_kozo = (float)( eff1_kozo[i] * eff2_kozo[j] * effijx );
 139                         tmpptr = localhom[i][j];
 140                         while( tmpptr )
 141                         {
 142 //                              fprintf( stderr, "start1 = %d\n", tmpptr->start1 );
 143 //                              fprintf( stderr, "end1   = %d\n", tmpptr->end1   );
 144 //                              fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] );
 145 //                              fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] );
 146                                 pt = seq1[i];
 147                                 tmpint = -1;
 148                                 while( *pt != 0 )
 149                                 {
 150                                         if( *pt++ != '-' ) tmpint++;
 151                                         if( tmpint == tmpptr->start1 ) break;
 152                                 }
 153                                 start1 = pt - seq1[i] - 1;
 154
 155                                 if( tmpptr->start1 == tmpptr->end1 ) end1 = start1;
 156                                 else
 157                                 {
 158 #if MACHIGAI
 159                                         while( *pt != 0 )
 160                                         {
 161 //                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
 162                                                 if( tmpint == tmpptr->end1 ) break;
 163                                                 if( *pt++ != '-' ) tmpint++;
 164                                         }
 165                                         end1 = pt - seq1[i] - 0;
 166 #else
 167                                         while( *pt != 0 )
 168                                         {
 169 //                                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] );
 170                                                 if( *pt++ != '-' ) tmpint++;
 171                                                 if( tmpint == tmpptr->end1 ) break;
 172                                         }
 173                                         end1 = pt - seq1[i] - 1;
 174 #endif
 175                                 }
 176
 177                                 pt = seq2[j];
 178                                 tmpint = -1;
 179                                 while( *pt != 0 )
 180                                 {
 181                                         if( *pt++ != '-' ) tmpint++;
 182                                         if( tmpint == tmpptr->start2 ) break;
 183                                 }
 184                                 start2 = pt - seq2[j] - 1;
 185                                 if( tmpptr->start2 == tmpptr->end2 ) end2 = start2;
 186                                 else
 187                                 {
 188 #if MACHIGAI
 189                                         while( *pt != 0 )
 190                                         {
 191                                                 if( tmpint == tmpptr->end2 ) break;
 192                                                 if( *pt++ != '-' ) tmpint++;
 193                                         }
 194                                         end2 = pt - seq2[j] - 0;
 195 #else
 196                                         while( *pt != 0 )
 197                                         {
 198                                                 if( *pt++ != '-' ) tmpint++;
 199                                                 if( tmpint == tmpptr->end2 ) break;
 200                                         }
 201                                         end2 = pt - seq2[j] - 1;
 202 #endif
 203                                 }
 204 //                              fprintf( stderr, "start1 = %d (%c), end1 = %d (%c), start2 = %d (%c), end2 = %d (%c)\n", start1, seq1[i][start1], end1, seq1[i][end1], start2, seq2[j][start2], end2, seq2[j][end2] );
 205 //                              fprintf( stderr, "step 0\n" );
 206                                 if( end1 - start1 != end2 - start2 )
 207                                 {
 208 //                                      fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
 209                                 }
 210
 211 #if 1
 212                                 k1 = start1; k2 = start2;
 213                                 pt1 = seq1[i] + k1;
 214                                 pt2 = seq2[j] + k2;
 215                                 while( *pt1 && *pt2 )
 216                                 {
 217                                         if( *pt1 != '-' && *pt2 != '-' )
 218                                         {
 219 // ½Å¤ß¤òÆó½Å¤Ë¤«¤±¤Ê¤¤¤è¤¦¤ËÃí°Õ¤·¤Æ²¼¤µ¤¤¡£
 220 //                                              impmtx[k1][k2] += tmpptr->wimportance * fastathreshold;
 221 //                                              impmtx[k1][k2] += tmpptr->importance * effij;
 222 //                                              impmtx[k1][k2] += tmpptr->fimportance * effij;
 223                                                 if( tmpptr->korh == 'k' )
 224                                                         impmtx[k1][k2] += tmpptr->fimportance * effij_kozo;
 225                                                 else
 226                                                         impmtx[k1][k2] += tmpptr->fimportance * effij;
 227
 228 //                                              fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij );
 229 //                                              fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 230 //                                              fprintf( stderr, "%d (%c) - %d (%c)  - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij );
 231                                                 k1++; k2++;
 232                                                 pt1++; pt2++;
 233                                         }
 234                                         else if( *pt1 != '-' && *pt2 == '-' )
 235                                         {
 236 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 237                                                 k2++; pt2++;
 238                                         }
 239                                         else if( *pt1 == '-' && *pt2 != '-' )
 240                                         {
 241 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 242                                                 k1++; pt1++;
 243                                         }
 244                                         else if( *pt1 == '-' && *pt2 == '-' )
 245                                         {
 246 //                                              fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 );
 247                                                 k1++; pt1++;
 248                                                 k2++; pt2++;
 249                                         }
 250                                         if( k1 > end1 || k2 > end2 ) break;
 251                                 }
 252 #else
 253                                 while( k1 <= end1 && k2 <= end2 )
 254                                 {
 255                                         fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 );
 256                                         if( !nocount1[k1] && !nocount2[k2] )
 257                                         {
 258                                                 impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j]  * fastathreshold;
 259                                                 fprintf( stderr, "marked\n" );
 260                                         }
 261                                         else
 262                                                 fprintf( stderr, "no count\n" );
 263                                         k1++; k2++;
 264                                 }
 265 #endif
 266                                 tmpptr = tmpptr->next;
 267                         }
 268                 }
 269         }
 270
 271 #if 0
 272         if( clus1 == 1 && clus2 == 1 )
 273         {
 274                 fprintf( stderr, "writing impmtx\n" );
 275                 fprintf( stderr, "\n" );
 276                 fprintf( stderr, "seq1[0] =  %s\n", seq1[0] );
 277                 fprintf( stderr, "seq2[0] =  %s\n", seq2[0] );
 278                 fprintf( stderr, "impmtx = \n" );
 279                 for( k2=0; k2<lgth2; k2++ )
 280                         fprintf( stderr, "%6.3f ", (double)k2 );
 281                 fprintf( stderr, "\n" );
 282                 for( k1=0; k1<lgth1; k1++ )
 283                 {
 284                         fprintf( stderr, "%d ", k1 );
 285                         for( k2=0; k2<30; k2++ )
 286                                 fprintf( stderr, "%2.1f ", impmtx[k1][k2] );
 287                         fprintf( stderr, "\n" );
 288                 }
 289 //              exit( 1 );
 290         }
 291 #endif
 292 }
 293
 294 #if 0
 295 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom )
 296 {
 297         int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2;
 298         static int impalloclen = 0;
 299         char *pt;
 300         int allgap;
 301         static char *nocount1 = NULL;
 302         static char *nocount2 = NULL;
 303
 304         if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 )
 305         {
 306                 if( impmtx ) FreeFloatMtx( impmtx );
 307                 if( nocount1 ) free( nocount1 );
 308                 if( nocount2 ) free( nocount2 );
 309                 impalloclen = MAX( lgth1, lgth2 ) + 2;
 310                 impmtx = AllocateFloatMtx( impalloclen, impalloclen );
 311                 nocount1 = AllocateCharVec( impalloclen );
 312                 nocount2 = AllocateCharVec( impalloclen );
 313         }
 314
 315         for( i=0; i<lgth1; i++ )
 316         {
 317                 for( j=0; j<clus1; j++ )
 318                         if( seq1[j][i] == '-' ) break;
 319                 if( j != clus1 ) nocount1[i] = 1;
 320                 else                     nocount1[i] = 0;
 321         }
 322         for( i=0; i<lgth2; i++ )
 323         {
 324                 for( j=0; j<clus2; j++ )
 325                         if( seq2[j][i] == '-' ) break;
 326                 if( j != clus2 ) nocount2[i] = 1;
 327                 else                     nocount2[i] = 0;
 328         }
 329
 330 #if 0
 331 fprintf( stderr, "nocount2 =\n" );
 332 for( i = 0; i<impalloclen; i++ )
 333 {
 334         fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] );
 335 }
 336 #endif
 337
 338         for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ )
 339                 impmtx[i][j] = 0;
 340         for( i=0; i<clus1; i++ )
 341         {
 342                 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] );
 343                 for( j=0; j<clus2; j++ )
 344                 {
 345                         fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 );
 346                         fprintf( stderr, "end1   = %d\n", localhom[i][j]->end1   );
 347                         fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] );
 348                         pt = seq1[i];
 349                         tmpint = -1;
 350                         while( *pt != 0 )
 351                         {
 352                                 if( *pt++ != '-' ) tmpint++;
 353                                 if( tmpint == localhom[i][j]->start1 ) break;
 354                         }
 355                         start1 = pt - seq1[i] - 1;
 356
 357                         while( *pt != 0 )
 358                         {
 359 //                              fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] );
 360                                 if( *pt++ != '-' ) tmpint++;
 361                                 if( tmpint == localhom[i][j]->end1 ) break;
 362                         }
 363                         end1 = pt - seq1[i] - 1;
 364
 365                         pt = seq2[j];
 366                         tmpint = -1;
 367                         while( *pt != 0 )
 368                         {
 369                                 if( *pt++ != '-' ) tmpint++;
 370                                 if( tmpint == localhom[i][j]->start2 ) break;
 371                         }
 372                         start2 = pt - seq2[j] - 1;
 373                         while( *pt != 0 )
 374                         {
 375                                 if( *pt++ != '-' ) tmpint++;
 376                                 if( tmpint == localhom[i][j]->end2 ) break;
 377                         }
 378                         end2 = pt - seq2[j] - 1;
 379 //                      fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 );
 380                         k1 = start1;
 381                         k2 = start2;
 382                         fprintf( stderr, "step 0\n" );
 383                         while( k1 <= end1 && k2 <= end2 )
 384                         {
 385 #if 0
 386                                 if( !nocount1[k1] && !nocount2[k2] )
 387                                         impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j];
 388                                 k1++; k2++;
 389 #else
 390                                 if( !nocount1[k1] && !nocount2[k2] )
 391                                         impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j];
 392                                 k1++; k2++;
 393 #endif
 394                         }
 395
 396                         dif = ( end1 - start1 ) - ( end2 - start2 );
 397                         fprintf( stderr, "dif = %d\n", dif );
 398                         if( dif > 0 )
 399                         {
 400                                 do
 401                                 {
 402                                         fprintf( stderr, "dif = %d\n", dif );
 403                                         k1 = start1;
 404                                         k2 = start2 - dif;
 405                                         while( k1 <= end1 && k2 <= end2 )
 406                                         {
 407                                                 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] )
 408                                                         impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
 409                                                 k1++; k2++;
 410                                         }
 411                                 }
 412                                 while( dif-- );
 413                         }
 414                         else
 415                         {
 416                                 do
 417                                 {
 418                                         k1 = start1 + dif;
 419                                         k2 = start2;
 420                                         while( k1 <= end1 )
 421                                         {
 422                                                 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] )
 423                                                         impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j];
 424                                                 k1++; k2++;
 425                                         }
 426                                 }
 427                                 while( dif++ );
 428                         }
 429                 }
 430         }
 431 #if 0
 432         fprintf( stderr, "impmtx = \n" );
 433         for( k2=0; k2<lgth2; k2++ )
 434                 fprintf( stderr, "%6.3f ", (double)k2 );
 435         fprintf( stderr, "\n" );
 436         for( k1=0; k1<lgth1; k1++ )
 437         {
 438                 fprintf( stderr, "%d", k1 );
 439                 for( k2=0; k2<lgth2; k2++ )
 440                         fprintf( stderr, "%6.3f ", impmtx[k1][k2] );
 441                 fprintf( stderr, "\n" );
 442         }
 443 #endif
 444 }
 445 #endif
 446
 447 static void match_calc( float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize )
 448 {
 449 #if FASTMATCHCALC
 450         int j, l;
 451         float scarr[26];
 452         float **cpmxpd = floatwork;
 453         int **cpmxpdn = intwork;
 454         float *matchpt, *cpmxpdpt, **cpmxpdptpt;
 455         int *cpmxpdnpt, **cpmxpdnptpt;
 456         if( initialize )
 457         {
 458                 int count = 0;
 459                 for( j=0; j<lgth2; j++ )
 460                 {
 461                         count = 0;
 462                         for( l=0; l<26; l++ )
 463                         {
 464                                 if( cpmx2[l][j] )
 465                                 {
 466                                         cpmxpd[j][count] = cpmx2[l][j];
 467                                         cpmxpdn[j][count] = l;
 468                                         count++;
 469                                 }
 470                         }
 471                         cpmxpdn[j][count] = -1;
 472                 }
 473         }
 474
 475         {
 476                 for( l=0; l<26; l++ )
 477                 {
 478                         scarr[l] = 0.0;
 479                         for( j=0; j<26; j++ )
 480 //                              scarr[l] += n_dis[j][l] * cpmx1[j][i1];
 481                                 scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1];
 482                 }
 483                 matchpt = match;
 484                 cpmxpdnptpt = cpmxpdn;
 485                 cpmxpdptpt = cpmxpd;
 486                 while( lgth2-- )
 487                 {
 488                         *matchpt = 0.0;
 489                         cpmxpdnpt = *cpmxpdnptpt++;
 490                         cpmxpdpt = *cpmxpdptpt++;
 491                         while( *cpmxpdnpt>-1 )
 492                                 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++;
 493                         matchpt++;
 494                 }
 495         }
 496 #else
 497         int j, k, l;
 498         float scarr[26];
 499         float **cpmxpd = floatwork;
 500         int **cpmxpdn = intwork;
 501 // simple
 502         if( initialize )
 503         {
 504                 int count = 0;
 505                 for( j=0; j<lgth2; j++ )
 506                 {
 507                         count = 0;
 508                         for( l=0; l<26; l++ )
 509                         {
 510                                 if( cpmx2[l][j] )
 511                                 {
 512                                         cpmxpd[count][j] = cpmx2[l][j];
 513                                         cpmxpdn[count][j] = l;
 514                                         count++;
 515                                 }
 516                         }
 517                         cpmxpdn[count][j] = -1;
 518                 }
 519         }
 520         for( l=0; l<26; l++ )
 521         {
 522                 scarr[l] = 0.0;
 523                 for( k=0; k<26; k++ )
 524                         scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1];
 525 //                      scarr[l] += n_dis[k][l] * cpmx1[k][i1];
 526         }
 527         for( j=0; j<lgth2; j++ )
 528         {
 529                 match[j] = 0.0;
 530                 for( k=0; cpmxpdn[k][j]>-1; k++ )
 531                         match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
 532         }
 533 #endif
 534 }
 535
 536 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
 537                                                 char **seq1, char **seq2,
 538                         char **mseq1, char **mseq2,
 539                         float **cpmx1, float **cpmx2,
 540                         int **ijp, int icyc, int jcyc )
 541 {
 542         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 543         float wm;
 544         char *gaptable1, *gt1bk;
 545         char *gaptable2, *gt2bk;
 546         lgth1 = strlen( seq1[0] );
 547         lgth2 = strlen( seq2[0] );
 548         gt1bk = AllocateCharVec( lgth1+lgth2+1 );
 549         gt2bk = AllocateCharVec( lgth1+lgth2+1 );
 550
 551 #if 0
 552         for( i=0; i<lgth1; i++ )
 553         {
 554                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 555         }
 556 #endif
 557
 558         if( outgap == 1 )
 559                 ;
 560         else
 561         {
 562                 wm = lastverticalw[0];
 563                 for( i=0; i<lgth1; i++ )
 564                 {
 565                         if( lastverticalw[i] >= wm )
 566                         {
 567                                 wm = lastverticalw[i];
 568                                 iin = i; jin = lgth2-1;
 569                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 570                         }
 571                 }
 572                 for( j=0; j<lgth2; j++ )
 573                 {
 574                         if( lasthorizontalw[j] >= wm )
 575                         {
 576                                 wm = lasthorizontalw[j];
 577                                 iin = lgth1-1; jin = j;
 578                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 579                         }
 580                 }
 581         }
 582
 583     for( i=0; i<lgth1+1; i++ )
 584     {
 585         ijp[i][0] = i + 1;
 586     }
 587     for( j=0; j<lgth2+1; j++ )
 588     {
 589         ijp[0][j] = -( j + 1 );
 590     }
 591
 592         gaptable1 = gt1bk + lgth1+lgth2;
 593         *gaptable1 = 0;
 594         gaptable2 = gt2bk + lgth1+lgth2;
 595         *gaptable2 = 0;
 596
 597         iin = lgth1; jin = lgth2;
 598         *impwmpt = 0.0;
 599         for( k=0; k<=lgth1+lgth2; k++ )
 600         {
 601                 if( ijp[iin][jin] < 0 )
 602                 {
 603                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 604                 }
 605                 else if( ijp[iin][jin] > 0 )
 606                 {
 607                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 608                 }
 609                 else
 610                 {
 611                         ifi = iin-1; jfi = jin-1;
 612                 }
 613                 l = iin - ifi;
 614                 while( --l )
 615                 {
 616                         *--gaptable1 = 'o';
 617                         *--gaptable2 = '-';
 618                         k++;
 619                 }
 620                 l= jin - jfi;
 621                 while( --l )
 622                 {
 623                         *--gaptable1 = '-';
 624                         *--gaptable2 = 'o';
 625                         k++;
 626                 }
 627                 if( iin == lgth1 || jin == lgth2 )
 628                         ;
 629                 else
 630                 {
 631                         *impwmpt += imp_match_out_sc( iin, jin );
 632
 633 //              fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 634                 }
 635                 if( iin <= 0 || jin <= 0 ) break;
 636                 *--gaptable1 = 'o';
 637                 *--gaptable2 = 'o';
 638                 k++;
 639                 iin = ifi; jin = jfi;
 640         }
 641
 642         for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
 643         for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
 644
 645         free( gt1bk );
 646         free( gt2bk );
 647 }
 648 static void Atracking_localhom_gapmap( float *impwmpt, float *lasthorizontalw, float *lastverticalw,
 649                                                 char **seq1, char **seq2,
 650                         char **mseq1, char **mseq2,
 651                         float **cpmx1, float **cpmx2,
 652                         int **ijp, int icyc, int jcyc,
 653                                                 int *gapmap1, int *gapmap2 )
 654 {
 655         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 656         float wm;
 657         char *gaptable1, *gt1bk;
 658         char *gaptable2, *gt2bk;
 659         lgth1 = strlen( seq1[0] );
 660         lgth2 = strlen( seq2[0] );
 661         gt1bk = AllocateCharVec( lgth1+lgth2+1 );
 662         gt2bk = AllocateCharVec( lgth1+lgth2+1 );
 663
 664 #if 0
 665         for( i=0; i<lgth1; i++ )
 666         {
 667                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 668         }
 669 #endif
 670
 671         if( outgap == 1 )
 672                 ;
 673         else
 674         {
 675                 wm = lastverticalw[0];
 676                 for( i=0; i<lgth1; i++ )
 677                 {
 678                         if( lastverticalw[i] >= wm )
 679                         {
 680                                 wm = lastverticalw[i];
 681                                 iin = i; jin = lgth2-1;
 682                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 683                         }
 684                 }
 685                 for( j=0; j<lgth2; j++ )
 686                 {
 687                         if( lasthorizontalw[j] >= wm )
 688                         {
 689                                 wm = lasthorizontalw[j];
 690                                 iin = lgth1-1; jin = j;
 691                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 692                         }
 693                 }
 694         }
 695
 696     for( i=0; i<lgth1+1; i++ )
 697     {
 698         ijp[i][0] = i + 1;
 699     }
 700     for( j=0; j<lgth2+1; j++ )
 701     {
 702         ijp[0][j] = -( j + 1 );
 703     }
 704
 705         gaptable1 = gt1bk + lgth1+lgth2;
 706         *gaptable1 = 0;
 707         gaptable2 = gt2bk + lgth1+lgth2;
 708         *gaptable2 = 0;
 709
 710         iin = lgth1; jin = lgth2;
 711         *impwmpt = 0.0;
 712         for( k=0; k<=lgth1+lgth2; k++ )
 713         {
 714                 if( ijp[iin][jin] < 0 )
 715                 {
 716                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 717                 }
 718                 else if( ijp[iin][jin] > 0 )
 719                 {
 720                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 721                 }
 722                 else
 723                 {
 724                         ifi = iin-1; jfi = jin-1;
 725                 }
 726                 l = iin - ifi;
 727                 while( --l )
 728                 {
 729                         *--gaptable1 = 'o';
 730                         *--gaptable2 = '-';
 731                         k++;
 732                 }
 733                 l= jin - jfi;
 734                 while( --l )
 735                 {
 736                         *--gaptable1 = '-';
 737                         *--gaptable2 = 'o';
 738                         k++;
 739                 }
 740                 if( iin == lgth1 || jin == lgth2 )
 741                         ;
 742                 else
 743                 {
 744                         *impwmpt += imp_match_out_sc( gapmap1[iin], gapmap2[jin] );
 745
 746 //              fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] );
 747                 }
 748                 if( iin <= 0 || jin <= 0 ) break;
 749                 *--gaptable1 = '-';
 750                 *--gaptable2 = '-';
 751                 k++;
 752                 iin = ifi; jin = jfi;
 753         }
 754         for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
 755         for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
 756
 757         free( gt1bk );
 758         free( gt2bk );
 759 }
 760 static float Atracking( float *lasthorizontalw, float *lastverticalw,
 761                                                 char **seq1, char **seq2,
 762                         char **mseq1, char **mseq2,
 763                         float **cpmx1, float **cpmx2,
 764                         int **ijp, int icyc, int jcyc )
 765 {
 766         int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k;
 767         float wm;
 768         char *gaptable1, *gt1bk;
 769         char *gaptable2, *gt2bk;
 770         lgth1 = strlen( seq1[0] );
 771         lgth2 = strlen( seq2[0] );
 772
 773         gt1bk = AllocateCharVec( lgth1+lgth2+1 );
 774         gt2bk = AllocateCharVec( lgth1+lgth2+1 );
 775
 776 #if 0
 777         for( i=0; i<lgth1; i++ )
 778         {
 779                 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
 780         }
 781 #endif
 782
 783         if( outgap == 1 )
 784                 ;
 785         else
 786         {
 787                 wm = lastverticalw[0];
 788                 for( i=0; i<lgth1; i++ )
 789                 {
 790                         if( lastverticalw[i] >= wm )
 791                         {
 792                                 wm = lastverticalw[i];
 793                                 iin = i; jin = lgth2-1;
 794                                 ijp[lgth1][lgth2] = +( lgth1 - i );
 795                         }
 796                 }
 797                 for( j=0; j<lgth2; j++ )
 798                 {
 799                         if( lasthorizontalw[j] >= wm )
 800                         {
 801                                 wm = lasthorizontalw[j];
 802                                 iin = lgth1-1; jin = j;
 803                                 ijp[lgth1][lgth2] = -( lgth2 - j );
 804                         }
 805                 }
 806         }
 807
 808     for( i=0; i<lgth1+1; i++ )
 809     {
 810         ijp[i][0] = i + 1;
 811     }
 812     for( j=0; j<lgth2+1; j++ )
 813     {
 814         ijp[0][j] = -( j + 1 );
 815     }
 816
 817         gaptable1 = gt1bk + lgth1+lgth2;
 818         *gaptable1 = 0;
 819         gaptable2 = gt2bk + lgth1+lgth2;
 820         *gaptable2 = 0;
 821
 822         iin = lgth1; jin = lgth2;
 823         for( k=0; k<=lgth1+lgth2; k++ )
 824         {
 825                 if( ijp[iin][jin] < 0 )
 826                 {
 827                         ifi = iin-1; jfi = jin+ijp[iin][jin];
 828                 }
 829                 else if( ijp[iin][jin] > 0 )
 830                 {
 831                         ifi = iin-ijp[iin][jin]; jfi = jin-1;
 832                 }
 833                 else
 834                 {
 835                         ifi = iin-1; jfi = jin-1;
 836                 }
 837                 l = iin - ifi;
 838                 while( --l )
 839                 {
 840                         *--gaptable1 = 'o';
 841                         *--gaptable2 = '-';
 842                         k++;
 843                 }
 844                 l= jin - jfi;
 845                 while( --l )
 846                 {
 847                         *--gaptable1 = '-';
 848                         *--gaptable2 = 'o';
 849                         k++;
 850                 }
 851                 if( iin <= 0 || jin <= 0 ) break;
 852                 *--gaptable1 = 'o';
 853                 *--gaptable2 = 'o';
 854                 k++;
 855                 iin = ifi; jin = jfi;
 856         }
 857
 858         for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 );
 859         for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 );
 860
 861         free( gt1bk );
 862         free( gt2bk );
 863
 864         return( 0.0 );
 865 }
 866
 867 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 )
 868 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
 869 {
 870 //      int k;
 871         register int i, j;
 872         int lasti, lastj;      /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
 873         int lgth1, lgth2;
 874         int resultlen;
 875         float wm = 0.0;   /* int ?????? */
 876         float g;
 877         float *currentw, *previousw;
 878 //      float fpenalty = (float)penalty;
 879 #if USE_PENALTY_EX
 880         float fpenalty_ex = (float)penalty_ex;
 881 #endif
 882 #if 1
 883         float *wtmp;
 884         int *ijppt;
 885         float *mjpt, *prept, *curpt;
 886         int *mpjpt;
 887 #endif
 888         static float mi, *m;
 889         static int **ijp;
 890         static int mpi, *mp;
 891         static float *w1, *w2;
 892         static float *match;
 893         static float *initverticalw;    /* kufuu sureba iranai */
 894         static float *lastverticalw;    /* kufuu sureba iranai */
 895         static char **mseq1;
 896         static char **mseq2;
 897         static char **mseq;
 898         static float *ogcp1;
 899         static float *ogcp2;
 900         static float *fgcp1;
 901         static float *fgcp2;
 902         static float **cpmx1;
 903         static float **cpmx2;
 904         static int **intwork;
 905         static float **floatwork;
 906         static int orlgth1 = 0, orlgth2 = 0;
 907         float fpenalty = (float)penalty;
 908         float *fgcp2pt;
 909         float *ogcp2pt;
 910         float fgcp1va;
 911         float ogcp1va;
 912
 913
 914
 915 #if 0
 916         fprintf( stderr, "####  eff in SA+++align\n" );
 917         fprintf( stderr, "####  seq1[0] = %s\n", seq1[0] );
 918         fprintf( stderr, "####  strlen( seq1[0] ) = %d\n", strlen( seq1[0] ) );
 919         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
 920 #endif
 921         if( orlgth1 == 0 )
 922         {
 923                 mseq1 = AllocateCharMtx( njob, 0 );
 924                 mseq2 = AllocateCharMtx( njob, 0 );
 925         }
 926
 927
 928         lgth1 = strlen( seq1[0] );
 929         lgth2 = strlen( seq2[0] );
 930 #if 0
 931         if( lgth1 == 0 || lgth2 == 0 )
 932         {
 933                 fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
 934         }
 935 #endif
 936
 937         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
 938         {
 939                 int ll1, ll2;
 940
 941                 if( orlgth1 > 0 && orlgth2 > 0 )
 942                 {
 943                         FreeFloatVec( w1 );
 944                         FreeFloatVec( w2 );
 945                         FreeFloatVec( match );
 946                         FreeFloatVec( initverticalw );
 947                         FreeFloatVec( lastverticalw );
 948
 949                         FreeFloatVec( m );
 950                         FreeIntVec( mp );
 951
 952                         FreeCharMtx( mseq );
 953
 954                         FreeFloatVec( ogcp1 );
 955                         FreeFloatVec( ogcp2 );
 956                         FreeFloatVec( fgcp1 );
 957                         FreeFloatVec( fgcp2 );
 958
 959
 960                         FreeFloatMtx( cpmx1 );
 961                         FreeFloatMtx( cpmx2 );
 962
 963                         FreeFloatMtx( floatwork );
 964                         FreeIntMtx( intwork );
 965                 }
 966
 967                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
 968                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
 969
 970 #if DEBUG
 971                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
 972 #endif
 973
 974                 w1 = AllocateFloatVec( ll2+2 );
 975                 w2 = AllocateFloatVec( ll2+2 );
 976                 match = AllocateFloatVec( ll2+2 );
 977
 978                 initverticalw = AllocateFloatVec( ll1+2 );
 979                 lastverticalw = AllocateFloatVec( ll1+2 );
 980
 981                 m = AllocateFloatVec( ll2+2 );
 982                 mp = AllocateIntVec( ll2+2 );
 983
 984                 mseq = AllocateCharMtx( njob, ll1+ll2 );
 985
 986                 ogcp1 = AllocateFloatVec( ll1+2 );
 987                 ogcp2 = AllocateFloatVec( ll2+2 );
 988                 fgcp1 = AllocateFloatVec( ll1+2 );
 989                 fgcp2 = AllocateFloatVec( ll2+2 );
 990
 991                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
 992                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
 993
 994 #if FASTMATCHCALC
 995                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
 996                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 27 );
 997 #else
 998                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
 999                 intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
1000 #endif
1001
1002 #if DEBUG
1003                 fprintf( stderr, "succeeded\n" );
1004 #endif
1005
1006                 orlgth1 = ll1 - 100;
1007                 orlgth2 = ll2 - 100;
1008         }
1009
1010
1011         for( i=0; i<icyc; i++ )
1012         {
1013                 mseq1[i] = mseq[i];
1014                 seq1[i][lgth1] = 0;
1015         }
1016         for( j=0; j<jcyc; j++ )
1017         {
1018                 mseq2[j] = mseq[icyc+j];
1019                 seq2[j][lgth2] = 0;
1020         }
1021
1022
1023         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1024         {
1025                 int ll1, ll2;
1026
1027                 if( commonAlloc1 && commonAlloc2 )
1028                 {
1029                         FreeIntMtx( commonIP );
1030                 }
1031
1032                 ll1 = MAX( orlgth1, commonAlloc1 );
1033                 ll2 = MAX( orlgth2, commonAlloc2 );
1034
1035 #if DEBUG
1036                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1037 #endif
1038
1039                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1040
1041 #if DEBUG
1042                 fprintf( stderr, "succeeded\n\n" );
1043 #endif
1044
1045                 commonAlloc1 = ll1;
1046                 commonAlloc2 = ll2;
1047         }
1048         ijp = commonIP;
1049
1050 #if 0
1051         {
1052                 float t = 0.0;
1053                 for( i=0; i<icyc; i++ )
1054                         t += eff1[i];
1055         fprintf( stderr, "## totaleff = %f\n", t );
1056         }
1057 #endif
1058
1059         cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc );
1060         cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc );
1061
1062         if( sgap1 )
1063         {
1064                 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 );
1065                 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 );
1066                 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 );
1067                 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 );
1068         }
1069         else
1070         {
1071                 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 );
1072                 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1073                 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1074                 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1075         }
1076
1077         for( i=0; i<lgth1; i++ )
1078         {
1079                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty;
1080                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty;
1081         }
1082         for( i=0; i<lgth2; i++ )
1083         {
1084                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty;
1085                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty;
1086         }
1087 #if 0
1088         for( i=0; i<lgth1; i++ )
1089                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1090 #endif
1091
1092         currentw = w1;
1093         previousw = w2;
1094
1095         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1096         if( localhom )
1097                 imp_match_out_vead_tate( initverticalw, 0, lgth1 ); // 060306
1098
1099         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1100         if( localhom )
1101                 imp_match_out_vead( currentw, 0, lgth2 ); // 060306
1102 #if 0 // -> tbfast.c
1103         if( localhom )
1104                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1105
1106 #endif
1107
1108         if( outgap == 1 )
1109         {
1110                 for( i=1; i<lgth1+1; i++ )
1111                 {
1112                         initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1113                 }
1114                 for( j=1; j<lgth2+1; j++ )
1115                 {
1116                         currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1117                 }
1118         }
1119 #if OUTGAP0TRY
1120         else
1121         {
1122                 for( j=1; j<lgth2+1; j++ )
1123                         currentw[j] -= offset * j / 2.0;
1124                 for( i=1; i<lgth1+1; i++ )
1125                         initverticalw[i] -= offset * i / 2.0;
1126         }
1127 #endif
1128
1129         for( j=1; j<lgth2+1; ++j )
1130         {
1131                 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1132         }
1133         if( lgth2 == 0 )
1134                 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari
1135         else
1136                 lastverticalw[0] = currentw[lgth2-1];
1137
1138         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1139
1140 #if XXXXXXX
1141 fprintf( stderr, "currentw = \n" );
1142 for( i=0; i<lgth1+1; i++ )
1143 {
1144         fprintf( stderr, "%5.2f ", currentw[i] );
1145 }
1146 fprintf( stderr, "\n" );
1147 fprintf( stderr, "initverticalw = \n" );
1148 for( i=0; i<lgth2+1; i++ )
1149 {
1150         fprintf( stderr, "%5.2f ", initverticalw[i] );
1151 }
1152 fprintf( stderr, "\n" );
1153 fprintf( stderr, "fcgp\n" );
1154 for( i=0; i<lgth1; i++ )
1155         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1156 for( i=0; i<lgth2; i++ )
1157         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1158 #endif
1159
1160         for( i=1; i<lasti; i++ )
1161         {
1162                 wtmp = previousw;
1163                 previousw = currentw;
1164                 currentw = wtmp;
1165
1166                 previousw[0] = initverticalw[i-1];
1167
1168                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1169 #if XXXXXXX
1170 fprintf( stderr, "\n" );
1171 fprintf( stderr, "i=%d\n", i );
1172 fprintf( stderr, "currentw = \n" );
1173 for( j=0; j<lgth2; j++ )
1174 {
1175         fprintf( stderr, "%5.2f ", currentw[j] );
1176 }
1177 fprintf( stderr, "\n" );
1178 #endif
1179                 if( localhom )
1180                 {
1181 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1182 #if  0
1183                         imp_match_out_vead( currentw, i, lgth2 );
1184 #else
1185                         imp_match_out_vead( currentw, i, lgth2 );
1186 #endif
1187                 }
1188 #if XXXXXXX
1189 fprintf( stderr, "\n" );
1190 fprintf( stderr, "i=%d\n", i );
1191 fprintf( stderr, "currentw = \n" );
1192 for( j=0; j<lgth2; j++ )
1193 {
1194         fprintf( stderr, "%5.2f ", currentw[j] );
1195 }
1196 fprintf( stderr, "\n" );
1197 #endif
1198                 currentw[0] = initverticalw[i];
1199
1200
1201                 mi = previousw[0] + ogcp2[1]; mpi = 0;
1202                 ijppt = ijp[i] + 1;
1203                 mjpt = m + 1;
1204                 prept = previousw;
1205                 curpt = currentw + 1;
1206                 mpjpt = mp + 1;
1207                 fgcp2pt = fgcp2;
1208                 ogcp2pt = ogcp2 + 1;
1209                 fgcp1va = fgcp1[i-1];
1210                 ogcp1va = ogcp1[i];
1211                 lastj = lgth2+1;
1212                 for( j=1; j<lastj; j++ )
1213                 {
1214                         wm = *prept;
1215                         *ijppt = 0;
1216
1217 #if 0
1218                         fprintf( stderr, "%5.0f->", wm );
1219 #endif
1220 #if 0
1221                         fprintf( stderr, "%5.0f?", g );
1222 #endif
1223                         if( (g=mi+*fgcp2pt) > wm )
1224                         {
1225                                 wm = g;
1226                                 *ijppt = -( j - mpi );
1227                         }
1228                         if( (g=*prept+*ogcp2pt) >= mi )
1229                         {
1230                                 mi = g;
1231                                 mpi = j-1;
1232                         }
1233 #if USE_PENALTY_EX
1234                         mi += fpenalty_ex;
1235 #endif
1236
1237 #if 0
1238                         fprintf( stderr, "%5.0f?", g );
1239 #endif
1240                         if( (g=*mjpt+fgcp1va) > wm )
1241                         {
1242                                 wm = g;
1243                                 *ijppt = +( i - *mpjpt );
1244                         }
1245                         if( (g=*prept+ogcp1va) >= *mjpt )
1246                         {
1247                                 *mjpt = g;
1248                                 *mpjpt = i-1;
1249                         }
1250 #if USE_PENALTY_EX
1251                         m[j] += fpenalty_ex;
1252 #endif
1253
1254 #if 0
1255                         fprintf( stderr, "%5.0f ", wm );
1256 #endif
1257                         *curpt++ += wm;
1258                         ijppt++;
1259                         mjpt++;
1260                         prept++;
1261                         mpjpt++;
1262                         fgcp2pt++;
1263                         ogcp2pt++;
1264                 }
1265                 lastverticalw[i] = currentw[lgth2-1];
1266         }
1267
1268 //      fprintf( stderr, "wm = %f\n", wm );
1269
1270 #if OUTGAP0TRY
1271         if( !outgap )
1272         {
1273                 for( j=1; j<lgth2+1; j++ )
1274                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1275                 for( i=1; i<lgth1+1; i++ )
1276                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1277         }
1278 #endif
1279
1280         /*
1281         fprintf( stderr, "\n" );
1282         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1283         fprintf( stderr, "#####\n" );
1284         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1285         fprintf( stderr, "====>" );
1286         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1287         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1288         */
1289         if( localhom )
1290         {
1291                 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1292         }
1293         else
1294                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1295
1296 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1297
1298         resultlen = strlen( mseq1[0] );
1299         if( alloclen < resultlen || resultlen > N )
1300         {
1301                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1302                 ErrorExit( "LENGTH OVER!\n" );
1303         }
1304
1305
1306         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1307         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1308         /*
1309         fprintf( stderr, "\n" );
1310         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1311         fprintf( stderr, "#####\n" );
1312         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1313         */
1314
1315 //      fprintf( stderr, "wm = %f\n", wm );
1316
1317         return( wm );
1318 }
1319
1320 float A__align_gapmap( char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, int *gapmap1, int *gapmap2 )
1321 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
1322 {
1323 //      int k;
1324         register int i, j;
1325         int lasti, lastj;     /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */
1326         int lgth1, lgth2;
1327         int resultlen;
1328         float wm = 0.0;   /* int ?????? */
1329         float g;
1330         float *currentw, *previousw;
1331 //      float fpenalty = (float)penalty;
1332 #if USE_PENALTY_EX
1333         float fpenalty_ex = (float)penalty_ex;
1334 #endif
1335 #if 1
1336         float *wtmp;
1337         int *ijppt;
1338         float *mjpt, *prept, *curpt;
1339         int *mpjpt;
1340 #endif
1341         static float mi, *m;
1342         static int **ijp;
1343         static int mpi, *mp;
1344         static float *w1, *w2;
1345         static float *match;
1346         static float *initverticalw;    /* kufuu sureba iranai */
1347         static float *lastverticalw;    /* kufuu sureba iranai */
1348         static char **mseq1;
1349         static char **mseq2;
1350         static char **mseq;
1351         static float *ogcp1;
1352         static float *ogcp2;
1353         static float *fgcp1;
1354         static float *fgcp2;
1355         static float **cpmx1;
1356         static float **cpmx2;
1357         static int **intwork;
1358         static float **floatwork;
1359         static int orlgth1 = 0, orlgth2 = 0;
1360         float *fgcp2pt;
1361         float *ogcp2pt;
1362         float fgcp1va;
1363         float ogcp1va;
1364
1365
1366 #if 0
1367         fprintf( stderr, "eff in SA+++align\n" );
1368         for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] );
1369 #endif
1370         if( orlgth1 == 0 )
1371         {
1372                 mseq1 = AllocateCharMtx( njob, 0 );
1373                 mseq2 = AllocateCharMtx( njob, 0 );
1374         }
1375
1376
1377         lgth1 = strlen( seq1[0] );
1378         lgth2 = strlen( seq2[0] );
1379 #if 0
1380         if( lgth1 == 0 || lgth2 == 0 )
1381         {
1382                 fprintf( stderr, "WARNING (Aalign_gapmap): lgth1=%d, lgth2=%d\n", lgth1, lgth2 );
1383         }
1384 #endif
1385
1386         if( lgth1 > orlgth1 || lgth2 > orlgth2 )
1387         {
1388                 int ll1, ll2;
1389
1390                 if( orlgth1 > 0 && orlgth2 > 0 )
1391                 {
1392                         FreeFloatVec( w1 );
1393                         FreeFloatVec( w2 );
1394                         FreeFloatVec( match );
1395                         FreeFloatVec( initverticalw );
1396                         FreeFloatVec( lastverticalw );
1397
1398                         FreeFloatVec( m );
1399                         FreeIntVec( mp );
1400
1401                         FreeCharMtx( mseq );
1402
1403                         FreeFloatVec( ogcp1 );
1404                         FreeFloatVec( ogcp2 );
1405                         FreeFloatVec( fgcp1 );
1406                         FreeFloatVec( fgcp2 );
1407
1408
1409                         FreeFloatMtx( cpmx1 );
1410                         FreeFloatMtx( cpmx2 );
1411
1412                         FreeFloatMtx( floatwork );
1413                         FreeIntMtx( intwork );
1414                 }
1415
1416                 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
1417                 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
1418
1419 #if DEBUG
1420                 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
1421 #endif
1422
1423                 w1 = AllocateFloatVec( ll2+2 );
1424                 w2 = AllocateFloatVec( ll2+2 );
1425                 match = AllocateFloatVec( ll2+2 );
1426
1427                 initverticalw = AllocateFloatVec( ll1+2 );
1428                 lastverticalw = AllocateFloatVec( ll1+2 );
1429
1430                 m = AllocateFloatVec( ll2+2 );
1431                 mp = AllocateIntVec( ll2+2 );
1432
1433                 mseq = AllocateCharMtx( njob, ll1+ll2 );
1434
1435                 ogcp1 = AllocateFloatVec( ll1+2 );
1436                 ogcp2 = AllocateFloatVec( ll2+2 );
1437                 fgcp1 = AllocateFloatVec( ll1+2 );
1438                 fgcp2 = AllocateFloatVec( ll2+2 );
1439
1440                 cpmx1 = AllocateFloatMtx( 26, ll1+2 );
1441                 cpmx2 = AllocateFloatMtx( 26, ll2+2 );
1442
1443 #if FASTMATCHCALC
1444                 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, 26 );
1445                 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, 26 );
1446 #else
1447                 floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
1448                 intwork = AllocateIntMtx( 27, MAX( ll1, ll2 )+2 );
1449 #endif
1450
1451 #if DEBUG
1452                 fprintf( stderr, "succeeded\n" );
1453 #endif
1454
1455                 orlgth1 = ll1 - 100;
1456                 orlgth2 = ll2 - 100;
1457         }
1458
1459
1460         for( i=0; i<icyc; i++ )
1461         {
1462                 mseq1[i] = mseq[i];
1463                 seq1[i][lgth1] = 0;
1464         }
1465         for( j=0; j<jcyc; j++ )
1466         {
1467                 mseq2[j] = mseq[icyc+j];
1468                 seq2[j][lgth2] = 0;
1469         }
1470
1471
1472         if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 )
1473         {
1474                 int ll1, ll2;
1475
1476                 if( commonAlloc1 && commonAlloc2 )
1477                 {
1478                         FreeIntMtx( commonIP );
1479                 }
1480
1481                 ll1 = MAX( orlgth1, commonAlloc1 );
1482                 ll2 = MAX( orlgth2, commonAlloc2 );
1483
1484 #if DEBUG
1485                 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
1486 #endif
1487
1488                 commonIP = AllocateIntMtx( ll1+10, ll2+10 );
1489
1490 #if DEBUG
1491                 fprintf( stderr, "succeeded\n\n" );
1492 #endif
1493
1494                 commonAlloc1 = ll1;
1495                 commonAlloc2 = ll2;
1496         }
1497         ijp = commonIP;
1498
1499         cpmx_calc_new( seq1, cpmx1, eff1, strlen( seq1[0] ), icyc );
1500         cpmx_calc_new( seq2, cpmx2, eff2, strlen( seq2[0] ), jcyc );
1501
1502         st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1  );
1503         st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 );
1504         st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 );
1505         st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 );
1506
1507         for( i=0; i<lgth1; i++ )
1508         {
1509                 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] );
1510                 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] );
1511         }
1512         for( i=0; i<lgth2; i++ )
1513         {
1514                 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] );
1515                 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] );
1516         }
1517 #if 0
1518         for( i=0; i<lgth1; i++ )
1519                 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
1520 #endif
1521
1522         currentw = w1;
1523         previousw = w2;
1524
1525
1526         match_calc( initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 );
1527         if( localhom )
1528                 imp_match_out_vead_tate_gapmap( initverticalw, gapmap2[0], lgth1, gapmap1 ); // 060306
1529
1530
1531         match_calc( currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 );
1532         if( localhom )
1533                 imp_match_out_vead_gapmap( currentw, gapmap1[0], lgth2, gapmap2 ); // 060306
1534 #if 0 // -> tbfast.c
1535         if( localhom )
1536                 imp_match_calc( currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 );
1537
1538 #endif
1539
1540         if( outgap == 1 )
1541         {
1542                 for( i=1; i<lgth1+1; i++ )
1543                 {
1544                         initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ;
1545                 }
1546                 for( j=1; j<lgth2+1; j++ )
1547                 {
1548                         currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ;
1549                 }
1550         }
1551 #if OUTGAP0TRY
1552         else
1553         {
1554                 for( j=1; j<lgth2+1; j++ )
1555                         currentw[j] -= offset * j / 2.0;
1556                 for( i=1; i<lgth1+1; i++ )
1557                         initverticalw[i] -= offset * i / 2.0;
1558         }
1559 #endif
1560
1561         for( j=1; j<lgth2+1; ++j )
1562         {
1563                 m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0;
1564         }
1565
1566         if( lgth2 == 0 )
1567                 lastverticalw[0] = 0.0;
1568         else
1569                 lastverticalw[0] = currentw[lgth2-1];
1570
1571         if( outgap ) lasti = lgth1+1; else lasti = lgth1;
1572
1573 #if XXXXXXX
1574 fprintf( stderr, "currentw = \n" );
1575 for( i=0; i<lgth1+1; i++ )
1576 {
1577         fprintf( stderr, "%5.2f ", currentw[i] );
1578 }
1579 fprintf( stderr, "\n" );
1580 fprintf( stderr, "initverticalw = \n" );
1581 for( i=0; i<lgth2+1; i++ )
1582 {
1583         fprintf( stderr, "%5.2f ", initverticalw[i] );
1584 }
1585 fprintf( stderr, "\n" );
1586 fprintf( stderr, "fcgp\n" );
1587 for( i=0; i<lgth1; i++ )
1588         fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] );
1589 for( i=0; i<lgth2; i++ )
1590         fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] );
1591 #endif
1592
1593         for( i=1; i<lasti; i++ )
1594         {
1595                 wtmp = previousw;
1596                 previousw = currentw;
1597                 currentw = wtmp;
1598
1599                 previousw[0] = initverticalw[i-1];
1600
1601                 match_calc( currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 );
1602 #if XXXXXXX
1603 fprintf( stderr, "\n" );
1604 fprintf( stderr, "i=%d\n", i );
1605 fprintf( stderr, "currentw = \n" );
1606 for( j=0; j<lgth2; j++ )
1607 {
1608         fprintf( stderr, "%5.2f ", currentw[j] );
1609 }
1610 fprintf( stderr, "\n" );
1611 #endif
1612                 if( localhom )
1613                 {
1614 //                      fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i );
1615 #if  0
1616                         imp_match_out_vead( currentw, i, lgth2 );
1617 #else
1618                         imp_match_out_vead_gapmap( currentw, gapmap1[i], lgth2, gapmap2 );
1619 #endif
1620                 }
1621 #if XXXXXXX
1622 fprintf( stderr, "\n" );
1623 fprintf( stderr, "i=%d\n", i );
1624 fprintf( stderr, "currentw = \n" );
1625 for( j=0; j<lgth2; j++ )
1626 {
1627         fprintf( stderr, "%5.2f ", currentw[j] );
1628 }
1629 fprintf( stderr, "\n" );
1630 #endif
1631                 currentw[0] = initverticalw[i];
1632
1633
1634                 mi = previousw[0] + ogcp2[1]; mpi = 0;
1635
1636                 ijppt = ijp[i] + 1;
1637                 mjpt = m + 1;
1638                 prept = previousw;
1639                 curpt = currentw + 1;
1640                 mpjpt = mp + 1;
1641                 fgcp2pt = fgcp2;
1642                 ogcp2pt = ogcp2 + 1;
1643                 fgcp1va = fgcp1[i-1];
1644                 ogcp1va = ogcp1[i];
1645                 lastj = lgth2+1;
1646                 for( j=1; j<lastj; j++ )
1647                 {
1648                         wm = *prept;
1649                         *ijppt = 0;
1650
1651 #if 0
1652                         fprintf( stderr, "%5.0f->", wm );
1653 #endif
1654                         g = mi + *fgcp2pt;
1655 #if 0
1656                         fprintf( stderr, "%5.0f?", g );
1657 #endif
1658                         if( g > wm )
1659                         {
1660                                 wm = g;
1661                                 *ijppt = -( j - mpi );
1662                         }
1663                         g = *prept + *ogcp2pt;
1664                         if( g >= mi )
1665                         {
1666                                 mi = g;
1667                                 mpi = j-1;
1668                         }
1669 #if USE_PENALTY_EX
1670                         mi += fpenalty_ex;
1671 #endif
1672
1673                         g = *mjpt + fgcp1va;
1674 #if 0
1675                         fprintf( stderr, "%5.0f?", g );
1676 #endif
1677                         if( g > wm )
1678                         {
1679                                 wm = g;
1680                                 *ijppt = +( i - *mpjpt );
1681                         }
1682                         g = *prept + ogcp1va;
1683                         if( g >= *mjpt )
1684                         {
1685                                 *mjpt = g;
1686                                 *mpjpt = i-1;
1687                         }
1688 #if USE_PENALTY_EX
1689                         m[j] += fpenalty_ex;
1690 #endif
1691
1692 #if 0
1693                         fprintf( stderr, "%5.0f ", wm );
1694 #endif
1695                         *curpt++ += wm;
1696                         ijppt++;
1697                         mjpt++;
1698                         prept++;
1699                         mpjpt++;
1700                         fgcp2pt++;
1701                         ogcp2pt++;
1702                 }
1703                 lastverticalw[i] = currentw[lgth2-1];
1704         }
1705
1706 #if OUTGAP0TRY
1707         if( !outgap )
1708         {
1709                 for( j=1; j<lgth2+1; j++ )
1710                         currentw[j] -= offset * ( lgth2 - j ) / 2.0;
1711                 for( i=1; i<lgth1+1; i++ )
1712                         lastverticalw[i] -= offset * ( lgth1 - i  / 2.0);
1713         }
1714 #endif
1715
1716         /*
1717         fprintf( stderr, "\n" );
1718         for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] );
1719         fprintf( stderr, "#####\n" );
1720         for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] );
1721         fprintf( stderr, "====>" );
1722         for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] );
1723         for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] );
1724         */
1725         if( localhom )
1726         {
1727                 Atracking_localhom_gapmap( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc, gapmap1, gapmap2 );
1728         }
1729         else
1730                 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijp, icyc, jcyc );
1731
1732 //      fprintf( stderr, "### impmatch = %f\n", *impmatch );
1733
1734         resultlen = strlen( mseq1[0] );
1735         if( alloclen < resultlen || resultlen > N )
1736         {
1737                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1738                 ErrorExit( "LENGTH OVER!\n" );
1739         }
1740
1741
1742         for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] );
1743         for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] );
1744         /*
1745         fprintf( stderr, "\n" );
1746         for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] );
1747         fprintf( stderr, "#####\n" );
1748         for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] );
1749         */
1750
1751 //      fprintf( stderr, "wm = %f\n", wm );
1752
1753         return( wm );
1754 }
1755
1756 float translate_and_Calign( char **mseq1, char **mseq2, double *effarr1, double *effarr2, int clus1, int clus2, int alloclen )
1757 {
1758     int i;
1759     float wm;
1760     char **result;
1761     char *seq = NULL, **aseq = NULL;
1762     double *effarr = NULL;
1763     int nseq = 0;
1764         int resultlen;
1765
1766     if     ( clus1 == 1 && clus2 != 1 )
1767     {
1768         seq = mseq1[0]; aseq = mseq2; effarr = effarr2; nseq = clus2+1;
1769 #if 0
1770                 printf( "effarr in transl... = \n" );
1771                 for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
1772 #endif
1773     }
1774     else if( clus1 != 1 && clus2 == 1 )
1775     {
1776         seq = mseq2[0]; aseq = mseq1; effarr = effarr1; nseq = clus1+1;
1777     }
1778     else ErrorExit( "ERROR in translate_and_Calign" );
1779
1780     result = Calignm1( &wm, aseq, seq, effarr, nseq-2, 0 );
1781
1782         resultlen = strlen( result[0] );
1783         if( alloclen < resultlen || resultlen > N )
1784         {
1785                 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
1786                 ErrorExit( "LENGTH OVER!\n" );
1787         }
1788     for( i=0; i<nseq-1; i++ ) strcpy( aseq[i], result[i] );
1789     strcpy( seq, result[nseq-1] );
1790
1791     return( 0.0 );
1792 }