Mercurial > repos > nick > duplex
comparison mafft/core/Salignmm.c @ 18:e4d75f9efb90 draft
planemo upload commit b'4303231da9e48b2719b4429a29b72421d24310f4\n'-dirty
| author | nick |
|---|---|
| date | Thu, 02 Feb 2017 18:44:31 -0500 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 17:836fa4fe9494 | 18:e4d75f9efb90 |
|---|---|
| 1 #include "mltaln.h" | |
| 2 #include "dp.h" | |
| 3 | |
| 4 #define MACHIGAI 0 | |
| 5 #define OUTGAP0TRY 0 | |
| 6 #define DEBUG 0 | |
| 7 #define XXXXXXX 0 | |
| 8 #define USE_PENALTY_EX 1 | |
| 9 #define FASTMATCHCALC 1 | |
| 10 #define MCD 0 | |
| 11 | |
| 12 | |
| 13 static TLS float **impmtx = NULL; | |
| 14 static TLS int impalloclen = 0; | |
| 15 float imp_match_out_sc( int i1, int j1 ) | |
| 16 { | |
| 17 // fprintf( stderr, "imp+match = %f\n", impmtx[i1][j1] * fastathreshold ); | |
| 18 // fprintf( stderr, "val = %f\n", impmtx[i1][j1] ); | |
| 19 return( impmtx[i1][j1] ); | |
| 20 } | |
| 21 | |
| 22 static void imp_match_out_vead_gapmap( float *imp, int i1, int lgth2, int *gapmap2 ) | |
| 23 { | |
| 24 #if FASTMATCHCALC | |
| 25 float *pt = impmtx[i1]; | |
| 26 int *gapmappt = gapmap2; | |
| 27 while( lgth2-- ) | |
| 28 *imp++ += pt[*gapmappt++]; | |
| 29 #else | |
| 30 int j; | |
| 31 float *pt = impmtx[i1]; | |
| 32 for( j=0; j<lgth2; j++ ) | |
| 33 *imp++ += pt[gapmap2[j]]; | |
| 34 #endif | |
| 35 } | |
| 36 | |
| 37 | |
| 38 static void imp_match_out_vead_tate_gapmap( float *imp, int j1, int lgth1, int *gapmap1 ) | |
| 39 { | |
| 40 #if FASTMATCHCALC | |
| 41 int *gapmappt = gapmap1; | |
| 42 while( lgth1-- ) | |
| 43 *imp++ += impmtx[*gapmappt++][j1]; | |
| 44 #else | |
| 45 int i; | |
| 46 for( i=0; i<lgth1; i++ ) | |
| 47 *imp++ += impmtx[gapmap1[i]][j1]; | |
| 48 #endif | |
| 49 } | |
| 50 | |
| 51 static void imp_match_out_vead( float *imp, int i1, int lgth2 ) | |
| 52 { | |
| 53 #if FASTMATCHCALC | |
| 54 float *pt = impmtx[i1]; | |
| 55 while( lgth2-- ) | |
| 56 *imp++ += *pt++; | |
| 57 #else | |
| 58 int j; | |
| 59 float *pt = impmtx[i1]; | |
| 60 for( j=0; j<lgth2; j++ ) | |
| 61 *imp++ += pt[j]; | |
| 62 #endif | |
| 63 } | |
| 64 static void imp_match_out_vead_tate( float *imp, int j1, int lgth1 ) | |
| 65 { | |
| 66 int i; | |
| 67 for( i=0; i<lgth1; i++ ) | |
| 68 *imp++ += impmtx[i][j1]; | |
| 69 } | |
| 70 | |
| 71 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 ) | |
| 72 { | |
| 73 foldrna( nseq1, nseq2, seq1, seq2, eff1, eff2, grouprna1, grouprna2, impmtx, gapmap1, gapmap2, pair ); | |
| 74 } | |
| 75 | |
| 76 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 ) | |
| 77 { | |
| 78 int i, j, k1, k2, tmpint, start1, start2, end1, end2; | |
| 79 float effij; | |
| 80 float effij_kozo; | |
| 81 double effijx; | |
| 82 char *pt, *pt1, *pt2; | |
| 83 static TLS char *nocount1 = NULL; | |
| 84 static TLS char *nocount2 = NULL; | |
| 85 LocalHom *tmpptr; | |
| 86 | |
| 87 if( seq1 == NULL ) | |
| 88 { | |
| 89 if( impmtx ) FreeFloatMtx( impmtx ); | |
| 90 impmtx = NULL; | |
| 91 if( nocount1 ) free( nocount1 ); | |
| 92 nocount1 = NULL; | |
| 93 if( nocount2 ) free( nocount2 ); | |
| 94 nocount2 = NULL; | |
| 95 | |
| 96 return; | |
| 97 } | |
| 98 | |
| 99 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 ) | |
| 100 { | |
| 101 if( impmtx ) FreeFloatMtx( impmtx ); | |
| 102 if( nocount1 ) free( nocount1 ); | |
| 103 if( nocount2 ) free( nocount2 ); | |
| 104 impalloclen = MAX( lgth1, lgth2 ) + 2; | |
| 105 impmtx = AllocateFloatMtx( impalloclen, impalloclen ); | |
| 106 nocount1 = AllocateCharVec( impalloclen ); | |
| 107 nocount2 = AllocateCharVec( impalloclen ); | |
| 108 } | |
| 109 | |
| 110 for( i=0; i<lgth1; i++ ) | |
| 111 { | |
| 112 for( j=0; j<clus1; j++ ) | |
| 113 if( seq1[j][i] == '-' ) break; | |
| 114 if( j != clus1 ) nocount1[i] = 1; | |
| 115 else nocount1[i] = 0; | |
| 116 } | |
| 117 for( i=0; i<lgth2; i++ ) | |
| 118 { | |
| 119 for( j=0; j<clus2; j++ ) | |
| 120 if( seq2[j][i] == '-' ) break; | |
| 121 if( j != clus2 ) nocount2[i] = 1; | |
| 122 else nocount2[i] = 0; | |
| 123 } | |
| 124 | |
| 125 #if 0 | |
| 126 fprintf( stderr, "nocount2 =\n" ); | |
| 127 for( i = 0; i<impalloclen; i++ ) | |
| 128 { | |
| 129 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] ); | |
| 130 } | |
| 131 #endif | |
| 132 | |
| 133 | |
| 134 | |
| 135 #if 0 | |
| 136 fprintf( stderr, "eff1 in _init_strict = \n" ); | |
| 137 for( i=0; i<clus1; i++ ) | |
| 138 fprintf( stderr, "eff1[] = %f\n", eff1[i] ); | |
| 139 for( i=0; i<clus2; i++ ) | |
| 140 fprintf( stderr, "eff2[] = %f\n", eff2[i] ); | |
| 141 #endif | |
| 142 | |
| 143 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ ) | |
| 144 impmtx[i][j] = 0.0; | |
| 145 effijx = fastathreshold; | |
| 146 for( i=0; i<clus1; i++ ) | |
| 147 { | |
| 148 for( j=0; j<clus2; j++ ) | |
| 149 { | |
| 150 effij = (float)( eff1[i] * eff2[j] * effijx ); | |
| 151 effij_kozo = (float)( eff1_kozo[i] * eff2_kozo[j] * effijx ); | |
| 152 tmpptr = localhom[i][j]; | |
| 153 while( tmpptr ) | |
| 154 { | |
| 155 // fprintf( stderr, "start1 = %d\n", tmpptr->start1 ); | |
| 156 // fprintf( stderr, "end1 = %d\n", tmpptr->end1 ); | |
| 157 // fprintf( stderr, "i = %d, seq1 = \n%s\n", i, seq1[i] ); | |
| 158 // fprintf( stderr, "j = %d, seq2 = \n%s\n", j, seq2[j] ); | |
| 159 pt = seq1[i]; | |
| 160 tmpint = -1; | |
| 161 while( *pt != 0 ) | |
| 162 { | |
| 163 if( *pt++ != '-' ) tmpint++; | |
| 164 if( tmpint == tmpptr->start1 ) break; | |
| 165 } | |
| 166 start1 = pt - seq1[i] - 1; | |
| 167 | |
| 168 if( tmpptr->start1 == tmpptr->end1 ) end1 = start1; | |
| 169 else | |
| 170 { | |
| 171 #if MACHIGAI | |
| 172 while( *pt != 0 ) | |
| 173 { | |
| 174 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] ); | |
| 175 if( tmpint == tmpptr->end1 ) break; | |
| 176 if( *pt++ != '-' ) tmpint++; | |
| 177 } | |
| 178 end1 = pt - seq1[i] - 0; | |
| 179 #else | |
| 180 while( *pt != 0 ) | |
| 181 { | |
| 182 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, tmpptr->end1, pt-seq1[i] ); | |
| 183 if( *pt++ != '-' ) tmpint++; | |
| 184 if( tmpint == tmpptr->end1 ) break; | |
| 185 } | |
| 186 end1 = pt - seq1[i] - 1; | |
| 187 #endif | |
| 188 } | |
| 189 | |
| 190 pt = seq2[j]; | |
| 191 tmpint = -1; | |
| 192 while( *pt != 0 ) | |
| 193 { | |
| 194 if( *pt++ != '-' ) tmpint++; | |
| 195 if( tmpint == tmpptr->start2 ) break; | |
| 196 } | |
| 197 start2 = pt - seq2[j] - 1; | |
| 198 if( tmpptr->start2 == tmpptr->end2 ) end2 = start2; | |
| 199 else | |
| 200 { | |
| 201 #if MACHIGAI | |
| 202 while( *pt != 0 ) | |
| 203 { | |
| 204 if( tmpint == tmpptr->end2 ) break; | |
| 205 if( *pt++ != '-' ) tmpint++; | |
| 206 } | |
| 207 end2 = pt - seq2[j] - 0; | |
| 208 #else | |
| 209 while( *pt != 0 ) | |
| 210 { | |
| 211 if( *pt++ != '-' ) tmpint++; | |
| 212 if( tmpint == tmpptr->end2 ) break; | |
| 213 } | |
| 214 end2 = pt - seq2[j] - 1; | |
| 215 #endif | |
| 216 } | |
| 217 // 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] ); | |
| 218 // fprintf( stderr, "step 0\n" ); | |
| 219 if( end1 - start1 != end2 - start2 ) | |
| 220 { | |
| 221 // fprintf( stderr, "CHUUI!!, start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 ); | |
| 222 } | |
| 223 | |
| 224 #if 1 | |
| 225 k1 = start1; k2 = start2; | |
| 226 pt1 = seq1[i] + k1; | |
| 227 pt2 = seq2[j] + k2; | |
| 228 while( *pt1 && *pt2 ) | |
| 229 { | |
| 230 if( *pt1 != '-' && *pt2 != '-' ) | |
| 231 { | |
| 232 // 重みを二重にかけないように注意して下さい。 | |
| 233 // impmtx[k1][k2] += tmpptr->wimportance * fastathreshold; | |
| 234 // impmtx[k1][k2] += tmpptr->importance * effij; | |
| 235 // impmtx[k1][k2] += tmpptr->fimportance * effij; | |
| 236 if( tmpptr->korh == 'k' ) | |
| 237 impmtx[k1][k2] += tmpptr->fimportance * effij_kozo; | |
| 238 else | |
| 239 impmtx[k1][k2] += tmpptr->fimportance * effij; | |
| 240 | |
| 241 // fprintf( stderr, "#### impmtx[k1][k2] = %f, tmpptr->fimportance=%f, effij=%f\n", impmtx[k1][k2], tmpptr->fimportance, effij ); | |
| 242 // fprintf( stderr, "mark, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); | |
| 243 // fprintf( stderr, "%d (%c) - %d (%c) - %f\n", k1, *pt1, k2, *pt2, tmpptr->fimportance * effij ); | |
| 244 k1++; k2++; | |
| 245 pt1++; pt2++; | |
| 246 } | |
| 247 else if( *pt1 != '-' && *pt2 == '-' ) | |
| 248 { | |
| 249 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); | |
| 250 k2++; pt2++; | |
| 251 } | |
| 252 else if( *pt1 == '-' && *pt2 != '-' ) | |
| 253 { | |
| 254 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); | |
| 255 k1++; pt1++; | |
| 256 } | |
| 257 else if( *pt1 == '-' && *pt2 == '-' ) | |
| 258 { | |
| 259 // fprintf( stderr, "skip, %d (%c) - %d (%c) \n", k1, *pt1, k2, *pt2 ); | |
| 260 k1++; pt1++; | |
| 261 k2++; pt2++; | |
| 262 } | |
| 263 if( k1 > end1 || k2 > end2 ) break; | |
| 264 } | |
| 265 #else | |
| 266 while( k1 <= end1 && k2 <= end2 ) | |
| 267 { | |
| 268 fprintf( stderr, "k1,k2=%d,%d - ", k1, k2 ); | |
| 269 if( !nocount1[k1] && !nocount2[k2] ) | |
| 270 { | |
| 271 impmtx[k1][k2] += tmpptr->wimportance * eff1[i] * eff2[j] * fastathreshold; | |
| 272 fprintf( stderr, "marked\n" ); | |
| 273 } | |
| 274 else | |
| 275 fprintf( stderr, "no count\n" ); | |
| 276 k1++; k2++; | |
| 277 } | |
| 278 #endif | |
| 279 tmpptr = tmpptr->next; | |
| 280 } | |
| 281 } | |
| 282 } | |
| 283 | |
| 284 #if 0 | |
| 285 if( clus1 == 1 && clus2 == 1 ) | |
| 286 { | |
| 287 fprintf( stderr, "writing impmtx\n" ); | |
| 288 fprintf( stderr, "\n" ); | |
| 289 fprintf( stderr, "seq1[0] = %s\n", seq1[0] ); | |
| 290 fprintf( stderr, "seq2[0] = %s\n", seq2[0] ); | |
| 291 fprintf( stderr, "impmtx = \n" ); | |
| 292 for( k2=0; k2<lgth2; k2++ ) | |
| 293 fprintf( stderr, "%6.3f ", (double)k2 ); | |
| 294 fprintf( stderr, "\n" ); | |
| 295 for( k1=0; k1<lgth1; k1++ ) | |
| 296 { | |
| 297 fprintf( stderr, "%d ", k1 ); | |
| 298 for( k2=0; k2<30; k2++ ) | |
| 299 fprintf( stderr, "%2.1f ", impmtx[k1][k2] ); | |
| 300 fprintf( stderr, "\n" ); | |
| 301 } | |
| 302 // exit( 1 ); | |
| 303 } | |
| 304 #endif | |
| 305 } | |
| 306 | |
| 307 #if 0 | |
| 308 void imp_match_init( float *imp, int clus1, int clus2, int lgth1, int lgth2, char **seq1, char **seq2, double *eff1, double *eff2, LocalHom ***localhom ) | |
| 309 { | |
| 310 int dif, i, j, k1, k2, tmpint, start1, start2, end1, end2; | |
| 311 static TLS int impalloclen = 0; | |
| 312 char *pt; | |
| 313 int allgap; | |
| 314 static TLS char *nocount1 = NULL; | |
| 315 static TLS char *nocount2 = NULL; | |
| 316 | |
| 317 if( impalloclen < lgth1 + 2 || impalloclen < lgth2 + 2 ) | |
| 318 { | |
| 319 if( impmtx ) FreeFloatMtx( impmtx ); | |
| 320 if( nocount1 ) free( nocount1 ); | |
| 321 if( nocount2 ) free( nocount2 ); | |
| 322 impalloclen = MAX( lgth1, lgth2 ) + 2; | |
| 323 impmtx = AllocateFloatMtx( impalloclen, impalloclen ); | |
| 324 nocount1 = AllocateCharVec( impalloclen ); | |
| 325 nocount2 = AllocateCharVec( impalloclen ); | |
| 326 } | |
| 327 | |
| 328 for( i=0; i<lgth1; i++ ) | |
| 329 { | |
| 330 for( j=0; j<clus1; j++ ) | |
| 331 if( seq1[j][i] == '-' ) break; | |
| 332 if( j != clus1 ) nocount1[i] = 1; | |
| 333 else nocount1[i] = 0; | |
| 334 } | |
| 335 for( i=0; i<lgth2; i++ ) | |
| 336 { | |
| 337 for( j=0; j<clus2; j++ ) | |
| 338 if( seq2[j][i] == '-' ) break; | |
| 339 if( j != clus2 ) nocount2[i] = 1; | |
| 340 else nocount2[i] = 0; | |
| 341 } | |
| 342 | |
| 343 #if 0 | |
| 344 fprintf( stderr, "nocount2 =\n" ); | |
| 345 for( i = 0; i<impalloclen; i++ ) | |
| 346 { | |
| 347 fprintf( stderr, "nocount2[%d] = %d (%c)\n", i, nocount2[i], seq2[0][i] ); | |
| 348 } | |
| 349 #endif | |
| 350 | |
| 351 for( i=0; i<lgth1; i++ ) for( j=0; j<lgth2; j++ ) | |
| 352 impmtx[i][j] = 0; | |
| 353 for( i=0; i<clus1; i++ ) | |
| 354 { | |
| 355 fprintf( stderr, "i = %d, seq1 = %s\n", i, seq1[i] ); | |
| 356 for( j=0; j<clus2; j++ ) | |
| 357 { | |
| 358 fprintf( stderr, "start1 = %d\n", localhom[i][j]->start1 ); | |
| 359 fprintf( stderr, "end1 = %d\n", localhom[i][j]->end1 ); | |
| 360 fprintf( stderr, "j = %d, seq2 = %s\n", j, seq2[j] ); | |
| 361 pt = seq1[i]; | |
| 362 tmpint = -1; | |
| 363 while( *pt != 0 ) | |
| 364 { | |
| 365 if( *pt++ != '-' ) tmpint++; | |
| 366 if( tmpint == localhom[i][j]->start1 ) break; | |
| 367 } | |
| 368 start1 = pt - seq1[i] - 1; | |
| 369 | |
| 370 while( *pt != 0 ) | |
| 371 { | |
| 372 // fprintf( stderr, "tmpint = %d, end1 = %d pos = %d\n", tmpint, localhom[i][j].end1, pt-seq1[i] ); | |
| 373 if( *pt++ != '-' ) tmpint++; | |
| 374 if( tmpint == localhom[i][j]->end1 ) break; | |
| 375 } | |
| 376 end1 = pt - seq1[i] - 1; | |
| 377 | |
| 378 pt = seq2[j]; | |
| 379 tmpint = -1; | |
| 380 while( *pt != 0 ) | |
| 381 { | |
| 382 if( *pt++ != '-' ) tmpint++; | |
| 383 if( tmpint == localhom[i][j]->start2 ) break; | |
| 384 } | |
| 385 start2 = pt - seq2[j] - 1; | |
| 386 while( *pt != 0 ) | |
| 387 { | |
| 388 if( *pt++ != '-' ) tmpint++; | |
| 389 if( tmpint == localhom[i][j]->end2 ) break; | |
| 390 } | |
| 391 end2 = pt - seq2[j] - 1; | |
| 392 // fprintf( stderr, "start1 = %d, end1 = %d, start2 = %d, end2 = %d\n", start1, end1, start2, end2 ); | |
| 393 k1 = start1; | |
| 394 k2 = start2; | |
| 395 fprintf( stderr, "step 0\n" ); | |
| 396 while( k1 <= end1 && k2 <= end2 ) | |
| 397 { | |
| 398 #if 0 | |
| 399 if( !nocount1[k1] && !nocount2[k2] ) | |
| 400 impmtx[k1][k2] += localhom[i][j].wimportance * eff1[i] * eff2[j]; | |
| 401 k1++; k2++; | |
| 402 #else | |
| 403 if( !nocount1[k1] && !nocount2[k2] ) | |
| 404 impmtx[k1][k2] += localhom[i][j]->wimportance * eff1[i] * eff2[j]; | |
| 405 k1++; k2++; | |
| 406 #endif | |
| 407 } | |
| 408 | |
| 409 dif = ( end1 - start1 ) - ( end2 - start2 ); | |
| 410 fprintf( stderr, "dif = %d\n", dif ); | |
| 411 if( dif > 0 ) | |
| 412 { | |
| 413 do | |
| 414 { | |
| 415 fprintf( stderr, "dif = %d\n", dif ); | |
| 416 k1 = start1; | |
| 417 k2 = start2 - dif; | |
| 418 while( k1 <= end1 && k2 <= end2 ) | |
| 419 { | |
| 420 if( 0 <= k2 && start2 <= k2 && !nocount1[k1] && !nocount2[k2] ) | |
| 421 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j]; | |
| 422 k1++; k2++; | |
| 423 } | |
| 424 } | |
| 425 while( dif-- ); | |
| 426 } | |
| 427 else | |
| 428 { | |
| 429 do | |
| 430 { | |
| 431 k1 = start1 + dif; | |
| 432 k2 = start2; | |
| 433 while( k1 <= end1 ) | |
| 434 { | |
| 435 if( k1 >= 0 && k1 >= start1 && !nocount1[k1] && !nocount2[k2] ) | |
| 436 impmtx[k1][k2] = localhom[i][j]->wimportance * eff1[i] * eff2[j]; | |
| 437 k1++; k2++; | |
| 438 } | |
| 439 } | |
| 440 while( dif++ ); | |
| 441 } | |
| 442 } | |
| 443 } | |
| 444 #if 0 | |
| 445 fprintf( stderr, "impmtx = \n" ); | |
| 446 for( k2=0; k2<lgth2; k2++ ) | |
| 447 fprintf( stderr, "%6.3f ", (double)k2 ); | |
| 448 fprintf( stderr, "\n" ); | |
| 449 for( k1=0; k1<lgth1; k1++ ) | |
| 450 { | |
| 451 fprintf( stderr, "%d", k1 ); | |
| 452 for( k2=0; k2<lgth2; k2++ ) | |
| 453 fprintf( stderr, "%6.3f ", impmtx[k1][k2] ); | |
| 454 fprintf( stderr, "\n" ); | |
| 455 } | |
| 456 #endif | |
| 457 } | |
| 458 #endif | |
| 459 | |
| 460 #if MCD | |
| 461 static void match_calc_dynamicmtx( double **pairoffset, double **n_dynamicmtx, float *match, int n1, char **seq1, double *eff1, int n2, char **seq2, double *eff2, int i1, int lgth2, float **floatwork, int **intwork, int initialize, int flip ) | |
| 462 { | |
| 463 // osoi! | |
| 464 int i, j, k; | |
| 465 int c1, c2; | |
| 466 double po; | |
| 467 // fprintf( stderr, "\nmatch_calc_dynamicmtx... %d", i1 ); | |
| 468 // fprintf( stderr, "\nseq1[0]=%s\n", seq1[0] ); | |
| 469 // fprintf( stderr, "\nseq2[0]=%s\n", seq2[0] ); | |
| 470 for( k=0; k<lgth2; k++ ) | |
| 471 { | |
| 472 match[k] = 0.0; | |
| 473 for( i=0; i<n1; i++ ) for( j=0; j<n2; j++ ) | |
| 474 { | |
| 475 if( flip ) po = pairoffset[j][i]; | |
| 476 else po = pairoffset[i][j]; | |
| 477 // if( k==0 ) fprintf( stderr, "pairoffset[%d][%d] = %f\n", i, j, po ); | |
| 478 c1 = amino_n[(int)seq1[i][i1]]; | |
| 479 c2 = amino_n[(int)seq2[j][k]]; | |
| 480 if( seq1[i][i1] == '-' || seq2[j][k] == '-' ) continue; | |
| 481 if( c1 < 0 || c2 < 0 ) continue; | |
| 482 // fprintf( stderr, "c1=%d, c2=%d\n", c1, c2 ); | |
| 483 if( flip ) | |
| 484 match[k] += ( n_dynamicmtx[c1][c2] + po * 600 ) * eff1[i] * eff2[j]; | |
| 485 else | |
| 486 match[k] += ( n_dynamicmtx[c1][c2] + po * 600 ) * eff1[i] * eff2[j]; | |
| 487 // fprintf( stderr, "match[k] = %f\n", match[k] ); | |
| 488 // fprintf( stderr, "pairoffset[i][j] = %f\n", pairoffset[i][j] ); | |
| 489 } | |
| 490 } | |
| 491 // fprintf( stderr, "done\n" ); | |
| 492 return; | |
| 493 } | |
| 494 #endif | |
| 495 | |
| 496 static void fillzero( float *s, int l ) | |
| 497 { | |
| 498 while( l-- ) *s++ = 0.0; | |
| 499 } | |
| 500 | |
| 501 static void match_calc_add( double **scoreingmtx, float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize ) | |
| 502 { | |
| 503 #if FASTMATCHCALC | |
| 504 // fprintf( stderr, "\nmatch_calc... %d", i1 ); | |
| 505 int j, l; | |
| 506 // float scarr[26]; | |
| 507 float **cpmxpd = floatwork; | |
| 508 int **cpmxpdn = intwork; | |
| 509 float *matchpt, *cpmxpdpt, **cpmxpdptpt; | |
| 510 int *cpmxpdnpt, **cpmxpdnptpt; | |
| 511 float *scarr; | |
| 512 scarr = calloc( nalphabets, sizeof( float ) ); | |
| 513 if( initialize ) | |
| 514 { | |
| 515 int count = 0; | |
| 516 for( j=0; j<lgth2; j++ ) | |
| 517 { | |
| 518 count = 0; | |
| 519 for( l=0; l<nalphabets; l++ ) | |
| 520 { | |
| 521 if( cpmx2[l][j] ) | |
| 522 { | |
| 523 cpmxpd[j][count] = cpmx2[l][j]; | |
| 524 cpmxpdn[j][count] = l; | |
| 525 count++; | |
| 526 } | |
| 527 } | |
| 528 cpmxpdn[j][count] = -1; | |
| 529 } | |
| 530 } | |
| 531 | |
| 532 { | |
| 533 for( l=0; l<nalphabets; l++ ) | |
| 534 { | |
| 535 scarr[l] = 0.0; | |
| 536 for( j=0; j<nalphabets; j++ ) | |
| 537 // scarr[l] += n_dis[j][l] * cpmx1[j][i1]; | |
| 538 // scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1]; | |
| 539 scarr[l] += scoreingmtx[j][l] * cpmx1[j][i1]; | |
| 540 } | |
| 541 matchpt = match; | |
| 542 cpmxpdnptpt = cpmxpdn; | |
| 543 cpmxpdptpt = cpmxpd; | |
| 544 while( lgth2-- ) | |
| 545 { | |
| 546 // *matchpt = 0.0; | |
| 547 cpmxpdnpt = *cpmxpdnptpt++; | |
| 548 cpmxpdpt = *cpmxpdptpt++; | |
| 549 while( *cpmxpdnpt>-1 ) | |
| 550 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++; | |
| 551 matchpt++; | |
| 552 } | |
| 553 } | |
| 554 free( scarr ); | |
| 555 // fprintf( stderr, "done\n" ); | |
| 556 #else | |
| 557 int j, k, l; | |
| 558 // float scarr[26]; | |
| 559 float **cpmxpd = floatwork; | |
| 560 int **cpmxpdn = intwork; | |
| 561 float *scarr; | |
| 562 scarr = calloc( nalphabets, sizeof( float ) ); | |
| 563 // simple | |
| 564 if( initialize ) | |
| 565 { | |
| 566 int count = 0; | |
| 567 for( j=0; j<lgth2; j++ ) | |
| 568 { | |
| 569 count = 0; | |
| 570 for( l=0; l<nalphabets; l++ ) | |
| 571 { | |
| 572 if( cpmx2[l][j] ) | |
| 573 { | |
| 574 cpmxpd[count][j] = cpmx2[l][j]; | |
| 575 cpmxpdn[count][j] = l; | |
| 576 count++; | |
| 577 } | |
| 578 } | |
| 579 cpmxpdn[count][j] = -1; | |
| 580 } | |
| 581 } | |
| 582 for( l=0; l<nalphabets; l++ ) | |
| 583 { | |
| 584 scarr[l] = 0.0; | |
| 585 for( k=0; k<nalphabets; k++ ) | |
| 586 // scarr[l] += n_dis[k][l] * cpmx1[k][i1]; | |
| 587 // scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1]; | |
| 588 scarr[l] += scoreingmtx[k][l] * cpmx1[k][i1]; | |
| 589 } | |
| 590 for( j=0; j<lgth2; j++ ) | |
| 591 { | |
| 592 match[j] = 0.0; | |
| 593 for( k=0; cpmxpdn[k][j]>-1; k++ ) | |
| 594 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j]; | |
| 595 } | |
| 596 free( scarr ); | |
| 597 #endif | |
| 598 } | |
| 599 static void match_calc( double **n_dynamicmtx, float *match, float **cpmx1, float **cpmx2, int i1, int lgth2, float **floatwork, int **intwork, int initialize ) | |
| 600 { | |
| 601 #if FASTMATCHCALC | |
| 602 // fprintf( stderr, "\nmatch_calc... %d", i1 ); | |
| 603 int j, l; | |
| 604 // float scarr[26]; | |
| 605 float **cpmxpd = floatwork; | |
| 606 int **cpmxpdn = intwork; | |
| 607 float *matchpt, *cpmxpdpt, **cpmxpdptpt; | |
| 608 int *cpmxpdnpt, **cpmxpdnptpt; | |
| 609 float *scarr; | |
| 610 scarr = calloc( nalphabets, sizeof( float ) ); | |
| 611 if( initialize ) | |
| 612 { | |
| 613 int count = 0; | |
| 614 for( j=0; j<lgth2; j++ ) | |
| 615 { | |
| 616 count = 0; | |
| 617 for( l=0; l<nalphabets; l++ ) | |
| 618 { | |
| 619 if( cpmx2[l][j] ) | |
| 620 { | |
| 621 cpmxpd[j][count] = cpmx2[l][j]; | |
| 622 cpmxpdn[j][count] = l; | |
| 623 count++; | |
| 624 } | |
| 625 } | |
| 626 cpmxpdn[j][count] = -1; | |
| 627 } | |
| 628 } | |
| 629 | |
| 630 { | |
| 631 for( l=0; l<nalphabets; l++ ) | |
| 632 { | |
| 633 scarr[l] = 0.0; | |
| 634 for( j=0; j<nalphabets; j++ ) | |
| 635 // scarr[l] += n_dis[j][l] * cpmx1[j][i1]; | |
| 636 // scarr[l] += n_dis_consweight_multi[j][l] * cpmx1[j][i1]; | |
| 637 scarr[l] += n_dynamicmtx[j][l] * cpmx1[j][i1]; | |
| 638 } | |
| 639 matchpt = match; | |
| 640 cpmxpdnptpt = cpmxpdn; | |
| 641 cpmxpdptpt = cpmxpd; | |
| 642 while( lgth2-- ) | |
| 643 { | |
| 644 *matchpt = 0.0; | |
| 645 cpmxpdnpt = *cpmxpdnptpt++; | |
| 646 cpmxpdpt = *cpmxpdptpt++; | |
| 647 while( *cpmxpdnpt>-1 ) | |
| 648 *matchpt += scarr[*cpmxpdnpt++] * *cpmxpdpt++; | |
| 649 matchpt++; | |
| 650 } | |
| 651 } | |
| 652 free( scarr ); | |
| 653 // fprintf( stderr, "done\n" ); | |
| 654 #else | |
| 655 int j, k, l; | |
| 656 // float scarr[26]; | |
| 657 float **cpmxpd = floatwork; | |
| 658 int **cpmxpdn = intwork; | |
| 659 float *scarr; | |
| 660 scarr = calloc( nalphabets, sizeof( float ) ); | |
| 661 // simple | |
| 662 if( initialize ) | |
| 663 { | |
| 664 int count = 0; | |
| 665 for( j=0; j<lgth2; j++ ) | |
| 666 { | |
| 667 count = 0; | |
| 668 for( l=0; l<nalphabets; l++ ) | |
| 669 { | |
| 670 if( cpmx2[l][j] ) | |
| 671 { | |
| 672 cpmxpd[count][j] = cpmx2[l][j]; | |
| 673 cpmxpdn[count][j] = l; | |
| 674 count++; | |
| 675 } | |
| 676 } | |
| 677 cpmxpdn[count][j] = -1; | |
| 678 } | |
| 679 } | |
| 680 for( l=0; l<nalphabets; l++ ) | |
| 681 { | |
| 682 scarr[l] = 0.0; | |
| 683 for( k=0; k<nalphabets; k++ ) | |
| 684 // scarr[l] += n_dis[k][l] * cpmx1[k][i1]; | |
| 685 // scarr[l] += n_dis_consweight_multi[k][l] * cpmx1[k][i1]; | |
| 686 scarr[l] += n_dynamicmtx[k][l] * cpmx1[k][i1]; | |
| 687 } | |
| 688 for( j=0; j<lgth2; j++ ) | |
| 689 { | |
| 690 match[j] = 0.0; | |
| 691 for( k=0; cpmxpdn[k][j]>-1; k++ ) | |
| 692 match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j]; | |
| 693 } | |
| 694 free( scarr ); | |
| 695 #endif | |
| 696 } | |
| 697 | |
| 698 static void Atracking_localhom( float *impwmpt, float *lasthorizontalw, float *lastverticalw, | |
| 699 char **seq1, char **seq2, | |
| 700 char **mseq1, char **mseq2, | |
| 701 int **ijp, int icyc, int jcyc, | |
| 702 int *warpis, int *warpjs, int warpbase ) | |
| 703 { | |
| 704 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk; | |
| 705 float wm; | |
| 706 char *gaptable1, *gt1bk; | |
| 707 char *gaptable2, *gt2bk; | |
| 708 lgth1 = strlen( seq1[0] ); | |
| 709 lgth2 = strlen( seq2[0] ); | |
| 710 gt1bk = AllocateCharVec( lgth1+lgth2+1 ); | |
| 711 gt2bk = AllocateCharVec( lgth1+lgth2+1 ); | |
| 712 | |
| 713 #if 0 | |
| 714 for( i=0; i<lgth1; i++ ) | |
| 715 { | |
| 716 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] ); | |
| 717 } | |
| 718 #endif | |
| 719 | |
| 720 if( outgap == 1 ) | |
| 721 ; | |
| 722 else | |
| 723 { | |
| 724 wm = lastverticalw[0]; | |
| 725 for( i=0; i<lgth1; i++ ) | |
| 726 { | |
| 727 if( lastverticalw[i] >= wm ) | |
| 728 { | |
| 729 wm = lastverticalw[i]; | |
| 730 iin = i; jin = lgth2-1; | |
| 731 ijp[lgth1][lgth2] = +( lgth1 - i ); | |
| 732 } | |
| 733 } | |
| 734 for( j=0; j<lgth2; j++ ) | |
| 735 { | |
| 736 if( lasthorizontalw[j] >= wm ) | |
| 737 { | |
| 738 wm = lasthorizontalw[j]; | |
| 739 iin = lgth1-1; jin = j; | |
| 740 ijp[lgth1][lgth2] = -( lgth2 - j ); | |
| 741 } | |
| 742 } | |
| 743 } | |
| 744 | |
| 745 for( i=0; i<lgth1+1; i++ ) | |
| 746 { | |
| 747 ijp[i][0] = i + 1; | |
| 748 } | |
| 749 for( j=0; j<lgth2+1; j++ ) | |
| 750 { | |
| 751 ijp[0][j] = -( j + 1 ); | |
| 752 } | |
| 753 | |
| 754 gaptable1 = gt1bk + lgth1+lgth2; | |
| 755 *gaptable1 = 0; | |
| 756 gaptable2 = gt2bk + lgth1+lgth2; | |
| 757 *gaptable2 = 0; | |
| 758 | |
| 759 iin = lgth1; jin = lgth2; | |
| 760 limk = lgth1+lgth2 + 1; | |
| 761 *impwmpt = 0.0; | |
| 762 for( k=0; k<limk; k++ ) | |
| 763 { | |
| 764 if( ijp[iin][jin] >= warpbase ) | |
| 765 { | |
| 766 ifi = warpis[ijp[iin][jin]-warpbase]; | |
| 767 jfi = warpjs[ijp[iin][jin]-warpbase]; | |
| 768 } | |
| 769 else if( ijp[iin][jin] < 0 ) | |
| 770 { | |
| 771 ifi = iin-1; jfi = jin+ijp[iin][jin]; | |
| 772 } | |
| 773 else if( ijp[iin][jin] > 0 ) | |
| 774 { | |
| 775 ifi = iin-ijp[iin][jin]; jfi = jin-1; | |
| 776 } | |
| 777 else | |
| 778 { | |
| 779 ifi = iin-1; jfi = jin-1; | |
| 780 } | |
| 781 if( ifi == -warpbase && jfi == -warpbase ) | |
| 782 { | |
| 783 l = iin; | |
| 784 while( --l >= 0 ) | |
| 785 { | |
| 786 *--gaptable1 = 'o'; | |
| 787 *--gaptable2 = '-'; | |
| 788 k++; | |
| 789 } | |
| 790 l= jin; | |
| 791 while( --l >= 0 ) | |
| 792 { | |
| 793 *--gaptable1 = '-'; | |
| 794 *--gaptable2 = 'o'; | |
| 795 } | |
| 796 break; | |
| 797 } | |
| 798 else | |
| 799 { | |
| 800 l = iin - ifi; | |
| 801 while( --l ) | |
| 802 { | |
| 803 *--gaptable1 = 'o'; | |
| 804 *--gaptable2 = '-'; | |
| 805 k++; | |
| 806 } | |
| 807 l= jin - jfi; | |
| 808 while( --l ) | |
| 809 { | |
| 810 *--gaptable1 = '-'; | |
| 811 *--gaptable2 = 'o'; | |
| 812 k++; | |
| 813 } | |
| 814 } | |
| 815 if( iin == lgth1 || jin == lgth2 ) | |
| 816 ; | |
| 817 else | |
| 818 { | |
| 819 *impwmpt += imp_match_out_sc( iin, jin ); | |
| 820 | |
| 821 // fprintf( stderr, "impwm = %f (iin=%d, jin=%d) seq1=%c, seq2=%c\n", *impwmpt, iin, jin, seq1[0][iin], seq2[0][jin] ); | |
| 822 } | |
| 823 if( iin <= 0 || jin <= 0 ) break; | |
| 824 *--gaptable1 = 'o'; | |
| 825 *--gaptable2 = 'o'; | |
| 826 k++; | |
| 827 iin = ifi; jin = jfi; | |
| 828 } | |
| 829 | |
| 830 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 ); | |
| 831 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 ); | |
| 832 | |
| 833 free( gt1bk ); | |
| 834 free( gt2bk ); | |
| 835 } | |
| 836 | |
| 837 static float Atracking( float *lasthorizontalw, float *lastverticalw, | |
| 838 char **seq1, char **seq2, | |
| 839 char **mseq1, char **mseq2, | |
| 840 int **ijp, int icyc, int jcyc, | |
| 841 int tailgp, | |
| 842 int *warpis, int *warpjs, int warpbase ) | |
| 843 { | |
| 844 int i, j, l, iin, jin, ifi, jfi, lgth1, lgth2, k, limk; | |
| 845 float wm; | |
| 846 char *gaptable1, *gt1bk; | |
| 847 char *gaptable2, *gt2bk; | |
| 848 lgth1 = strlen( seq1[0] ); | |
| 849 lgth2 = strlen( seq2[0] ); | |
| 850 | |
| 851 gt1bk = AllocateCharVec( lgth1+lgth2+1 ); | |
| 852 gt2bk = AllocateCharVec( lgth1+lgth2+1 ); | |
| 853 | |
| 854 #if 0 | |
| 855 for( i=0; i<lgth1; i++ ) | |
| 856 { | |
| 857 fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] ); | |
| 858 } | |
| 859 #endif | |
| 860 | |
| 861 if( tailgp == 1 ) | |
| 862 ; | |
| 863 else | |
| 864 { | |
| 865 wm = lastverticalw[0]; | |
| 866 for( i=0; i<lgth1; i++ ) | |
| 867 { | |
| 868 if( lastverticalw[i] >= wm ) | |
| 869 { | |
| 870 wm = lastverticalw[i]; | |
| 871 iin = i; jin = lgth2-1; | |
| 872 ijp[lgth1][lgth2] = +( lgth1 - i ); | |
| 873 } | |
| 874 } | |
| 875 for( j=0; j<lgth2; j++ ) | |
| 876 { | |
| 877 if( lasthorizontalw[j] >= wm ) | |
| 878 { | |
| 879 wm = lasthorizontalw[j]; | |
| 880 iin = lgth1-1; jin = j; | |
| 881 ijp[lgth1][lgth2] = -( lgth2 - j ); | |
| 882 } | |
| 883 } | |
| 884 } | |
| 885 | |
| 886 for( i=0; i<lgth1+1; i++ ) | |
| 887 { | |
| 888 ijp[i][0] = i + 1; | |
| 889 } | |
| 890 for( j=0; j<lgth2+1; j++ ) | |
| 891 { | |
| 892 ijp[0][j] = -( j + 1 ); | |
| 893 } | |
| 894 | |
| 895 gaptable1 = gt1bk + lgth1+lgth2; | |
| 896 *gaptable1 = 0; | |
| 897 gaptable2 = gt2bk + lgth1+lgth2; | |
| 898 *gaptable2 = 0; | |
| 899 | |
| 900 iin = lgth1; jin = lgth2; | |
| 901 limk = lgth1+lgth2 + 1; | |
| 902 for( k=0; k<limk; k++ ) | |
| 903 { | |
| 904 if( ijp[iin][jin] >= warpbase ) | |
| 905 { | |
| 906 ifi = warpis[ijp[iin][jin]-warpbase]; | |
| 907 jfi = warpjs[ijp[iin][jin]-warpbase]; | |
| 908 } | |
| 909 else if( ijp[iin][jin] < 0 ) | |
| 910 { | |
| 911 ifi = iin-1; jfi = jin+ijp[iin][jin]; | |
| 912 } | |
| 913 else if( ijp[iin][jin] > 0 ) | |
| 914 { | |
| 915 ifi = iin-ijp[iin][jin]; jfi = jin-1; | |
| 916 } | |
| 917 else | |
| 918 { | |
| 919 ifi = iin-1; jfi = jin-1; | |
| 920 } | |
| 921 | |
| 922 if( ifi == -warpbase && jfi == -warpbase ) | |
| 923 { | |
| 924 l = iin; | |
| 925 while( --l >= 0 ) | |
| 926 { | |
| 927 *--gaptable1 = 'o'; | |
| 928 *--gaptable2 = '-'; | |
| 929 k++; | |
| 930 } | |
| 931 l= jin; | |
| 932 while( --l >= 0 ) | |
| 933 { | |
| 934 *--gaptable1 = '-'; | |
| 935 *--gaptable2 = 'o'; | |
| 936 } | |
| 937 break; | |
| 938 } | |
| 939 else | |
| 940 { | |
| 941 l = iin - ifi; | |
| 942 while( --l ) | |
| 943 { | |
| 944 *--gaptable1 = 'o'; | |
| 945 *--gaptable2 = '-'; | |
| 946 k++; | |
| 947 } | |
| 948 l= jin - jfi; | |
| 949 while( --l ) | |
| 950 { | |
| 951 *--gaptable1 = '-'; | |
| 952 *--gaptable2 = 'o'; | |
| 953 k++; | |
| 954 } | |
| 955 } | |
| 956 if( iin <= 0 || jin <= 0 ) break; | |
| 957 *--gaptable1 = 'o'; | |
| 958 *--gaptable2 = 'o'; | |
| 959 k++; | |
| 960 iin = ifi; jin = jfi; | |
| 961 } | |
| 962 | |
| 963 for( i=0; i<icyc; i++ ) gapireru( mseq1[i], seq1[i], gaptable1 ); | |
| 964 for( j=0; j<jcyc; j++ ) gapireru( mseq2[j], seq2[j], gaptable2 ); | |
| 965 | |
| 966 free( gt1bk ); | |
| 967 free( gt2bk ); | |
| 968 | |
| 969 return( 0.0 ); | |
| 970 } | |
| 971 | |
| 972 float A__align( double **n_dynamicmtx, char **seq1, char **seq2, double *eff1, double *eff2, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp ) | |
| 973 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */ | |
| 974 { | |
| 975 // int k; | |
| 976 register int i, j; | |
| 977 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */ | |
| 978 int lgth1, lgth2; | |
| 979 int resultlen; | |
| 980 float wm = 0.0; /* int ?????? */ | |
| 981 float g; | |
| 982 float *currentw, *previousw; | |
| 983 // float fpenalty = (float)penalty; | |
| 984 #if USE_PENALTY_EX | |
| 985 float fpenalty_ex = (float)penalty_ex; | |
| 986 #endif | |
| 987 #if 1 | |
| 988 float *wtmp; | |
| 989 int *ijppt; | |
| 990 float *mjpt, *prept, *curpt; | |
| 991 int *mpjpt; | |
| 992 #endif | |
| 993 static TLS float mi, *m; | |
| 994 static TLS int **ijp; | |
| 995 static TLS int mpi, *mp; | |
| 996 static TLS float *w1, *w2; | |
| 997 static TLS float *match; | |
| 998 static TLS float *initverticalw; /* kufuu sureba iranai */ | |
| 999 static TLS float *lastverticalw; /* kufuu sureba iranai */ | |
| 1000 static TLS char **mseq1; | |
| 1001 static TLS char **mseq2; | |
| 1002 static TLS char **mseq; | |
| 1003 static TLS float *ogcp1; | |
| 1004 static TLS float *ogcp2; | |
| 1005 static TLS float *fgcp1; | |
| 1006 static TLS float *fgcp2; | |
| 1007 static TLS float **cpmx1; | |
| 1008 static TLS float **cpmx2; | |
| 1009 static TLS int **intwork; | |
| 1010 static TLS float **floatwork; | |
| 1011 static TLS int orlgth1 = 0, orlgth2 = 0; | |
| 1012 static TLS float *gapfreq1; | |
| 1013 static TLS float *gapfreq2; | |
| 1014 float fpenalty = (float)penalty; | |
| 1015 float fpenalty_shift = (float)penalty_shift; | |
| 1016 float *fgcp2pt; | |
| 1017 float *ogcp2pt; | |
| 1018 float fgcp1va; | |
| 1019 float ogcp1va; | |
| 1020 float *gf2pt; | |
| 1021 float *gf2ptpre; | |
| 1022 float gf1va; | |
| 1023 float gf1vapre; | |
| 1024 float headgapfreq1; | |
| 1025 float headgapfreq2; | |
| 1026 | |
| 1027 int *warpis = NULL; | |
| 1028 int *warpjs = NULL; | |
| 1029 int *warpi = NULL; | |
| 1030 int *warpj = NULL; | |
| 1031 int *prevwarpi = NULL; | |
| 1032 int *prevwarpj = NULL; | |
| 1033 float *wmrecords = NULL; | |
| 1034 float *prevwmrecords = NULL; | |
| 1035 int warpn = 0; | |
| 1036 int warpbase; | |
| 1037 float curm = 0.0; | |
| 1038 float *wmrecordspt, *wmrecords1pt, *prevwmrecordspt; | |
| 1039 int *warpipt, *warpjpt; | |
| 1040 | |
| 1041 | |
| 1042 if( seq1 == NULL ) | |
| 1043 { | |
| 1044 if( orlgth1 ) | |
| 1045 { | |
| 1046 // fprintf( stderr, "## Freeing local arrays in A__align\n" ); | |
| 1047 orlgth1 = 0; | |
| 1048 orlgth2 = 0; | |
| 1049 | |
| 1050 imp_match_init_strict( NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0 ); | |
| 1051 | |
| 1052 free( mseq1 ); | |
| 1053 free( mseq2 ); | |
| 1054 FreeFloatVec( w1 ); | |
| 1055 FreeFloatVec( w2 ); | |
| 1056 FreeFloatVec( match ); | |
| 1057 FreeFloatVec( initverticalw ); | |
| 1058 FreeFloatVec( lastverticalw ); | |
| 1059 | |
| 1060 FreeFloatVec( m ); | |
| 1061 FreeIntVec( mp ); | |
| 1062 | |
| 1063 FreeCharMtx( mseq ); | |
| 1064 | |
| 1065 FreeFloatVec( ogcp1 ); | |
| 1066 FreeFloatVec( ogcp2 ); | |
| 1067 FreeFloatVec( fgcp1 ); | |
| 1068 FreeFloatVec( fgcp2 ); | |
| 1069 | |
| 1070 | |
| 1071 FreeFloatMtx( cpmx1 ); | |
| 1072 FreeFloatMtx( cpmx2 ); | |
| 1073 | |
| 1074 FreeFloatVec( gapfreq1 ); | |
| 1075 FreeFloatVec( gapfreq2 ); | |
| 1076 | |
| 1077 FreeFloatMtx( floatwork ); | |
| 1078 FreeIntMtx( intwork ); | |
| 1079 | |
| 1080 } | |
| 1081 else | |
| 1082 { | |
| 1083 // fprintf( stderr, "## Not allocated\n" ); | |
| 1084 } | |
| 1085 return( 0.0 ); | |
| 1086 } | |
| 1087 | |
| 1088 | |
| 1089 lgth1 = strlen( seq1[0] ); | |
| 1090 lgth2 = strlen( seq2[0] ); | |
| 1091 #if 0 | |
| 1092 if( lgth1 == 0 || lgth2 == 0 ) | |
| 1093 { | |
| 1094 fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 ); | |
| 1095 } | |
| 1096 #endif | |
| 1097 if( lgth1 == 0 && lgth2 == 0 ) | |
| 1098 return( 0.0 ); | |
| 1099 | |
| 1100 if( lgth1 == 0 ) | |
| 1101 { | |
| 1102 for( i=0; i<icyc; i++ ) | |
| 1103 { | |
| 1104 j = lgth2; | |
| 1105 seq1[i][j] = 0; | |
| 1106 while( j ) seq1[i][--j] = *newgapstr; | |
| 1107 // fprintf( stderr, "seq1[i] = %s\n", seq1[i] ); | |
| 1108 } | |
| 1109 return( 0.0 ); | |
| 1110 } | |
| 1111 | |
| 1112 if( lgth2 == 0 ) | |
| 1113 { | |
| 1114 for( i=0; i<jcyc; i++ ) | |
| 1115 { | |
| 1116 j = lgth1; | |
| 1117 seq2[i][j] = 0; | |
| 1118 while( j ) seq2[i][--j] = *newgapstr; | |
| 1119 // fprintf( stderr, "seq2[i] = %s\n", seq2[i] ); | |
| 1120 } | |
| 1121 return( 0.0 ); | |
| 1122 } | |
| 1123 | |
| 1124 warpbase = lgth1 + lgth2; | |
| 1125 warpis = NULL; | |
| 1126 warpjs = NULL; | |
| 1127 warpn = 0; | |
| 1128 | |
| 1129 | |
| 1130 | |
| 1131 if( trywarp ) | |
| 1132 { | |
| 1133 // fprintf( stderr, "IN A__align, penalty_shift = %d\n", penalty_shift ); | |
| 1134 if( headgp == 0 || tailgp == 0 ) | |
| 1135 { | |
| 1136 fprintf( stderr, "At present, headgp and tailgp must be 1 to allow shift.\n" ); | |
| 1137 exit( 1 ); | |
| 1138 } | |
| 1139 wmrecords = AllocateFloatVec( lgth2+1 ); | |
| 1140 warpi = AllocateIntVec( lgth2+1 ); | |
| 1141 warpj = AllocateIntVec( lgth2+1 ); | |
| 1142 prevwmrecords = AllocateFloatVec( lgth2+1 ); | |
| 1143 prevwarpi = AllocateIntVec( lgth2+1 ); | |
| 1144 prevwarpj = AllocateIntVec( lgth2+1 ); | |
| 1145 for( i=0; i<lgth2+1; i++ ) wmrecords[i] = 0.0; | |
| 1146 for( i=0; i<lgth2+1; i++ ) prevwmrecords[i] = 0.0; | |
| 1147 for( i=0; i<lgth2+1; i++ ) prevwarpi[i] = -warpbase; | |
| 1148 for( i=0; i<lgth2+1; i++ ) prevwarpj[i] = -warpbase; | |
| 1149 for( i=0; i<lgth2+1; i++ ) warpi[i] = -warpbase; | |
| 1150 for( i=0; i<lgth2+1; i++ ) warpj[i] = -warpbase; | |
| 1151 } | |
| 1152 | |
| 1153 | |
| 1154 #if 0 | |
| 1155 fprintf( stderr, "#### eff in SA+++align\n" ); | |
| 1156 fprintf( stderr, "#### seq1[0] = %s\n", seq1[0] ); | |
| 1157 fprintf( stderr, "#### strlen( seq1[0] ) = %d\n", strlen( seq1[0] ) ); | |
| 1158 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] ); | |
| 1159 fprintf( stderr, "#### seq2[0] = %s\n", seq2[0] ); | |
| 1160 fprintf( stderr, "#### strlen( seq2[0] ) = %d\n", strlen( seq2[0] ) ); | |
| 1161 for( i=0; i<jcyc; i++ ) fprintf( stderr, "eff2[%d] = %f\n", i, eff2[i] ); | |
| 1162 #endif | |
| 1163 if( orlgth1 == 0 ) | |
| 1164 { | |
| 1165 mseq1 = AllocateCharMtx( njob, 0 ); | |
| 1166 mseq2 = AllocateCharMtx( njob, 0 ); | |
| 1167 } | |
| 1168 | |
| 1169 | |
| 1170 | |
| 1171 if( lgth1 > orlgth1 || lgth2 > orlgth2 ) | |
| 1172 { | |
| 1173 int ll1, ll2; | |
| 1174 | |
| 1175 | |
| 1176 if( orlgth1 > 0 && orlgth2 > 0 ) | |
| 1177 { | |
| 1178 FreeFloatVec( w1 ); | |
| 1179 FreeFloatVec( w2 ); | |
| 1180 FreeFloatVec( match ); | |
| 1181 FreeFloatVec( initverticalw ); | |
| 1182 FreeFloatVec( lastverticalw ); | |
| 1183 | |
| 1184 FreeFloatVec( m ); | |
| 1185 FreeIntVec( mp ); | |
| 1186 | |
| 1187 FreeCharMtx( mseq ); | |
| 1188 | |
| 1189 FreeFloatVec( ogcp1 ); | |
| 1190 FreeFloatVec( ogcp2 ); | |
| 1191 FreeFloatVec( fgcp1 ); | |
| 1192 FreeFloatVec( fgcp2 ); | |
| 1193 | |
| 1194 | |
| 1195 FreeFloatMtx( cpmx1 ); | |
| 1196 FreeFloatMtx( cpmx2 ); | |
| 1197 | |
| 1198 FreeFloatVec( gapfreq1 ); | |
| 1199 FreeFloatVec( gapfreq2 ); | |
| 1200 | |
| 1201 FreeFloatMtx( floatwork ); | |
| 1202 FreeIntMtx( intwork ); | |
| 1203 } | |
| 1204 | |
| 1205 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100; | |
| 1206 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100; | |
| 1207 | |
| 1208 #if DEBUG | |
| 1209 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 ); | |
| 1210 #endif | |
| 1211 | |
| 1212 w1 = AllocateFloatVec( ll2+2 ); | |
| 1213 w2 = AllocateFloatVec( ll2+2 ); | |
| 1214 match = AllocateFloatVec( ll2+2 ); | |
| 1215 | |
| 1216 initverticalw = AllocateFloatVec( ll1+2 ); | |
| 1217 lastverticalw = AllocateFloatVec( ll1+2 ); | |
| 1218 | |
| 1219 m = AllocateFloatVec( ll2+2 ); | |
| 1220 mp = AllocateIntVec( ll2+2 ); | |
| 1221 | |
| 1222 mseq = AllocateCharMtx( njob, ll1+ll2 ); | |
| 1223 | |
| 1224 ogcp1 = AllocateFloatVec( ll1+2 ); | |
| 1225 ogcp2 = AllocateFloatVec( ll2+2 ); | |
| 1226 fgcp1 = AllocateFloatVec( ll1+2 ); | |
| 1227 fgcp2 = AllocateFloatVec( ll2+2 ); | |
| 1228 | |
| 1229 cpmx1 = AllocateFloatMtx( nalphabets, ll1+2 ); | |
| 1230 cpmx2 = AllocateFloatMtx( nalphabets, ll2+2 ); | |
| 1231 | |
| 1232 gapfreq1 = AllocateFloatVec( ll1+2 ); | |
| 1233 gapfreq2 = AllocateFloatVec( ll2+2 ); | |
| 1234 | |
| 1235 #if FASTMATCHCALC | |
| 1236 floatwork = AllocateFloatMtx( MAX( ll1, ll2 )+2, nalphabets ); | |
| 1237 intwork = AllocateIntMtx( MAX( ll1, ll2 )+2, nalphabets+1 ); | |
| 1238 #else | |
| 1239 floatwork = AllocateFloatMtx( nalphabets, MAX( ll1, ll2 )+2 ); | |
| 1240 intwork = AllocateIntMtx( nalphabets, MAX( ll1, ll2 )+2 ); | |
| 1241 #endif | |
| 1242 | |
| 1243 #if DEBUG | |
| 1244 fprintf( stderr, "succeeded\n" ); | |
| 1245 #endif | |
| 1246 | |
| 1247 orlgth1 = ll1 - 100; | |
| 1248 orlgth2 = ll2 - 100; | |
| 1249 } | |
| 1250 | |
| 1251 | |
| 1252 for( i=0; i<icyc; i++ ) | |
| 1253 { | |
| 1254 mseq1[i] = mseq[i]; | |
| 1255 seq1[i][lgth1] = 0; | |
| 1256 } | |
| 1257 for( j=0; j<jcyc; j++ ) | |
| 1258 { | |
| 1259 mseq2[j] = mseq[icyc+j]; | |
| 1260 seq2[j][lgth2] = 0; | |
| 1261 } | |
| 1262 | |
| 1263 | |
| 1264 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 ) | |
| 1265 { | |
| 1266 int ll1, ll2; | |
| 1267 | |
| 1268 if( commonAlloc1 && commonAlloc2 ) | |
| 1269 { | |
| 1270 FreeIntMtx( commonIP ); | |
| 1271 } | |
| 1272 | |
| 1273 ll1 = MAX( orlgth1, commonAlloc1 ); | |
| 1274 ll2 = MAX( orlgth2, commonAlloc2 ); | |
| 1275 | |
| 1276 #if DEBUG | |
| 1277 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 ); | |
| 1278 #endif | |
| 1279 | |
| 1280 commonIP = AllocateIntMtx( ll1+10, ll2+10 ); | |
| 1281 | |
| 1282 #if DEBUG | |
| 1283 fprintf( stderr, "succeeded\n\n" ); | |
| 1284 #endif | |
| 1285 | |
| 1286 commonAlloc1 = ll1; | |
| 1287 commonAlloc2 = ll2; | |
| 1288 } | |
| 1289 ijp = commonIP; | |
| 1290 | |
| 1291 #if 0 | |
| 1292 { | |
| 1293 float t = 0.0; | |
| 1294 for( i=0; i<icyc; i++ ) | |
| 1295 t += eff1[i]; | |
| 1296 fprintf( stderr, "## totaleff = %f\n", t ); | |
| 1297 } | |
| 1298 #endif | |
| 1299 | |
| 1300 cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc ); | |
| 1301 cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc ); | |
| 1302 | |
| 1303 if( sgap1 ) | |
| 1304 { | |
| 1305 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 ); | |
| 1306 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 ); | |
| 1307 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 ); | |
| 1308 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 ); | |
| 1309 outgapcount( &headgapfreq1, icyc, sgap1, eff1 ); | |
| 1310 outgapcount( &headgapfreq2, jcyc, sgap2, eff2 ); | |
| 1311 outgapcount( gapfreq1+lgth1, icyc, egap1, eff1 ); | |
| 1312 outgapcount( gapfreq2+lgth2, jcyc, egap2, eff2 ); | |
| 1313 } | |
| 1314 else | |
| 1315 { | |
| 1316 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 ); | |
| 1317 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 ); | |
| 1318 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 ); | |
| 1319 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 ); | |
| 1320 headgapfreq1 = 0.0; | |
| 1321 headgapfreq2 = 0.0; | |
| 1322 gapfreq1[lgth1] = 0.0; | |
| 1323 gapfreq2[lgth2] = 0.0; | |
| 1324 } | |
| 1325 | |
| 1326 if( legacygapcost == 0 ) | |
| 1327 { | |
| 1328 gapcountf( gapfreq1, seq1, icyc, eff1, lgth1 ); | |
| 1329 gapcountf( gapfreq2, seq2, jcyc, eff2, lgth2 ); | |
| 1330 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0 - gapfreq1[i]; | |
| 1331 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0 - gapfreq2[i]; | |
| 1332 headgapfreq1 = 1.0 - headgapfreq1; | |
| 1333 headgapfreq2 = 1.0 - headgapfreq2; | |
| 1334 } | |
| 1335 else | |
| 1336 { | |
| 1337 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0; | |
| 1338 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0; | |
| 1339 headgapfreq1 = 1.0; | |
| 1340 headgapfreq2 = 1.0; | |
| 1341 } | |
| 1342 | |
| 1343 #if 0 | |
| 1344 fprintf( stderr, "\ngapfreq1[] =" ); | |
| 1345 for( i=0; i<lgth1; i++ ) fprintf( stderr, "%5.2f ", gapfreq1[i] ); | |
| 1346 fprintf( stderr, "\n" ); | |
| 1347 | |
| 1348 fprintf( stderr, "\ngapfreq2[] =" ); | |
| 1349 for( i=0; i<lgth2; i++ ) fprintf( stderr, "%5.2f ", gapfreq2[i] ); | |
| 1350 fprintf( stderr, "\n" ); | |
| 1351 #endif | |
| 1352 | |
| 1353 | |
| 1354 for( i=0; i<lgth1; i++ ) | |
| 1355 { | |
| 1356 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty * ( gapfreq1[i] ); | |
| 1357 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty * ( gapfreq1[i] ); | |
| 1358 } | |
| 1359 | |
| 1360 for( i=0; i<lgth2; i++ ) | |
| 1361 { | |
| 1362 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty * ( gapfreq2[i] ); | |
| 1363 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty * ( gapfreq2[i] ); | |
| 1364 } | |
| 1365 #if 0 | |
| 1366 for( i=0; i<lgth1; i++ ) | |
| 1367 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] ); | |
| 1368 #endif | |
| 1369 | |
| 1370 currentw = w1; | |
| 1371 previousw = w2; | |
| 1372 | |
| 1373 match_calc( n_dynamicmtx, initverticalw, cpmx2, cpmx1, 0, lgth1, floatwork, intwork, 1 ); | |
| 1374 if( localhom ) | |
| 1375 imp_match_out_vead_tate( initverticalw, 0, lgth1 ); // 060306 | |
| 1376 | |
| 1377 match_calc( n_dynamicmtx, currentw, cpmx1, cpmx2, 0, lgth2, floatwork, intwork, 1 ); | |
| 1378 if( localhom ) | |
| 1379 imp_match_out_vead( currentw, 0, lgth2 ); // 060306 | |
| 1380 #if 0 // -> tbfast.c | |
| 1381 if( localhom ) | |
| 1382 imp_match_calc( n_dynamicmtx, currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 ); | |
| 1383 | |
| 1384 #endif | |
| 1385 | |
| 1386 if( headgp == 1 ) | |
| 1387 { | |
| 1388 for( i=1; i<lgth1+1; i++ ) | |
| 1389 { | |
| 1390 // initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ; | |
| 1391 initverticalw[i] += ( ogcp1[0] * headgapfreq2 + fgcp1[i-1] * gapfreq2[0] ) ; | |
| 1392 } | |
| 1393 for( j=1; j<lgth2+1; j++ ) | |
| 1394 { | |
| 1395 // currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ; | |
| 1396 currentw[j] += ( ogcp2[0] * headgapfreq1 + fgcp2[j-1] * gapfreq1[0] ) ; | |
| 1397 } | |
| 1398 } | |
| 1399 #if OUTGAP0TRY | |
| 1400 else | |
| 1401 { | |
| 1402 fprintf( stderr, "offset = %d\n", offset ); | |
| 1403 for( j=1; j<lgth2+1; j++ ) | |
| 1404 currentw[j] -= offset * j / 2.0; | |
| 1405 for( i=1; i<lgth1+1; i++ ) | |
| 1406 initverticalw[i] -= offset * i / 2.0; | |
| 1407 } | |
| 1408 #endif | |
| 1409 #if 0 | |
| 1410 fprintf( stderr, "\n " ); | |
| 1411 for( j=0; j<lgth2+1; j++ ) fprintf( stderr, " %c ", seq2[0][j] ); | |
| 1412 fprintf( stderr, "\n%c ", seq1[0][0] ); | |
| 1413 for( j=0; j<lgth2+1; j++ ) | |
| 1414 { | |
| 1415 fprintf( stderr, "%5.0f ", currentw[j] ); | |
| 1416 } | |
| 1417 fprintf( stderr, "\n" ); | |
| 1418 #endif | |
| 1419 | |
| 1420 | |
| 1421 for( j=1; j<lgth2+1; ++j ) | |
| 1422 { | |
| 1423 // m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0; | |
| 1424 m[j] = currentw[j-1] + ogcp1[1] * gapfreq2[j-1]; mp[j] = 0;; | |
| 1425 } | |
| 1426 if( lgth2 == 0 ) | |
| 1427 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari | |
| 1428 else | |
| 1429 lastverticalw[0] = currentw[lgth2-1]; | |
| 1430 | |
| 1431 if( tailgp ) lasti = lgth1+1; else lasti = lgth1; | |
| 1432 lastj = lgth2+1; | |
| 1433 | |
| 1434 #if XXXXXXX | |
| 1435 fprintf( stderr, "currentw = \n" ); | |
| 1436 for( i=0; i<lgth1+1; i++ ) | |
| 1437 { | |
| 1438 fprintf( stderr, "%5.2f ", currentw[i] ); | |
| 1439 } | |
| 1440 fprintf( stderr, "\n" ); | |
| 1441 fprintf( stderr, "initverticalw = \n" ); | |
| 1442 for( i=0; i<lgth2+1; i++ ) | |
| 1443 { | |
| 1444 fprintf( stderr, "%5.2f ", initverticalw[i] ); | |
| 1445 } | |
| 1446 fprintf( stderr, "\n" ); | |
| 1447 fprintf( stderr, "fcgp\n" ); | |
| 1448 for( i=0; i<lgth1; i++ ) | |
| 1449 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] ); | |
| 1450 for( i=0; i<lgth2; i++ ) | |
| 1451 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] ); | |
| 1452 #endif | |
| 1453 | |
| 1454 for( i=1; i<lasti; i++ ) | |
| 1455 { | |
| 1456 | |
| 1457 #ifdef enablemultithread | |
| 1458 // fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref ); | |
| 1459 if( chudanpt && *chudanpt != chudanref ) | |
| 1460 { | |
| 1461 // fprintf( stderr, "\n\n## CHUUDAN!!! S\n" ); | |
| 1462 *chudanres = 1; | |
| 1463 return( -1.0 ); | |
| 1464 } | |
| 1465 #endif | |
| 1466 wtmp = previousw; | |
| 1467 previousw = currentw; | |
| 1468 currentw = wtmp; | |
| 1469 | |
| 1470 previousw[0] = initverticalw[i-1]; | |
| 1471 | |
| 1472 match_calc( n_dynamicmtx, currentw, cpmx1, cpmx2, i, lgth2, floatwork, intwork, 0 ); | |
| 1473 #if XXXXXXX | |
| 1474 fprintf( stderr, "\n" ); | |
| 1475 fprintf( stderr, "i=%d\n", i ); | |
| 1476 fprintf( stderr, "currentw = \n" ); | |
| 1477 for( j=0; j<lgth2; j++ ) | |
| 1478 { | |
| 1479 fprintf( stderr, "%5.2f ", currentw[j] ); | |
| 1480 } | |
| 1481 fprintf( stderr, "\n" ); | |
| 1482 #endif | |
| 1483 if( localhom ) | |
| 1484 { | |
| 1485 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i ); | |
| 1486 #if 0 | |
| 1487 imp_match_out_vead( currentw, i, lgth2 ); | |
| 1488 #else | |
| 1489 imp_match_out_vead( currentw, i, lgth2 ); | |
| 1490 #endif | |
| 1491 } | |
| 1492 #if XXXXXXX | |
| 1493 fprintf( stderr, "\n" ); | |
| 1494 fprintf( stderr, "i=%d\n", i ); | |
| 1495 fprintf( stderr, "currentw = \n" ); | |
| 1496 for( j=0; j<lgth2; j++ ) | |
| 1497 { | |
| 1498 fprintf( stderr, "%5.2f ", currentw[j] ); | |
| 1499 } | |
| 1500 fprintf( stderr, "\n" ); | |
| 1501 #endif | |
| 1502 currentw[0] = initverticalw[i]; | |
| 1503 | |
| 1504 #if 0 | |
| 1505 fprintf( stderr, "%c ", seq1[0][i] ); | |
| 1506 for( j=0; j<lgth2+1; j++ ) | |
| 1507 { | |
| 1508 fprintf( stderr, "%5.0f ", currentw[j] ); | |
| 1509 } | |
| 1510 fprintf( stderr, "\n" ); | |
| 1511 #endif | |
| 1512 | |
| 1513 // mi = previousw[0] + ogcp2[1]; mpi = 0; | |
| 1514 mi = previousw[0] + ogcp2[1] * gapfreq1[i-1]; mpi=0; | |
| 1515 ijppt = ijp[i] + 1; | |
| 1516 mjpt = m + 1; | |
| 1517 prept = previousw; | |
| 1518 curpt = currentw + 1; | |
| 1519 mpjpt = mp + 1; | |
| 1520 fgcp2pt = fgcp2; | |
| 1521 ogcp2pt = ogcp2 + 1; | |
| 1522 fgcp1va = fgcp1[i-1]; | |
| 1523 ogcp1va = ogcp1[i]; | |
| 1524 gf1va = gapfreq1[i]; | |
| 1525 gf1vapre = gapfreq1[i-1]; | |
| 1526 gf2pt = gapfreq2+1; | |
| 1527 gf2ptpre = gapfreq2; | |
| 1528 | |
| 1529 if( trywarp ) | |
| 1530 { | |
| 1531 prevwmrecordspt = prevwmrecords; | |
| 1532 wmrecordspt = wmrecords+1; | |
| 1533 wmrecords1pt = wmrecords; | |
| 1534 warpipt = warpi + 1; | |
| 1535 warpjpt = warpj + 1; | |
| 1536 } | |
| 1537 | |
| 1538 | |
| 1539 for( j=1; j<lastj; j++ ) | |
| 1540 { | |
| 1541 #ifdef xxxenablemultithread | |
| 1542 // fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref ); | |
| 1543 if( chudanpt && *chudanpt != chudanref ) | |
| 1544 { | |
| 1545 // fprintf( stderr, "\n\n## CHUUDAN!!! S\n" ); | |
| 1546 *chudanres = 1; | |
| 1547 return( -1.0 ); | |
| 1548 } | |
| 1549 #endif | |
| 1550 wm = *prept; | |
| 1551 *ijppt = 0; | |
| 1552 | |
| 1553 #if 0 | |
| 1554 fprintf( stderr, "\n i=%d, j=%d %c, %c", i, j, seq1[0][i], seq2[0][j] ); | |
| 1555 fprintf( stderr, "%5.0f->", wm ); | |
| 1556 fprintf( stderr, "%5.0f? (penal=%5.2f)", g=mi+*fgcp2pt*(1.0-gapfreq1[i]), *fgcp2pt*(1.0-gapfreq1[i]) ); | |
| 1557 #endif | |
| 1558 if( (g=mi+*fgcp2pt*gf1va) > wm ) | |
| 1559 { | |
| 1560 wm = g; | |
| 1561 *ijppt = -( j - mpi ); | |
| 1562 // fprintf( stderr, "Jump to %d (%c)!", mpi, seq2[0][mpi] ); | |
| 1563 } | |
| 1564 if( (g=*prept+*ogcp2pt*gf1vapre) >= mi ) | |
| 1565 { | |
| 1566 mi = g; | |
| 1567 mpi = j-1; | |
| 1568 } | |
| 1569 #if USE_PENALTY_EX | |
| 1570 mi += fpenalty_ex; | |
| 1571 #endif | |
| 1572 | |
| 1573 #if 0 | |
| 1574 fprintf( stderr, "%5.0f->", wm ); | |
| 1575 fprintf( stderr, "%5.0f? (penal=%5.2f)", g=*mjpt+fgcp1va*(1.0-gapfreq2[j]), fgcp1va*(1.0-gapfreq2[j]) ); | |
| 1576 #endif | |
| 1577 if( (g=*mjpt+ fgcp1va* *gf2pt) > wm ) | |
| 1578 { | |
| 1579 wm = g; | |
| 1580 *ijppt = +( i - *mpjpt ); | |
| 1581 // fprintf( stderr, "Jump to %d (%c)!", *mpjpt, seq1[0][*mpjpt] ); | |
| 1582 } | |
| 1583 if( (g=*prept+ ogcp1va* *gf2ptpre) >= *mjpt ) | |
| 1584 { | |
| 1585 *mjpt = g; | |
| 1586 *mpjpt = i-1; | |
| 1587 } | |
| 1588 #if USE_PENALTY_EX | |
| 1589 m[j] += fpenalty_ex; | |
| 1590 #endif | |
| 1591 | |
| 1592 | |
| 1593 if( trywarp ) | |
| 1594 { | |
| 1595 #if USE_PENALTY_EX | |
| 1596 if( ( g=*prevwmrecordspt++ + fpenalty_shift + fpenalty_ex * ( i - prevwarpi[j-1] + j - prevwarpj[j-1] ) ) > wm ) // naka ha osokute kamawanai | |
| 1597 #else | |
| 1598 if( ( g=*prevwmrecordspt++ + fpenalty_shift ) > wm ) // naka ha osokute kamawanai | |
| 1599 #endif | |
| 1600 { | |
| 1601 // fprintf( stderr, "WARP!!\n" ); | |
| 1602 if( warpn && prevwarpi[j-1] == warpis[warpn-1] && prevwarpj[j-1] == warpjs[warpn-1] ) | |
| 1603 { | |
| 1604 *ijppt = warpbase + warpn - 1; | |
| 1605 } | |
| 1606 else | |
| 1607 { | |
| 1608 *ijppt = warpbase + warpn; | |
| 1609 warpis = realloc( warpis, sizeof(int) * ( warpn+1 ) ); | |
| 1610 warpjs = realloc( warpjs, sizeof(int) * ( warpn+1 ) ); | |
| 1611 warpis[warpn] = prevwarpi[j-1]; | |
| 1612 warpjs[warpn] = prevwarpj[j-1]; | |
| 1613 warpn++; | |
| 1614 } | |
| 1615 wm = g; | |
| 1616 } | |
| 1617 | |
| 1618 #if 0 | |
| 1619 fprintf( stderr, "%5.0f ", wm ); | |
| 1620 #endif | |
| 1621 curm = *curpt + wm; | |
| 1622 | |
| 1623 if( *wmrecords1pt > *wmrecordspt ) | |
| 1624 { | |
| 1625 *wmrecordspt = *wmrecords1pt; | |
| 1626 *warpipt = *(warpipt-1); | |
| 1627 *warpjpt = *(warpjpt-1); | |
| 1628 } | |
| 1629 if( curm > *wmrecordspt ) | |
| 1630 { | |
| 1631 *wmrecordspt = curm; | |
| 1632 *warpipt = i; | |
| 1633 *warpjpt = j; | |
| 1634 } | |
| 1635 wmrecordspt++; | |
| 1636 wmrecords1pt++; | |
| 1637 warpipt++; | |
| 1638 warpjpt++; | |
| 1639 } | |
| 1640 | |
| 1641 *curpt++ += wm; | |
| 1642 ijppt++; | |
| 1643 mjpt++; | |
| 1644 prept++; | |
| 1645 mpjpt++; | |
| 1646 fgcp2pt++; | |
| 1647 ogcp2pt++; | |
| 1648 gf2ptpre++; | |
| 1649 gf2pt++; | |
| 1650 } | |
| 1651 lastverticalw[i] = currentw[lgth2-1]; | |
| 1652 | |
| 1653 if( trywarp ) | |
| 1654 { | |
| 1655 fltncpy( prevwmrecords, wmrecords, lastj ); | |
| 1656 intncpy( prevwarpi, warpi, lastj ); | |
| 1657 intncpy( prevwarpj, warpj, lastj ); | |
| 1658 } | |
| 1659 } | |
| 1660 | |
| 1661 if( trywarp ) | |
| 1662 { | |
| 1663 // fprintf( stderr, "wm = %f\n", wm ); | |
| 1664 // fprintf( stderr, "warpn = %d\n", warpn ); | |
| 1665 free( wmrecords ); | |
| 1666 free( prevwmrecords ); | |
| 1667 free( warpi ); | |
| 1668 free( warpj ); | |
| 1669 free( prevwarpi ); | |
| 1670 free( prevwarpj ); | |
| 1671 } | |
| 1672 | |
| 1673 #if OUTGAP0TRY | |
| 1674 if( !outgap ) | |
| 1675 { | |
| 1676 for( j=1; j<lgth2+1; j++ ) | |
| 1677 currentw[j] -= offset * ( lgth2 - j ) / 2.0; | |
| 1678 for( i=1; i<lgth1+1; i++ ) | |
| 1679 lastverticalw[i] -= offset * ( lgth1 - i / 2.0); | |
| 1680 } | |
| 1681 #endif | |
| 1682 | |
| 1683 /* | |
| 1684 fprintf( stderr, "\n" ); | |
| 1685 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] ); | |
| 1686 fprintf( stderr, "#####\n" ); | |
| 1687 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] ); | |
| 1688 fprintf( stderr, "====>" ); | |
| 1689 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] ); | |
| 1690 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] ); | |
| 1691 */ | |
| 1692 if( localhom ) | |
| 1693 { | |
| 1694 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, warpis, warpjs, warpbase ); | |
| 1695 } | |
| 1696 else | |
| 1697 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, tailgp, warpis, warpjs, warpbase ); | |
| 1698 | |
| 1699 if( warpis ) free( warpis ); | |
| 1700 if( warpjs ) free( warpjs ); | |
| 1701 | |
| 1702 // fprintf( stderr, "### impmatch = %f\n", *impmatch ); | |
| 1703 | |
| 1704 resultlen = strlen( mseq1[0] ); | |
| 1705 if( alloclen < resultlen || resultlen > N ) | |
| 1706 { | |
| 1707 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N ); | |
| 1708 ErrorExit( "LENGTH OVER!\n" ); | |
| 1709 } | |
| 1710 | |
| 1711 | |
| 1712 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] ); | |
| 1713 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] ); | |
| 1714 #if 0 | |
| 1715 fprintf( stderr, "\n" ); | |
| 1716 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] ); | |
| 1717 fprintf( stderr, "#####\n" ); | |
| 1718 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] ); | |
| 1719 #endif | |
| 1720 | |
| 1721 // fprintf( stderr, "wm = %f\n", wm ); | |
| 1722 | |
| 1723 return( wm ); | |
| 1724 } | |
| 1725 | |
| 1726 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 ) | |
| 1727 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */ | |
| 1728 { | |
| 1729 fprintf( stderr, "Unexpected error. Please contact kazutaka.katoh@aist.go.jp\n" ); | |
| 1730 exit( 1 ); | |
| 1731 } | |
| 1732 | |
| 1733 | |
| 1734 float A__align_variousdist( double **pairoffset, double ***matrices, double **n_dynamicmtx, char **seq1, char **seq2, double *eff1, double *eff2, double **eff1s, double **eff2s, int icyc, int jcyc, int alloclen, LocalHom ***localhom, float *impmatch, char *sgap1, char *sgap2, char *egap1, char *egap2, int *chudanpt, int chudanref, int *chudanres, int headgp, int tailgp ) | |
| 1735 /* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */ | |
| 1736 { | |
| 1737 | |
| 1738 // int k; | |
| 1739 register int i, j, c; | |
| 1740 int lasti, lastj; /* outgap == 0 -> lgth1, outgap == 1 -> lgth1+1 */ | |
| 1741 int lgth1, lgth2; | |
| 1742 int resultlen; | |
| 1743 float wm = 0.0; /* int ?????? */ | |
| 1744 float g; | |
| 1745 float *currentw, *previousw; | |
| 1746 // float fpenalty = (float)penalty; | |
| 1747 #if USE_PENALTY_EX | |
| 1748 float fpenalty_ex = (float)penalty_ex; | |
| 1749 #endif | |
| 1750 #if 1 | |
| 1751 float *wtmp; | |
| 1752 int *ijppt; | |
| 1753 float *mjpt, *prept, *curpt; | |
| 1754 int *mpjpt; | |
| 1755 #endif | |
| 1756 static TLS float mi, *m; | |
| 1757 static TLS int **ijp; | |
| 1758 static TLS int mpi, *mp; | |
| 1759 static TLS float *w1, *w2; | |
| 1760 static TLS float *match; | |
| 1761 static TLS float *initverticalw; /* kufuu sureba iranai */ | |
| 1762 static TLS float *lastverticalw; /* kufuu sureba iranai */ | |
| 1763 static TLS char **mseq1; | |
| 1764 static TLS char **mseq2; | |
| 1765 static TLS char **mseq; | |
| 1766 static TLS float *ogcp1; | |
| 1767 static TLS float *ogcp2; | |
| 1768 static TLS float *fgcp1; | |
| 1769 static TLS float *fgcp2; | |
| 1770 static TLS float ***cpmx1s; | |
| 1771 static TLS float ***cpmx2s; | |
| 1772 static TLS int ***intwork; | |
| 1773 static TLS float ***floatwork; | |
| 1774 static TLS int orlgth1 = 0, orlgth2 = 0; | |
| 1775 static TLS float *gapfreq1; | |
| 1776 static TLS float *gapfreq2; | |
| 1777 float fpenalty = (float)penalty; | |
| 1778 float fpenalty_shift = (float)penalty_shift; | |
| 1779 float *fgcp2pt; | |
| 1780 float *ogcp2pt; | |
| 1781 float fgcp1va; | |
| 1782 float ogcp1va; | |
| 1783 float *gf2pt; | |
| 1784 float *gf2ptpre; | |
| 1785 float gf1va; | |
| 1786 float gf1vapre; | |
| 1787 float headgapfreq1; | |
| 1788 float headgapfreq2; | |
| 1789 | |
| 1790 int *warpis = NULL; | |
| 1791 int *warpjs = NULL; | |
| 1792 int *warpi = NULL; | |
| 1793 int *warpj = NULL; | |
| 1794 int *prevwarpi = NULL; | |
| 1795 int *prevwarpj = NULL; | |
| 1796 float *wmrecords = NULL; | |
| 1797 float *prevwmrecords = NULL; | |
| 1798 int warpn = 0; | |
| 1799 int warpbase; | |
| 1800 float curm = 0.0; | |
| 1801 float *wmrecordspt, *wmrecords1pt, *prevwmrecordspt; | |
| 1802 int *warpipt, *warpjpt; | |
| 1803 | |
| 1804 | |
| 1805 if( seq1 == NULL ) | |
| 1806 { | |
| 1807 if( orlgth1 ) | |
| 1808 { | |
| 1809 // fprintf( stderr, "## Freeing local arrays in A__align\n" ); | |
| 1810 orlgth1 = 0; | |
| 1811 orlgth2 = 0; | |
| 1812 | |
| 1813 imp_match_init_strict( NULL, 0, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0 ); | |
| 1814 | |
| 1815 free( mseq1 ); | |
| 1816 free( mseq2 ); | |
| 1817 FreeFloatVec( w1 ); | |
| 1818 FreeFloatVec( w2 ); | |
| 1819 FreeFloatVec( match ); | |
| 1820 FreeFloatVec( initverticalw ); | |
| 1821 FreeFloatVec( lastverticalw ); | |
| 1822 | |
| 1823 FreeFloatVec( m ); | |
| 1824 FreeIntVec( mp ); | |
| 1825 | |
| 1826 FreeCharMtx( mseq ); | |
| 1827 | |
| 1828 FreeFloatVec( ogcp1 ); | |
| 1829 FreeFloatVec( ogcp2 ); | |
| 1830 FreeFloatVec( fgcp1 ); | |
| 1831 FreeFloatVec( fgcp2 ); | |
| 1832 | |
| 1833 | |
| 1834 FreeFloatCub( cpmx1s ); | |
| 1835 FreeFloatCub( cpmx2s ); | |
| 1836 | |
| 1837 FreeFloatVec( gapfreq1 ); | |
| 1838 FreeFloatVec( gapfreq2 ); | |
| 1839 | |
| 1840 FreeFloatCub( floatwork ); | |
| 1841 FreeIntCub( intwork ); | |
| 1842 | |
| 1843 } | |
| 1844 else | |
| 1845 { | |
| 1846 // fprintf( stderr, "## Not allocated\n" ); | |
| 1847 } | |
| 1848 return( 0.0 ); | |
| 1849 } | |
| 1850 | |
| 1851 | |
| 1852 | |
| 1853 // fprintf( stderr, "\npairoffset = \n" ); | |
| 1854 // for( i=0; i<icyc; i ++ ) for( j=0; j<jcyc; j ++ ) | |
| 1855 // fprintf( stderr, "%d-%d, %f\n", i, j, pairoffset[i][j] ); | |
| 1856 | |
| 1857 lgth1 = strlen( seq1[0] ); | |
| 1858 lgth2 = strlen( seq2[0] ); | |
| 1859 #if 0 | |
| 1860 if( lgth1 == 0 || lgth2 == 0 ) | |
| 1861 { | |
| 1862 fprintf( stderr, "WARNING (Aalignmm): lgth1=%d, lgth2=%d\n", lgth1, lgth2 ); | |
| 1863 } | |
| 1864 #endif | |
| 1865 if( lgth1 == 0 && lgth2 == 0 ) | |
| 1866 return( 0.0 ); | |
| 1867 | |
| 1868 if( lgth1 == 0 ) | |
| 1869 { | |
| 1870 for( i=0; i<icyc; i++ ) | |
| 1871 { | |
| 1872 j = lgth2; | |
| 1873 seq1[i][j] = 0; | |
| 1874 while( j ) seq1[i][--j] = *newgapstr; | |
| 1875 // fprintf( stderr, "seq1[i] = %s\n", seq1[i] ); | |
| 1876 } | |
| 1877 return( 0.0 ); | |
| 1878 } | |
| 1879 | |
| 1880 if( lgth2 == 0 ) | |
| 1881 { | |
| 1882 for( i=0; i<jcyc; i++ ) | |
| 1883 { | |
| 1884 j = lgth1; | |
| 1885 seq2[i][j] = 0; | |
| 1886 while( j ) seq2[i][--j] = *newgapstr; | |
| 1887 // fprintf( stderr, "seq2[i] = %s\n", seq2[i] ); | |
| 1888 } | |
| 1889 return( 0.0 ); | |
| 1890 } | |
| 1891 | |
| 1892 warpbase = lgth1 + lgth2; | |
| 1893 warpis = NULL; | |
| 1894 warpjs = NULL; | |
| 1895 warpn = 0; | |
| 1896 | |
| 1897 if( trywarp ) | |
| 1898 { | |
| 1899 // fprintf( stderr, "In A__align_variousdist !!!!!\n" ); | |
| 1900 if( headgp == 0 || tailgp == 0 ) | |
| 1901 { | |
| 1902 fprintf( stderr, "At present, headgp and tailgp must be 1.\n" ); | |
| 1903 exit( 1 ); | |
| 1904 } | |
| 1905 wmrecords = AllocateFloatVec( lgth2+1 ); | |
| 1906 warpi = AllocateIntVec( lgth2+1 ); | |
| 1907 warpj = AllocateIntVec( lgth2+1 ); | |
| 1908 prevwmrecords = AllocateFloatVec( lgth2+1 ); | |
| 1909 prevwarpi = AllocateIntVec( lgth2+1 ); | |
| 1910 prevwarpj = AllocateIntVec( lgth2+1 ); | |
| 1911 for( i=0; i<lgth2+1; i++ ) wmrecords[i] = 0.0; | |
| 1912 for( i=0; i<lgth2+1; i++ ) prevwmrecords[i] = 0.0; | |
| 1913 for( i=0; i<lgth2+1; i++ ) prevwarpi[i] = -warpbase; | |
| 1914 for( i=0; i<lgth2+1; i++ ) prevwarpj[i] = -warpbase; | |
| 1915 for( i=0; i<lgth2+1; i++ ) warpi[i] = -warpbase; | |
| 1916 for( i=0; i<lgth2+1; i++ ) warpj[i] = -warpbase; | |
| 1917 } | |
| 1918 | |
| 1919 #if 0 | |
| 1920 fprintf( stderr, "#### eff in SA+++align\n" ); | |
| 1921 fprintf( stderr, "#### seq1[0] = %s\n", seq1[0] ); | |
| 1922 fprintf( stderr, "#### strlen( seq1[0] ) = %d\n", strlen( seq1[0] ) ); | |
| 1923 for( i=0; i<icyc; i++ ) fprintf( stderr, "eff1[%d] = %f\n", i, eff1[i] ); | |
| 1924 fprintf( stderr, "#### seq2[0] = %s\n", seq2[0] ); | |
| 1925 fprintf( stderr, "#### strlen( seq2[0] ) = %d\n", strlen( seq2[0] ) ); | |
| 1926 for( i=0; i<jcyc; i++ ) fprintf( stderr, "eff2[%d] = %f\n", i, eff2[i] ); | |
| 1927 #endif | |
| 1928 if( orlgth1 == 0 ) | |
| 1929 { | |
| 1930 mseq1 = AllocateCharMtx( njob, 0 ); | |
| 1931 mseq2 = AllocateCharMtx( njob, 0 ); | |
| 1932 } | |
| 1933 | |
| 1934 | |
| 1935 | |
| 1936 if( lgth1 > orlgth1 || lgth2 > orlgth2 ) | |
| 1937 { | |
| 1938 int ll1, ll2; | |
| 1939 | |
| 1940 | |
| 1941 if( orlgth1 > 0 && orlgth2 > 0 ) | |
| 1942 { | |
| 1943 FreeFloatVec( w1 ); | |
| 1944 FreeFloatVec( w2 ); | |
| 1945 FreeFloatVec( match ); | |
| 1946 FreeFloatVec( initverticalw ); | |
| 1947 FreeFloatVec( lastverticalw ); | |
| 1948 | |
| 1949 FreeFloatVec( m ); | |
| 1950 FreeIntVec( mp ); | |
| 1951 | |
| 1952 FreeCharMtx( mseq ); | |
| 1953 | |
| 1954 FreeFloatVec( ogcp1 ); | |
| 1955 FreeFloatVec( ogcp2 ); | |
| 1956 FreeFloatVec( fgcp1 ); | |
| 1957 FreeFloatVec( fgcp2 ); | |
| 1958 | |
| 1959 | |
| 1960 FreeFloatCub( cpmx1s ); | |
| 1961 FreeFloatCub( cpmx2s ); | |
| 1962 | |
| 1963 FreeFloatVec( gapfreq1 ); | |
| 1964 FreeFloatVec( gapfreq2 ); | |
| 1965 | |
| 1966 FreeFloatCub( floatwork ); | |
| 1967 FreeIntCub( intwork ); | |
| 1968 } | |
| 1969 | |
| 1970 ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100; | |
| 1971 ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100; | |
| 1972 | |
| 1973 #if DEBUG | |
| 1974 fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 ); | |
| 1975 #endif | |
| 1976 | |
| 1977 w1 = AllocateFloatVec( ll2+2 ); | |
| 1978 w2 = AllocateFloatVec( ll2+2 ); | |
| 1979 match = AllocateFloatVec( ll2+2 ); | |
| 1980 | |
| 1981 initverticalw = AllocateFloatVec( ll1+2 ); | |
| 1982 lastverticalw = AllocateFloatVec( ll1+2 ); | |
| 1983 | |
| 1984 m = AllocateFloatVec( ll2+2 ); | |
| 1985 mp = AllocateIntVec( ll2+2 ); | |
| 1986 | |
| 1987 mseq = AllocateCharMtx( njob, ll1+ll2 ); | |
| 1988 | |
| 1989 ogcp1 = AllocateFloatVec( ll1+2 ); | |
| 1990 ogcp2 = AllocateFloatVec( ll2+2 ); | |
| 1991 fgcp1 = AllocateFloatVec( ll1+2 ); | |
| 1992 fgcp2 = AllocateFloatVec( ll2+2 ); | |
| 1993 | |
| 1994 cpmx1s = AllocateFloatCub( maxdistclass, nalphabets, ll1+2 ); | |
| 1995 cpmx2s = AllocateFloatCub( maxdistclass, nalphabets, ll2+2 ); | |
| 1996 | |
| 1997 gapfreq1 = AllocateFloatVec( ll1+2 ); | |
| 1998 gapfreq2 = AllocateFloatVec( ll2+2 ); | |
| 1999 | |
| 2000 floatwork = AllocateFloatCub( maxdistclass, MAX( ll1, ll2 )+2, nalphabets ); | |
| 2001 intwork = AllocateIntCub( maxdistclass, MAX( ll1, ll2 )+2, nalphabets+1 ); | |
| 2002 | |
| 2003 #if DEBUG | |
| 2004 fprintf( stderr, "succeeded\n" ); | |
| 2005 #endif | |
| 2006 | |
| 2007 orlgth1 = ll1 - 100; | |
| 2008 orlgth2 = ll2 - 100; | |
| 2009 } | |
| 2010 | |
| 2011 | |
| 2012 for( i=0; i<icyc; i++ ) | |
| 2013 { | |
| 2014 mseq1[i] = mseq[i]; | |
| 2015 seq1[i][lgth1] = 0; | |
| 2016 } | |
| 2017 for( j=0; j<jcyc; j++ ) | |
| 2018 { | |
| 2019 mseq2[j] = mseq[icyc+j]; | |
| 2020 seq2[j][lgth2] = 0; | |
| 2021 } | |
| 2022 | |
| 2023 | |
| 2024 if( orlgth1 > commonAlloc1 || orlgth2 > commonAlloc2 ) | |
| 2025 { | |
| 2026 int ll1, ll2; | |
| 2027 | |
| 2028 if( commonAlloc1 && commonAlloc2 ) | |
| 2029 { | |
| 2030 FreeIntMtx( commonIP ); | |
| 2031 } | |
| 2032 | |
| 2033 ll1 = MAX( orlgth1, commonAlloc1 ); | |
| 2034 ll2 = MAX( orlgth2, commonAlloc2 ); | |
| 2035 | |
| 2036 #if DEBUG | |
| 2037 fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 ); | |
| 2038 #endif | |
| 2039 | |
| 2040 commonIP = AllocateIntMtx( ll1+10, ll2+10 ); | |
| 2041 | |
| 2042 #if DEBUG | |
| 2043 fprintf( stderr, "succeeded\n\n" ); | |
| 2044 #endif | |
| 2045 | |
| 2046 commonAlloc1 = ll1; | |
| 2047 commonAlloc2 = ll2; | |
| 2048 } | |
| 2049 ijp = commonIP; | |
| 2050 | |
| 2051 #if 0 | |
| 2052 { | |
| 2053 float t = 0.0; | |
| 2054 for( i=0; i<icyc; i++ ) | |
| 2055 t += eff1[i]; | |
| 2056 fprintf( stderr, "## totaleff = %f\n", t ); | |
| 2057 } | |
| 2058 #endif | |
| 2059 | |
| 2060 // cpmx_calc_new( seq1, cpmx1, eff1, lgth1, icyc ); | |
| 2061 // cpmx_calc_new( seq2, cpmx2, eff2, lgth2, jcyc ); | |
| 2062 for( c=0; c<maxdistclass; c++ ) | |
| 2063 { | |
| 2064 // fprintf( stderr, "c=%d\n", c ); | |
| 2065 cpmx_calc_new( seq1, cpmx1s[c], eff1s[c], lgth1, icyc ); | |
| 2066 cpmx_calc_new( seq2, cpmx2s[c], eff2s[c], lgth2, jcyc ); | |
| 2067 } | |
| 2068 | |
| 2069 if( sgap1 ) | |
| 2070 { | |
| 2071 new_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1, sgap1 ); | |
| 2072 new_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2, sgap2 ); | |
| 2073 new_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1, egap1 ); | |
| 2074 new_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2, egap2 ); | |
| 2075 outgapcount( &headgapfreq1, icyc, sgap1, eff1 ); | |
| 2076 outgapcount( &headgapfreq2, jcyc, sgap2, eff2 ); | |
| 2077 outgapcount( gapfreq1+lgth1, icyc, egap1, eff1 ); | |
| 2078 outgapcount( gapfreq2+lgth2, jcyc, egap2, eff2 ); | |
| 2079 } | |
| 2080 else | |
| 2081 { | |
| 2082 st_OpeningGapCount( ogcp1, icyc, seq1, eff1, lgth1 ); | |
| 2083 st_OpeningGapCount( ogcp2, jcyc, seq2, eff2, lgth2 ); | |
| 2084 st_FinalGapCount( fgcp1, icyc, seq1, eff1, lgth1 ); | |
| 2085 st_FinalGapCount( fgcp2, jcyc, seq2, eff2, lgth2 ); | |
| 2086 headgapfreq1 = 0.0; | |
| 2087 headgapfreq2 = 0.0; | |
| 2088 gapfreq1[lgth1] = 0.0; | |
| 2089 gapfreq2[lgth2] = 0.0; | |
| 2090 } | |
| 2091 | |
| 2092 if( legacygapcost == 0 ) | |
| 2093 { | |
| 2094 gapcountf( gapfreq1, seq1, icyc, eff1, lgth1 ); | |
| 2095 gapcountf( gapfreq2, seq2, jcyc, eff2, lgth2 ); | |
| 2096 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0 - gapfreq1[i]; | |
| 2097 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0 - gapfreq2[i]; | |
| 2098 headgapfreq1 = 1.0 - headgapfreq1; | |
| 2099 headgapfreq2 = 1.0 - headgapfreq2; | |
| 2100 } | |
| 2101 else | |
| 2102 { | |
| 2103 for( i=0; i<lgth1+1; i++ ) gapfreq1[i] = 1.0; | |
| 2104 for( i=0; i<lgth2+1; i++ ) gapfreq2[i] = 1.0; | |
| 2105 headgapfreq1 = 1.0; | |
| 2106 headgapfreq2 = 1.0; | |
| 2107 } | |
| 2108 | |
| 2109 #if 0 | |
| 2110 fprintf( stderr, "\ngapfreq1[] =" ); | |
| 2111 for( i=0; i<lgth1; i++ ) fprintf( stderr, "%5.2f ", gapfreq1[i] ); | |
| 2112 fprintf( stderr, "\n" ); | |
| 2113 | |
| 2114 fprintf( stderr, "\ngapfreq2[] =" ); | |
| 2115 for( i=0; i<lgth2; i++ ) fprintf( stderr, "%5.2f ", gapfreq2[i] ); | |
| 2116 fprintf( stderr, "\n" ); | |
| 2117 #endif | |
| 2118 | |
| 2119 | |
| 2120 for( i=0; i<lgth1; i++ ) | |
| 2121 { | |
| 2122 ogcp1[i] = 0.5 * ( 1.0 - ogcp1[i] ) * fpenalty * ( gapfreq1[i] ); | |
| 2123 fgcp1[i] = 0.5 * ( 1.0 - fgcp1[i] ) * fpenalty * ( gapfreq1[i] ); | |
| 2124 } | |
| 2125 | |
| 2126 for( i=0; i<lgth2; i++ ) | |
| 2127 { | |
| 2128 ogcp2[i] = 0.5 * ( 1.0 - ogcp2[i] ) * fpenalty * ( gapfreq2[i] ); | |
| 2129 fgcp2[i] = 0.5 * ( 1.0 - fgcp2[i] ) * fpenalty * ( gapfreq2[i] ); | |
| 2130 } | |
| 2131 #if 0 | |
| 2132 for( i=0; i<lgth1; i++ ) | |
| 2133 fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] ); | |
| 2134 #endif | |
| 2135 | |
| 2136 currentw = w1; | |
| 2137 previousw = w2; | |
| 2138 | |
| 2139 // for( i=0; i<icyc; i++ ) fprintf( stderr, "seq1[i] = %s\n", seq1[i] ); | |
| 2140 // for( j=0; j<jcyc; j++ ) fprintf( stderr, "seq2[j] = %s\n", seq2[j] ); | |
| 2141 | |
| 2142 #if MCD | |
| 2143 match_calc_dynamicmtx( pairoffset, n_dynamicmtx, initverticalw, jcyc, seq2, eff2, icyc, seq1, eff1, 0, lgth1, *floatwork, *intwork, 1, 1 ); | |
| 2144 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "%d - %f\n", i, initverticalw[i] ); | |
| 2145 #else | |
| 2146 fillzero( initverticalw, lgth1 ); | |
| 2147 for( c=0; c<maxdistclass; c++ ) | |
| 2148 { | |
| 2149 // fprintf( stderr, "c=%d matrices[c][W][W] = %f\n", c, matrices[c][amino_n['W']][amino_n['W']] ); | |
| 2150 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "seq1[i] = %c, cpmx1s[c][3][%d] = %f\n", seq1[0][i], i, cpmx1s[c][3][i] ); | |
| 2151 // for( i=0; i<lgth2; i++ ) fprintf( stderr, "seq2[i] = %c, cpmx2s[c][3][%d] = %f\n", seq2[0][i], i, cpmx2s[c][3][i] ); | |
| 2152 match_calc_add( matrices[c], initverticalw, cpmx2s[c], cpmx1s[c], 0, lgth1, floatwork[c], intwork[c], 1 ); | |
| 2153 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "c=%d, %d - %f\n", c, i, initverticalw[i] ); | |
| 2154 } | |
| 2155 // for( i=0; i<lgth1; i++ ) fprintf( stderr, "%d - %f\n", i, initverticalw[i] ); | |
| 2156 #endif | |
| 2157 | |
| 2158 // exit( 1 ); | |
| 2159 | |
| 2160 if( localhom ) | |
| 2161 imp_match_out_vead_tate( initverticalw, 0, lgth1 ); // 060306 | |
| 2162 | |
| 2163 #if MCD | |
| 2164 match_calc_dynamicmtx( pairoffset, n_dynamicmtx, currentw, icyc, seq1, eff1, jcyc, seq2, eff2, 0, lgth2, *floatwork, *intwork, 1, 0 ); | |
| 2165 // for( i=0; i<lgth2; i++ ) fprintf( stderr, "%d - %f\n", i, currentw[i] ); | |
| 2166 // exit( 1 ); | |
| 2167 #else | |
| 2168 fillzero( currentw, lgth2 ); | |
| 2169 for( c=0; c<maxdistclass; c++ ) | |
| 2170 match_calc_add( matrices[c], currentw, cpmx1s[c], cpmx2s[c], 0, lgth2, floatwork[c], intwork[c], 1 ); | |
| 2171 // for( i=0; i<lgth2; i++ ) fprintf( stderr, "%d - %f\n", i, currentw[i] ); | |
| 2172 // exit( 1 ); | |
| 2173 #endif | |
| 2174 | |
| 2175 if( localhom ) | |
| 2176 imp_match_out_vead( currentw, 0, lgth2 ); // 060306 | |
| 2177 #if 0 // -> tbfast.c | |
| 2178 if( localhom ) | |
| 2179 imp_match_calc( n_dynamicmtx, currentw, icyc, jcyc, lgth1, lgth2, seq1, seq2, eff1, eff2, localhom, 1, 0 ); | |
| 2180 | |
| 2181 #endif | |
| 2182 | |
| 2183 if( headgp == 1 ) | |
| 2184 { | |
| 2185 for( i=1; i<lgth1+1; i++ ) | |
| 2186 { | |
| 2187 // initverticalw[i] += ( ogcp1[0] + fgcp1[i-1] ) ; | |
| 2188 initverticalw[i] += ( ogcp1[0] * headgapfreq2 + fgcp1[i-1] * gapfreq2[0] ) ; | |
| 2189 } | |
| 2190 for( j=1; j<lgth2+1; j++ ) | |
| 2191 { | |
| 2192 // currentw[j] += ( ogcp2[0] + fgcp2[j-1] ) ; | |
| 2193 currentw[j] += ( ogcp2[0] * headgapfreq1 + fgcp2[j-1] * gapfreq1[0] ) ; | |
| 2194 } | |
| 2195 } | |
| 2196 #if OUTGAP0TRY | |
| 2197 else | |
| 2198 { | |
| 2199 fprintf( stderr, "offset = %d\n", offset ); | |
| 2200 for( j=1; j<lgth2+1; j++ ) | |
| 2201 currentw[j] -= offset * j / 2.0; | |
| 2202 for( i=1; i<lgth1+1; i++ ) | |
| 2203 initverticalw[i] -= offset * i / 2.0; | |
| 2204 } | |
| 2205 #endif | |
| 2206 #if 0 | |
| 2207 fprintf( stderr, "\n " ); | |
| 2208 for( j=0; j<lgth2+1; j++ ) fprintf( stderr, " %c ", seq2[0][j] ); | |
| 2209 fprintf( stderr, "\n%c ", seq1[0][0] ); | |
| 2210 for( j=0; j<lgth2+1; j++ ) | |
| 2211 { | |
| 2212 fprintf( stderr, "%5.0f ", currentw[j] ); | |
| 2213 } | |
| 2214 fprintf( stderr, "\n" ); | |
| 2215 #endif | |
| 2216 | |
| 2217 | |
| 2218 for( j=1; j<lgth2+1; ++j ) | |
| 2219 { | |
| 2220 // m[j] = currentw[j-1] + ogcp1[1]; mp[j] = 0; | |
| 2221 m[j] = currentw[j-1] + ogcp1[1] * gapfreq2[j-1]; mp[j] = 0;; | |
| 2222 } | |
| 2223 if( lgth2 == 0 ) | |
| 2224 lastverticalw[0] = 0.0; // Falign kara yobaretatoki kounarukanousei ari | |
| 2225 else | |
| 2226 lastverticalw[0] = currentw[lgth2-1]; | |
| 2227 | |
| 2228 if( tailgp ) lasti = lgth1+1; else lasti = lgth1; | |
| 2229 lastj = lgth2+1; | |
| 2230 | |
| 2231 #if XXXXXXX | |
| 2232 fprintf( stderr, "currentw = \n" ); | |
| 2233 for( i=0; i<lgth1+1; i++ ) | |
| 2234 { | |
| 2235 fprintf( stderr, "%5.2f ", currentw[i] ); | |
| 2236 } | |
| 2237 fprintf( stderr, "\n" ); | |
| 2238 fprintf( stderr, "initverticalw = \n" ); | |
| 2239 for( i=0; i<lgth2+1; i++ ) | |
| 2240 { | |
| 2241 fprintf( stderr, "%5.2f ", initverticalw[i] ); | |
| 2242 } | |
| 2243 fprintf( stderr, "\n" ); | |
| 2244 fprintf( stderr, "fcgp\n" ); | |
| 2245 for( i=0; i<lgth1; i++ ) | |
| 2246 fprintf( stderr, "fgcp1[%d]=%f\n", i, ogcp1[i] ); | |
| 2247 for( i=0; i<lgth2; i++ ) | |
| 2248 fprintf( stderr, "fgcp2[%d]=%f\n", i, ogcp2[i] ); | |
| 2249 #endif | |
| 2250 | |
| 2251 for( i=1; i<lasti; i++ ) | |
| 2252 { | |
| 2253 | |
| 2254 #ifdef enablemultithread | |
| 2255 // fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref ); | |
| 2256 if( chudanpt && *chudanpt != chudanref ) | |
| 2257 { | |
| 2258 // fprintf( stderr, "\n\n## CHUUDAN!!! S\n" ); | |
| 2259 *chudanres = 1; | |
| 2260 return( -1.0 ); | |
| 2261 } | |
| 2262 #endif | |
| 2263 wtmp = previousw; | |
| 2264 previousw = currentw; | |
| 2265 currentw = wtmp; | |
| 2266 | |
| 2267 previousw[0] = initverticalw[i-1]; | |
| 2268 | |
| 2269 #if MCD | |
| 2270 match_calc_dynamicmtx( pairoffset, n_dynamicmtx, currentw, icyc, seq1, eff1, jcyc, seq2, eff2, i, lgth2, *floatwork, *intwork, 0, 0 ); | |
| 2271 #else | |
| 2272 fillzero( currentw, lgth2 ); | |
| 2273 for( c=0; c<maxdistclass; c++ ) | |
| 2274 match_calc_add( matrices[c], currentw, cpmx1s[c], cpmx2s[c], i, lgth2, floatwork[c], intwork[c], 0 ); | |
| 2275 #endif | |
| 2276 #if 0 | |
| 2277 if( i == 1 ) | |
| 2278 { | |
| 2279 fprintf( stderr, "\n" ); | |
| 2280 for( j=0; j<lgth2; j++ ) fprintf( stderr, "%d - %f\n", j, currentw[j] ); | |
| 2281 exit( 1 ); | |
| 2282 } | |
| 2283 #endif | |
| 2284 | |
| 2285 #if XXXXXXX | |
| 2286 fprintf( stderr, "\n" ); | |
| 2287 fprintf( stderr, "i=%d\n", i ); | |
| 2288 fprintf( stderr, "currentw = \n" ); | |
| 2289 for( j=0; j<lgth2; j++ ) | |
| 2290 { | |
| 2291 fprintf( stderr, "%5.2f ", currentw[j] ); | |
| 2292 } | |
| 2293 fprintf( stderr, "\n" ); | |
| 2294 #endif | |
| 2295 if( localhom ) | |
| 2296 { | |
| 2297 // fprintf( stderr, "Calling imp_match_calc (o) lgth = %d, i = %d\n", lgth1, i ); | |
| 2298 #if 0 | |
| 2299 imp_match_out_vead( currentw, i, lgth2 ); | |
| 2300 #else | |
| 2301 imp_match_out_vead( currentw, i, lgth2 ); | |
| 2302 #endif | |
| 2303 } | |
| 2304 #if XXXXXXX | |
| 2305 fprintf( stderr, "\n" ); | |
| 2306 fprintf( stderr, "i=%d\n", i ); | |
| 2307 fprintf( stderr, "currentw = \n" ); | |
| 2308 for( j=0; j<lgth2; j++ ) | |
| 2309 { | |
| 2310 fprintf( stderr, "%5.2f ", currentw[j] ); | |
| 2311 } | |
| 2312 fprintf( stderr, "\n" ); | |
| 2313 #endif | |
| 2314 currentw[0] = initverticalw[i]; | |
| 2315 | |
| 2316 #if 0 | |
| 2317 fprintf( stderr, "%c ", seq1[0][i] ); | |
| 2318 for( j=0; j<lgth2+1; j++ ) | |
| 2319 { | |
| 2320 fprintf( stderr, "%5.0f ", currentw[j] ); | |
| 2321 } | |
| 2322 fprintf( stderr, "\n" ); | |
| 2323 #endif | |
| 2324 | |
| 2325 // mi = previousw[0] + ogcp2[1]; mpi = 0; | |
| 2326 mi = previousw[0] + ogcp2[1] * gapfreq1[i-1]; mpi=0; | |
| 2327 ijppt = ijp[i] + 1; | |
| 2328 mjpt = m + 1; | |
| 2329 prept = previousw; | |
| 2330 curpt = currentw + 1; | |
| 2331 mpjpt = mp + 1; | |
| 2332 fgcp2pt = fgcp2; | |
| 2333 ogcp2pt = ogcp2 + 1; | |
| 2334 fgcp1va = fgcp1[i-1]; | |
| 2335 ogcp1va = ogcp1[i]; | |
| 2336 gf1va = gapfreq1[i]; | |
| 2337 gf1vapre = gapfreq1[i-1]; | |
| 2338 gf2pt = gapfreq2+1; | |
| 2339 gf2ptpre = gapfreq2; | |
| 2340 | |
| 2341 if( trywarp ) | |
| 2342 { | |
| 2343 prevwmrecordspt = prevwmrecords; | |
| 2344 wmrecordspt = wmrecords+1; | |
| 2345 wmrecords1pt = wmrecords; | |
| 2346 warpipt = warpi + 1; | |
| 2347 warpjpt = warpj + 1; | |
| 2348 } | |
| 2349 | |
| 2350 for( j=1; j<lastj; j++ ) | |
| 2351 { | |
| 2352 #ifdef xxxenablemultithread | |
| 2353 // fprintf( stderr, "chudan = %d, %d\n", *chudanpt, chudanref ); | |
| 2354 if( chudanpt && *chudanpt != chudanref ) | |
| 2355 { | |
| 2356 // fprintf( stderr, "\n\n## CHUUDAN!!! S\n" ); | |
| 2357 *chudanres = 1; | |
| 2358 return( -1.0 ); | |
| 2359 } | |
| 2360 #endif | |
| 2361 wm = *prept; | |
| 2362 *ijppt = 0; | |
| 2363 | |
| 2364 #if 0 | |
| 2365 fprintf( stderr, "\n i=%d, j=%d %c, %c", i, j, seq1[0][i], seq2[0][j] ); | |
| 2366 fprintf( stderr, "%5.0f->", wm ); | |
| 2367 fprintf( stderr, "%5.0f? (penal=%5.2f)", g=mi+*fgcp2pt*(1.0-gapfreq1[i]), *fgcp2pt*(1.0-gapfreq1[i]) ); | |
| 2368 #endif | |
| 2369 if( (g=mi+*fgcp2pt*gf1va) > wm ) | |
| 2370 { | |
| 2371 wm = g; | |
| 2372 *ijppt = -( j - mpi ); | |
| 2373 // fprintf( stderr, "Jump to %d (%c)!", mpi, seq2[0][mpi] ); | |
| 2374 } | |
| 2375 if( (g=*prept+*ogcp2pt*gf1vapre) >= mi ) | |
| 2376 { | |
| 2377 mi = g; | |
| 2378 mpi = j-1; | |
| 2379 } | |
| 2380 #if USE_PENALTY_EX | |
| 2381 mi += fpenalty_ex; | |
| 2382 #endif | |
| 2383 | |
| 2384 #if 0 | |
| 2385 fprintf( stderr, "%5.0f->", wm ); | |
| 2386 fprintf( stderr, "%5.0f? (penal=%5.2f)", g=*mjpt+fgcp1va*(1.0-gapfreq2[j]), fgcp1va*(1.0-gapfreq2[j]) ); | |
| 2387 #endif | |
| 2388 if( (g=*mjpt+ fgcp1va* *gf2pt) > wm ) | |
| 2389 { | |
| 2390 wm = g; | |
| 2391 *ijppt = +( i - *mpjpt ); | |
| 2392 // fprintf( stderr, "Jump to %d (%c)!", *mpjpt, seq1[0][*mpjpt] ); | |
| 2393 } | |
| 2394 if( (g=*prept+ ogcp1va* *gf2ptpre) >= *mjpt ) | |
| 2395 { | |
| 2396 *mjpt = g; | |
| 2397 *mpjpt = i-1; | |
| 2398 } | |
| 2399 #if USE_PENALTY_EX | |
| 2400 m[j] += fpenalty_ex; | |
| 2401 #endif | |
| 2402 if( trywarp ) | |
| 2403 { | |
| 2404 #if USE_PENALTY_EX | |
| 2405 if( ( g=*prevwmrecordspt++ + fpenalty_shift + fpenalty_ex * ( i - prevwarpi[j-1] + j - prevwarpj[j-1] ) ) > wm ) // naka ha osokute kamawanai | |
| 2406 #else | |
| 2407 if( ( g=*prevwmrecordspt++ + fpenalty_shift ) > wm ) // naka ha osokute kamawanai | |
| 2408 #endif | |
| 2409 { | |
| 2410 if( warpn && prevwarpi[j-1] == warpis[warpn-1] && prevwarpj[j-1] == warpjs[warpn-1] ) | |
| 2411 { | |
| 2412 *ijppt = warpbase + warpn - 1; | |
| 2413 } | |
| 2414 else | |
| 2415 { | |
| 2416 *ijppt = warpbase + warpn; | |
| 2417 warpis = realloc( warpis, sizeof(int) * ( warpn+1 ) ); | |
| 2418 warpjs = realloc( warpjs, sizeof(int) * ( warpn+1 ) ); | |
| 2419 warpis[warpn] = prevwarpi[j-1]; | |
| 2420 warpjs[warpn] = prevwarpj[j-1]; | |
| 2421 warpn++; | |
| 2422 } | |
| 2423 wm = g; | |
| 2424 } | |
| 2425 curm = *curpt + wm; | |
| 2426 | |
| 2427 if( *wmrecords1pt > *wmrecordspt ) | |
| 2428 { | |
| 2429 *wmrecordspt = *wmrecords1pt; | |
| 2430 *warpipt = *(warpipt-1); | |
| 2431 *warpjpt = *(warpjpt-1); | |
| 2432 } | |
| 2433 if( curm > *wmrecordspt ) | |
| 2434 { | |
| 2435 *wmrecordspt = curm; | |
| 2436 *warpipt = i; | |
| 2437 *warpjpt = j; | |
| 2438 } | |
| 2439 wmrecordspt++; | |
| 2440 wmrecords1pt++; | |
| 2441 warpipt++; | |
| 2442 warpjpt++; | |
| 2443 } | |
| 2444 | |
| 2445 #if 0 | |
| 2446 fprintf( stderr, "%5.0f ", wm ); | |
| 2447 #endif | |
| 2448 *curpt++ += wm; | |
| 2449 ijppt++; | |
| 2450 mjpt++; | |
| 2451 prept++; | |
| 2452 mpjpt++; | |
| 2453 fgcp2pt++; | |
| 2454 ogcp2pt++; | |
| 2455 gf2ptpre++; | |
| 2456 gf2pt++; | |
| 2457 | |
| 2458 } | |
| 2459 lastverticalw[i] = currentw[lgth2-1]; | |
| 2460 | |
| 2461 if( trywarp ) | |
| 2462 { | |
| 2463 fltncpy( prevwmrecords, wmrecords, lastj ); | |
| 2464 intncpy( prevwarpi, warpi, lastj ); | |
| 2465 intncpy( prevwarpj, warpj, lastj ); | |
| 2466 } | |
| 2467 } | |
| 2468 if( trywarp ) | |
| 2469 { | |
| 2470 // fprintf( stderr, "wm = %f\n", wm ); | |
| 2471 // fprintf( stderr, "warpn = %d\n", warpn ); | |
| 2472 free( wmrecords ); | |
| 2473 free( prevwmrecords ); | |
| 2474 free( warpi ); | |
| 2475 free( warpj ); | |
| 2476 free( prevwarpi ); | |
| 2477 free( prevwarpj ); | |
| 2478 } | |
| 2479 | |
| 2480 | |
| 2481 #if OUTGAP0TRY | |
| 2482 if( !outgap ) | |
| 2483 { | |
| 2484 for( j=1; j<lgth2+1; j++ ) | |
| 2485 currentw[j] -= offset * ( lgth2 - j ) / 2.0; | |
| 2486 for( i=1; i<lgth1+1; i++ ) | |
| 2487 lastverticalw[i] -= offset * ( lgth1 - i / 2.0); | |
| 2488 } | |
| 2489 #endif | |
| 2490 | |
| 2491 /* | |
| 2492 fprintf( stderr, "\n" ); | |
| 2493 for( i=0; i<icyc; i++ ) fprintf( stderr,"%s\n", seq1[i] ); | |
| 2494 fprintf( stderr, "#####\n" ); | |
| 2495 for( j=0; j<jcyc; j++ ) fprintf( stderr,"%s\n", seq2[j] ); | |
| 2496 fprintf( stderr, "====>" ); | |
| 2497 for( i=0; i<icyc; i++ ) strcpy( mseq1[i], seq1[i] ); | |
| 2498 for( j=0; j<jcyc; j++ ) strcpy( mseq2[j], seq2[j] ); | |
| 2499 */ | |
| 2500 if( localhom ) | |
| 2501 { | |
| 2502 Atracking_localhom( impmatch, currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, warpis, warpjs, warpbase ); | |
| 2503 } | |
| 2504 else | |
| 2505 Atracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, ijp, icyc, jcyc, tailgp, warpis, warpjs, warpbase ); | |
| 2506 | |
| 2507 if( warpis ) free( warpis ); | |
| 2508 if( warpjs ) free( warpjs ); | |
| 2509 | |
| 2510 // fprintf( stderr, "### impmatch = %f\n", *impmatch ); | |
| 2511 | |
| 2512 resultlen = strlen( mseq1[0] ); | |
| 2513 if( alloclen < resultlen || resultlen > N ) | |
| 2514 { | |
| 2515 fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N ); | |
| 2516 ErrorExit( "LENGTH OVER!\n" ); | |
| 2517 } | |
| 2518 | |
| 2519 | |
| 2520 for( i=0; i<icyc; i++ ) strcpy( seq1[i], mseq1[i] ); | |
| 2521 for( j=0; j<jcyc; j++ ) strcpy( seq2[j], mseq2[j] ); | |
| 2522 #if 0 | |
| 2523 fprintf( stderr, "\n" ); | |
| 2524 for( i=0; i<icyc; i++ ) fprintf( stderr, "%s\n", mseq1[i] ); | |
| 2525 fprintf( stderr, "#####\n" ); | |
| 2526 for( j=0; j<jcyc; j++ ) fprintf( stderr, "%s\n", mseq2[j] ); | |
| 2527 #endif | |
| 2528 | |
| 2529 // fprintf( stderr, "wm = %f\n", wm ); | |
| 2530 | |
| 2531 return( wm ); | |
| 2532 } |