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