Context Navigation

Qalignmm.c

Visit:

Last change on this file was 10371, checked in by aboeckma, 11 years ago
updated mafft version. Added extensions (no svn ignore, yet)
File size: 54.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: tags/arb-6.0/GDE/MAFFT/mafft-7.055-with-extensions/core/Qalignmm.c

Download in other formats: