Context Navigation

genalign11.c

Visit:

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

Line
1	#include "mltaln.h"
2	#include "dp.h"
3
4	#define DEBUG 0
5	#define DEBUG2 0
6	#define XXXXXXX 0
7	#define USE_PENALTY_EX 1
8
9	static TLS int localstop;
10
11	#if 1
12	static void match_calc( float match, char s1, char *s2, int i1, int lgth2 )
13	{
14	char tmpc = s1[0][i1];
15	char *seq2 = s2[0];
16
17	while( lgth2-- )
18	match++ = amino_dis[(int)tmpc][(int)seq2++];
19	}
20	#else
21	static void match_calc( float match, char s1, char *s2, int i1, int lgth2 )
22	{
23	int j;
24
25	for( j=0; j<lgth2; j++ )
26	match[j] = amino_dis[(s1)[i1]][(s2)[j]];
27	}
28	#endif
29
30	#if 0
31	static void match_calc_bk( float match, float cpmx1, float cpmx2, int i1, int lgth2, float floatwork, int *intwork, int initialize )
32	{
33	int j, k, l;
34	float scarr[26];
35	float **cpmxpd = floatwork;
36	int **cpmxpdn = intwork;
37	int count = 0;
38
39	if( initialize )
40	{
41	for( j=0; j<lgth2; j++ )
42	{
43	count = 0;
44	for( l=0; l<26; l++ )
45	{
46	if( cpmx2[l][j] )
47	{
48	cpmxpd[count][j] = cpmx2[l][j];
49	cpmxpdn[count][j] = l;
50	count++;
51	}
52	}
53	cpmxpdn[count][j] = -1;
54	}
55	}
56
57	for( l=0; l<26; l++ )
58	{
59	scarr[l] = 0.0;
60	for( k=0; k<26; k++ )
61	scarr[l] += n_dis[k][l] * cpmx1[k][i1];
62	}
63	#if 0
64	{
65	float fpt, fptpt, fpt2;
66	int ipt, *iptpt;
67	fpt2 = match;
68	iptpt = cpmxpdn;
69	fptpt = cpmxpd;
70	while( lgth2-- )
71	{
72	*fpt2 = 0.0;
73	ipt=iptpt,fpt=fptpt;
74	while( *ipt > -1 )
75	fpt2 += scarr[ipt++] * *fpt++;
76	fpt2++,iptpt++,fptpt++;
77	}
78	}
79	#else
80	for( j=0; j<lgth2; j++ )
81	{
82	match[j] = 0.0;
83	for( k=0; cpmxpdn[k][j]>-1; k++ )
84	match[j] += scarr[cpmxpdn[k][j]] * cpmxpd[k][j];
85	}
86	#endif
87	}
88	#endif
89
90	static float gentracking( float lasthorizontalw, float lastverticalw,
91	char seq1, char seq2,
92	char mseq1, char mseq2,
93	float cpmx1, float cpmx2,
94	int ijpi, int ijpj, int off1pt, int off2pt, int endi, int endj )
95	{
96	int i, j, l, iin, jin, lgth1, lgth2, k, limk;
97	int ifi=0, jfi=0; // by D.Mathog
98	// char gap[] = "-";
99	char *gap;
100	gap = newgapstr;
101	lgth1 = strlen( seq1[0] );
102	lgth2 = strlen( seq2[0] );
103
104	#if 0
105	for( i=0; i<lgth1; i++ )
106	{
107	fprintf( stderr, "lastverticalw[%d] = %f\n", i, lastverticalw[i] );
108	}
109	#endif
110
111	for( i=0; i<lgth1+1; i++ )
112	{
113	ijpi[i][0] = localstop;
114	ijpj[i][0] = localstop;
115	}
116	for( j=0; j<lgth2+1; j++ )
117	{
118	ijpi[0][j] = localstop;
119	ijpj[0][j] = localstop;
120	}
121
122	mseq1[0] += lgth1+lgth2;
123	*mseq1[0] = 0;
124	mseq2[0] += lgth1+lgth2;
125	*mseq2[0] = 0;
126	iin = endi; jin = endj;
127	limk = lgth1+lgth2;
128	for( k=0; k<=limk; k++ )
129	{
130
131	ifi = ( ijpi[iin][jin] );
132	jfi = ( ijpj[iin][jin] );
133	l = iin - ifi;
134	// if( ijpi[iin][jin] < 0 \|\| ijpj[iin][jin] < 0 )
135	// {
136	// fprintf( stderr, "skip! %d-%d\n", ijpi[iin][jin], ijpj[iin][jin] );
137	// fprintf( stderr, "1: %c-%c\n", seq1[0][iin], seq1[0][ifi] );
138	// fprintf( stderr, "2: %c-%c\n", seq2[0][jin], seq2[0][jfi] );
139	// }
140	while( --l )
141	{
142	*--mseq1[0] = seq1[0][ifi+l];
143	--mseq2[0] = gap;
144	k++;
145	}
146	l= jin - jfi;
147	while( --l )
148	{
149	--mseq1[0] = gap;
150	*--mseq2[0] = seq2[0][jfi+l];
151	k++;
152	}
153
154	if( iin <= 0 \|\| jin <= 0 ) break;
155	*--mseq1[0] = seq1[0][ifi];
156	*--mseq2[0] = seq2[0][jfi];
157
158	if( ijpi[ifi][jfi] == localstop ) break;
159	if( ijpj[ifi][jfi] == localstop ) break;
160	k++;
161	iin = ifi; jin = jfi;
162	}
163	if( ifi == -1 ) off1pt = 0; else off1pt = ifi;
164	if( jfi == -1 ) off2pt = 0; else off2pt = jfi;
165
166	// fprintf( stderr, "ifn = %d, jfn = %d\n", ifi, jfi );
167
168
169	return( 0.0 );
170	}
171
172
173	float genL__align11( char seq1, char seq2, int alloclen, int off1pt, int off2pt )
174	/* score no keisan no sai motokaraaru gap no atukai ni mondai ga aru */
175	{
176	// int k;
177	register int i, j;
178	int lasti, lastj;
179	int lgth1, lgth2;
180	int resultlen;
181	float wm = 0.0; /* int ?????? */
182	float g;
183	float currentw, previousw;
184	#if 1
185	float *wtmp;
186	int *ijpipt;
187	int *ijpjpt;
188	float mjpt, Mjpt, prept, curpt;
189	int mpjpt, Mpjpt;
190	#endif
191	static TLS float mi, *m;
192	static TLS float Mi, *largeM;
193	static TLS int **ijpi;
194	static TLS int **ijpj;
195	static TLS int mpi, *mp;
196	static TLS int Mpi, *Mp;
197	static TLS float w1, w2;
198	static TLS float *match;
199	static TLS float initverticalw; / kufuu sureba iranai */
200	static TLS float lastverticalw; / kufuu sureba iranai */
201	static TLS char **mseq1;
202	static TLS char **mseq2;
203	static TLS char **mseq;
204	static TLS float **cpmx1;
205	static TLS float **cpmx2;
206	static TLS int **intwork;
207	static TLS float **floatwork;
208	static TLS int orlgth1 = 0, orlgth2 = 0;
209	float maxwm;
210	float tbk;
211	int tbki, tbkj;
212	int endali, endalj;
213	// float localthr = 0.0;
214	// float localthr2 = 0.0;
215	float fpenalty = (float)penalty;
216	float fpenalty_OP = (float)penalty_OP;
217	float fpenalty_ex = (float)penalty_ex;
218	// float fpenalty_EX = (float)penalty_EX;
219	float foffset = (float)offset;
220	float localthr = -foffset;
221	float localthr2 = -foffset;
222
223	if( seq1 == NULL )
224	{
225	if( orlgth1 > 0 && orlgth2 > 0 )
226	{
227	orlgth1 = 0;
228	orlgth2 = 0;
229	free( mseq1 );
230	free( mseq2 );
231	FreeFloatVec( w1 );
232	FreeFloatVec( w2 );
233	FreeFloatVec( match );
234	FreeFloatVec( initverticalw );
235	FreeFloatVec( lastverticalw );
236
237	FreeFloatVec( m );
238	FreeIntVec( mp );
239	free( largeM );
240	free( Mp );
241
242	FreeCharMtx( mseq );
243
244	FreeFloatMtx( cpmx1 );
245	FreeFloatMtx( cpmx2 );
246
247	FreeFloatMtx( floatwork );
248	FreeIntMtx( intwork );
249
250	}
251	return( 0.0 );
252	}
253
254
255
256	// fprintf( stderr, "@@@@@@@@@@@@@ penalty_OP = %f, penalty_EX = %f, pelanty = %f\n", fpenalty_OP, fpenalty_EX, fpenalty );
257
258	if( orlgth1 == 0 )
259	{
260	mseq1 = AllocateCharMtx( njob, 0 );
261	mseq2 = AllocateCharMtx( njob, 0 );
262	}
263
264
265	lgth1 = strlen( seq1[0] );
266	lgth2 = strlen( seq2[0] );
267
268	if( lgth1 > orlgth1 \|\| lgth2 > orlgth2 )
269	{
270	int ll1, ll2;
271
272	if( orlgth1 > 0 && orlgth2 > 0 )
273	{
274	FreeFloatVec( w1 );
275	FreeFloatVec( w2 );
276	FreeFloatVec( match );
277	FreeFloatVec( initverticalw );
278	FreeFloatVec( lastverticalw );
279
280	FreeFloatVec( m );
281	FreeIntVec( mp );
282	FreeFloatVec( largeM );
283	FreeIntVec( Mp );
284
285	FreeCharMtx( mseq );
286
287	FreeFloatMtx( cpmx1 );
288	FreeFloatMtx( cpmx2 );
289
290	FreeFloatMtx( floatwork );
291	FreeIntMtx( intwork );
292	}
293
294	ll1 = MAX( (int)(1.3*lgth1), orlgth1 ) + 100;
295	ll2 = MAX( (int)(1.3*lgth2), orlgth2 ) + 100;
296
297	#if DEBUG
298	fprintf( stderr, "\ntrying to allocate (%d+%d)xn matrices ... ", ll1, ll2 );
299	#endif
300
301	w1 = AllocateFloatVec( ll2+2 );
302	w2 = AllocateFloatVec( ll2+2 );
303	match = AllocateFloatVec( ll2+2 );
304
305	initverticalw = AllocateFloatVec( ll1+2 );
306	lastverticalw = AllocateFloatVec( ll1+2 );
307
308	m = AllocateFloatVec( ll2+2 );
309	mp = AllocateIntVec( ll2+2 );
310	largeM = AllocateFloatVec( ll2+2 );
311	Mp = AllocateIntVec( ll2+2 );
312
313	mseq = AllocateCharMtx( njob, ll1+ll2 );
314
315	cpmx1 = AllocateFloatMtx( 26, ll1+2 );
316	cpmx2 = AllocateFloatMtx( 26, ll2+2 );
317
318	floatwork = AllocateFloatMtx( 26, MAX( ll1, ll2 )+2 );
319	intwork = AllocateIntMtx( 26, MAX( ll1, ll2 )+2 );
320
321	#if DEBUG
322	fprintf( stderr, "succeeded\n" );
323	#endif
324
325	orlgth1 = ll1 - 100;
326	orlgth2 = ll2 - 100;
327	}
328
329
330	mseq1[0] = mseq[0];
331	mseq2[0] = mseq[1];
332
333
334	if( orlgth1 > commonAlloc1 \|\| orlgth2 > commonAlloc2 )
335	{
336	int ll1, ll2;
337
338	if( commonAlloc1 && commonAlloc2 )
339	{
340	FreeIntMtx( commonIP );
341	FreeIntMtx( commonJP );
342	}
343
344	ll1 = MAX( orlgth1, commonAlloc1 );
345	ll2 = MAX( orlgth2, commonAlloc2 );
346
347	#if DEBUG
348	fprintf( stderr, "\n\ntrying to allocate %dx%d matrices ... ", ll1+1, ll2+1 );
349	#endif
350
351	commonIP = AllocateIntMtx( ll1+10, ll2+10 );
352	commonJP = AllocateIntMtx( ll1+10, ll2+10 );
353
354	#if DEBUG
355	fprintf( stderr, "succeeded\n\n" );
356	#endif
357
358	commonAlloc1 = ll1;
359	commonAlloc2 = ll2;
360	}
361	ijpi = commonIP;
362	ijpj = commonJP;
363
364
365	#if 0
366	for( i=0; i<lgth1; i++ )
367	fprintf( stderr, "ogcp1[%d]=%f\n", i, ogcp1[i] );
368	#endif
369
370	currentw = w1;
371	previousw = w2;
372
373	match_calc( initverticalw, seq2, seq1, 0, lgth1 );
374
375	match_calc( currentw, seq1, seq2, 0, lgth2 );
376
377
378	lasti = lgth2+1;
379	for( j=1; j<lasti; ++j )
380	{
381	m[j] = currentw[j-1]; mp[j] = 0;
382	largeM[j] = currentw[j-1]; Mp[j] = 0;
383	}
384
385	lastverticalw[0] = currentw[lgth2-1];
386
387
388	#if 0
389	fprintf( stderr, "currentw = \n" );
390	for( i=0; i<lgth1+1; i++ )
391	{
392	fprintf( stderr, "%5.2f ", currentw[i] );
393	}
394	fprintf( stderr, "\n" );
395	fprintf( stderr, "initverticalw = \n" );
396	for( i=0; i<lgth2+1; i++ )
397	{
398	fprintf( stderr, "%5.2f ", initverticalw[i] );
399	}
400	fprintf( stderr, "\n" );
401	#endif
402	#if DEBUG2
403	fprintf( stderr, "\n" );
404	fprintf( stderr, " " );
405	for( j=0; j<lgth2+1; j++ )
406	fprintf( stderr, "%c ", seq2[0][j] );
407	fprintf( stderr, "\n" );
408	#endif
409
410	localstop = lgth1+lgth2+1;
411	maxwm = -999999999.9;
412	endali = endalj = 0;
413	#if DEBUG2
414	fprintf( stderr, "\n" );
415	fprintf( stderr, "%c ", seq1[0][0] );
416
417	for( j=0; j<lgth2+1; j++ )
418	fprintf( stderr, "%5.0f ", currentw[j] );
419	fprintf( stderr, "\n" );
420	#endif
421
422	lasti = lgth1+1;
423	for( i=1; i<lasti; i++ )
424	{
425	wtmp = previousw;
426	previousw = currentw;
427	currentw = wtmp;
428
429	previousw[0] = initverticalw[i-1];
430
431	match_calc( currentw, seq1, seq2, i, lgth2 );
432	#if DEBUG2
433	fprintf( stderr, "%c ", seq1[0][i] );
434	fprintf( stderr, "%5.0f ", currentw[0] );
435	#endif
436
437	#if XXXXXXX
438	fprintf( stderr, "\n" );
439	fprintf( stderr, "i=%d\n", i );
440	fprintf( stderr, "currentw = \n" );
441	for( j=0; j<lgth2; j++ )
442	{
443	fprintf( stderr, "%5.2f ", currentw[j] );
444	}
445	fprintf( stderr, "\n" );
446	#endif
447	#if XXXXXXX
448	fprintf( stderr, "\n" );
449	fprintf( stderr, "i=%d\n", i );
450	fprintf( stderr, "currentw = \n" );
451	for( j=0; j<lgth2; j++ )
452	{
453	fprintf( stderr, "%5.2f ", currentw[j] );
454	}
455	fprintf( stderr, "\n" );
456	#endif
457	currentw[0] = initverticalw[i];
458
459	mi = previousw[0]; mpi = 0;
460	Mi = previousw[0]; Mpi = 0;
461
462	#if 0
463	if( mi < localthr ) mi = localthr2;
464	#endif
465
466	ijpipt = ijpi[i] + 1;
467	ijpjpt = ijpj[i] + 1;
468	mjpt = m + 1;
469	Mjpt = largeM + 1;
470	prept = previousw;
471	curpt = currentw + 1;
472	mpjpt = mp + 1;
473	Mpjpt = Mp + 1;
474	tbk = -999999.9;
475	tbki = 0;
476	tbkj = 0;
477	lastj = lgth2+1;
478	for( j=1; j<lastj; j++ )
479	{
480	wm = *prept;
481	*ijpipt = i-1;
482	*ijpjpt = j-1;
483
484
485	// fprintf( stderr, "i,j=%d,%d %c-%c\n", i, j, seq1[0][i], seq2[0][j] );
486	// fprintf( stderr, "wm=%f\n", wm );
487	#if 0
488	fprintf( stderr, "%5.0f->", wm );
489	#endif
490	g = mi + fpenalty;
491	#if 0
492	fprintf( stderr, "%5.0f?", g );
493	#endif
494	if( g > wm )
495	{
496	wm = g;
497	// *ijpipt = i - 1;
498	*ijpjpt = mpi;
499	}
500	g = *prept;
501	if( g > mi )
502	{
503	mi = g;
504	mpi = j-1;
505	}
506
507	#if USE_PENALTY_EX
508	mi += fpenalty_ex;
509	#endif
510
511	#if 0
512	fprintf( stderr, "%5.0f->", wm );
513	#endif
514	g = *mjpt + fpenalty;
515	#if 0
516	fprintf( stderr, "m%5.0f?", g );
517	#endif
518	if( g > wm )
519	{
520	wm = g;
521	ijpipt = mpjpt;
522	*ijpjpt = j - 1; //IRU!
523	}
524	g = *prept;
525	if( g > *mjpt )
526	{
527	*mjpt = g;
528	*mpjpt = i-1;
529	}
530	#if USE_PENALTY_EX
531	*mjpt += fpenalty_ex;
532	#endif
533
534
535	g = tbk + fpenalty_OP;
536	// g = tbk;
537	if( g > wm )
538	{
539	wm = g;
540	*ijpipt = tbki;
541	*ijpjpt = tbkj;
542	// fprintf( stderr, "hit! i%d, j%d, ijpi = %d, ijpj = %d\n", i, j, ijpipt, ijpjpt );
543	}
544	// g = Mi;
545	if( Mi > tbk )
546	{
547	tbk = Mi; //error desu.
548	tbki = i-1;
549	tbkj = Mpi;
550	}
551	// g = *Mjpt;
552	if( *Mjpt > tbk )
553	{
554	tbk = *Mjpt;
555	tbki = *Mpjpt;
556	tbkj = j-1;
557	}
558	// tbk += fpenalty_EX;// + foffset;
559
560	// g = *prept;
561	if( prept > Mjpt )
562	{
563	Mjpt = prept;
564	*Mpjpt = i-1;
565	}
566	// *Mjpt += fpenalty_EX;// + foffset;
567
568	// g = *prept;
569	if( *prept > Mi )
570	{
571	Mi = *prept;
572	Mpi = j-1;
573	}
574	// Mi += fpenalty_EX;// + foffset;
575
576
577	// fprintf( stderr, "wm=%f, tbk=%f(%c-%c), mi=%f, mjpt=%f\n", wm, tbk, seq1[0][tbki], seq2[0][tbkj], mi, mjpt );
578	// fprintf( stderr, "ijp = %c,%c\n", seq1[0][abs(ijpipt)], seq2[0][abs(ijpjpt)] );
579
580
581	if( maxwm < wm )
582	{
583	maxwm = wm;
584	endali = i;
585	endalj = j;
586	}
587	#if 1
588	if( wm < localthr )
589	{
590	// fprintf( stderr, "stop i=%d, j=%d, curpt=%f\n", i, j, *curpt );
591	*ijpipt = localstop;
592	// *ijpjpt = localstop;
593	wm = localthr2;
594	}
595	#endif
596	#if 0
597	fprintf( stderr, "%5.0f ", *curpt );
598	#endif
599	#if DEBUG2
600	fprintf( stderr, "%5.0f ", wm );
601	// fprintf( stderr, "%c-%c ijppt = %d, localstop = %d\n", seq1[0][i], seq2[0][j], ijppt, localstop );
602	#endif
603
604	*curpt += wm;
605	ijpipt++;
606	ijpjpt++;
607	mjpt++;
608	Mjpt++;
609	prept++;
610	mpjpt++;
611	Mpjpt++;
612	curpt++;
613	}
614	#if DEBUG2
615	fprintf( stderr, "\n" );
616	#endif
617
618	lastverticalw[i] = currentw[lgth2-1];
619	}
620
621
622	#if DEBUG2
623	fprintf( stderr, "maxwm = %f\n", maxwm );
624	fprintf( stderr, "endali = %d\n", endali );
625	fprintf( stderr, "endalj = %d\n", endalj );
626	#endif
627
628	if( ijpi[endali][endalj] == localstop ) // && ijpj[endali][endalj] == localstop )
629	{
630	strcpy( seq1[0], "" );
631	strcpy( seq2[0], "" );
632	off1pt = off2pt = 0;
633	return( 0.0 );
634	}
635
636
637	gentracking( currentw, lastverticalw, seq1, seq2, mseq1, mseq2, cpmx1, cpmx2, ijpi, ijpj, off1pt, off2pt, endali, endalj );
638
639	// fprintf( stderr, "### impmatch = %f\n", *impmatch );
640
641	resultlen = strlen( mseq1[0] );
642	if( alloclen < resultlen \|\| resultlen > N )
643	{
644	fprintf( stderr, "alloclen=%d, resultlen=%d, N=%d\n", alloclen, resultlen, N );
645	ErrorExit( "LENGTH OVER!\n" );
646	}
647
648
649	strcpy( seq1[0], mseq1[0] );
650	strcpy( seq2[0], mseq2[0] );
651
652	#if 0
653	fprintf( stderr, "\n" );
654	fprintf( stderr, ">\n%s\n", mseq1[0] );
655	fprintf( stderr, ">\n%s\n", mseq2[0] );
656	#endif
657
658
659	return( maxwm );
660	}
661

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

Context Navigation

source: branches/lib/GDE/MAFFT/mafft-7.055-with-extensions/core/genalign11.c

Download in other formats: