Context Navigation

iteration.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: 9.1 KB

Line
1	/* iteration ( algorithm C ) */
2	#include "mltaln.h"
3
4	#define DEBUG 0
5
6	static void Writeoptions( FILE *fp )
7	{
8	if( scoremtx == 1 )
9	fprintf( fp, "Dayhoff( ... )\n" );
10	else if( scoremtx == -1 )
11	fprintf( fp, "DNA\n" );
12	else if( scoremtx == 2 )
13	fprintf( fp, "Miyata-Yasunaga\n" );
14	else
15	fprintf( fp, "JTT %dPAM\n", pamN );
16
17	if( scoremtx == 0 )
18	fprintf( fp, "Gap Penalty = %+d, %+d\n", penalty, offset );
19	else
20	fprintf( fp, "Gap Penalty = %+d\n", penalty );
21
22	fprintf( fp, "marginal score to search : best - %f\n", cut );
23	if( scmtd == 3 )
24	fprintf( fp, "score of rnd or sco\n" );
25	else if( scmtd == 4 )
26	fprintf( fp, "score = sigma( score for a pair of homologous amino acids ) / ( number of amino acids pairs )\n" );
27	else if( scmtd == 5 )
28	fprintf( fp, "score : SP\n" );
29	if( mix )
30	fprintf( fp, "?\n" );
31	else
32	{
33	if( weight == 2 )
34	fprintf( fp, "weighted, geta2 = %f\n", geta2 );
35	else if( weight == 3 )
36	{
37	if( scmtd == 4 )
38	fprintf( fp, "reversely weighted in function 'align', unweighted in function 'score_calc'\n" );
39	else
40	fprintf( fp, "weighted like ClustalW," );
41	}
42	else
43	fprintf( fp, "unweighted\n" );
44	}
45	if( weight && utree )
46	{
47	fprintf( fp, "using tree defined by the file hat2 with simplified UPG method\n" );
48	}
49	if( weight && !utree )
50	fprintf( fp, "using temporary tree by simplified UPG method\n" );
51	fprintf( fp, "Algorithm %c\n", alg );
52	}
53
54
55
56
57	char *align0( float wm, char *aseq, char seq, double effarr[M], int icyc, int ex )
58	{
59	char **result;
60
61	if( alg == 'B' )
62	{
63	ErrorExit( "Sorry!" );
64	/*
65	if( outgap == 0 )
66	{
67	result = alignm1_o( wm, aseq, seq, scmx, effarr, icyc, ex );
68	}
69	if( outgap == 1 )
70	{
71	result = alignm1( wm, aseq, seq, scmx, effarr, icyc, ex );
72	}
73	*/
74	}
75	else if( alg == 'C' )
76	{
77	result = Calignm1( wm, aseq, seq, effarr, icyc, ex );
78	}
79	return( result );
80	}
81
82
83	double score_m_1_0( char aseq, int locnjob, int s, double eff, double effarr[M] )
84	{
85	double x;
86
87	if( alg == 'B' )
88	{
89	ErrorExit( "Sorry!" );
90	}
91	if( alg == 'C' )
92	{
93	x = Cscore_m_1( aseq, locnjob, s, eff );
94	}
95	fprintf( stderr, "in score_m_1_0 %f\n", x );
96	return( x );
97	}
98
99	int iteration( int locnjob, char name[M][B], int nlen[M], char aseq, char bseq, int *topol, double len, double **eff )
100	{
101	double tscore, mscore;
102	int identity;
103	static char mseq1, *mseq2 = NULL;
104	static char **result;
105	int i, l;
106	static double effarr[M];
107	int s;
108	int sss[2];
109	char ou;
110	int alloclen;
111	int resultlen;
112	int nlenmax0 = nlenmax;
113	FILE *prep;
114	char sai[M];
115	char sai1[M];
116	char sai2[M];
117	#if 0
118	double his[2][M][MAXITERATION/locnjob+1];
119	#else
120	double ***his;
121	#endif
122	int cyc[2];
123	char shindou = 0;
124	float wm;
125	int returnvalue;
126
127	for( i=0; i<locnjob; i++ )
128	{
129	sai[i] = 1;
130	sai1[i] = 1;
131	sai2[i] = 2;
132	}
133	sai[locnjob] = sai1[locnjob] = sai2[locnjob] = 0;
134
135
136	Writeoptions( trap_g );
137
138	his = AllocateDoubleCub( 2, M, MAXITERATION/locnjob+1 );
139
140	if( mseq2 == NULL )
141	{
142	alloclen = nlenmax * 2.0;
143	AllocateTmpSeqs( &mseq2, &mseq1, alloclen );
144	}
145
146	if( !tbitr && !tbweight )
147	{
148	writePre( locnjob, name, nlen, aseq, 0 );
149
150	#if 0
151	prep = fopen( "best", "w" );
152	Write( prep, locnjob, name, nlen, aseq );
153	fclose( prep );
154	#endif
155	}
156
157
158
159
160	treeconstruction( aseq, locnjob, topol, len, eff );
161	tscore = score_calc0( aseq, locnjob, eff, 0 );
162
163	#if DEBUG
164	fprintf( stderr, "eff mtx in iteration\n" );
165	for( i=0; i<locnjob; i++ )
166	{
167	for( j=0; j<locnjob; j++ )
168	{
169	fprintf( stderr, "%5.3f ", eff[i][j] );
170	}
171	fprintf( stderr, "\n" );
172	}
173	#endif
174
175	fprintf( stderr, "\n" );
176	if( disp )
177	{
178	fprintf( stderr, "aseq(initial)\n" );
179	display( aseq, locnjob );
180	}
181	fprintf( stderr, "initial score = %f \n", tscore );
182	fprintf( stderr, "\n" );
183	for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
184	mscore = tscore;
185	srand( time(NULL) );
186
187	sss[1] = 0;
188	sss[0] = locnjob-1;
189	/*
190	sss[0] = (int)( (float)locnjob/2.0 );
191	*/
192	ou = 1;
193	cyc[0] = 0; cyc[1] = 0;
194
195	for( s=-1, l=0; l<MAXITERATION; l++ )
196	{
197	int ss;
198	double tmpscore, tmpscore1;
199
200	if( strlen( aseq[0] ) > nlenmax )
201	nlenmax = strlen( aseq[0] );
202
203	/*
204	s = ( int )( rnd() * locnjob );
205	s++;
206	if( s == locnjob ) s = 0;
207	ou = 0;
208	*/
209	if( ou == 0 )
210	{
211	ou = 1;
212	s = sss[0];
213	/*
214	sss[0]++;
215	if( sss[0] == locnjob )
216	{
217	sss[0] = 0;
218	cyc[0]++;
219	}
220	*/
221	sss[0]--;
222	if( sss[0] == -1 )
223	{
224	sss[0] = locnjob-1;
225	cyc[0]++;
226	}
227	}
228	else
229	{
230	ou = 0;
231	s = sss[1];
232	sss[1]++;
233	if( sss[1] == locnjob )
234	{
235	sss[1] = 0;
236	cyc[1]++;
237	}
238	}
239	fprintf( trap_g, "%d ", weight );
240
241	/*
242	for( i=0, count=0; i<strlen( aseq[s] ); i++ )
243	{
244	if( aseq[s][i] != '-' )
245	{
246	mseq1[count] = aseq[s][i];
247	count++;
248	}
249	}
250	mseq1[count] = 0;
251	*/
252	gappick0( mseq1, aseq[s] );
253
254	if( checkC )
255	tmpscore = score_m_1_0( aseq, locnjob, s, eff, effarr );
256
257	gappick( locnjob, s, aseq, mseq2, eff, effarr );
258
259	result = align0( &wm, mseq2, mseq1, effarr, locnjob-2, s );
260	resultlen = strlen( result[0] );
261	if( resultlen > alloclen )
262	{
263	if( resultlen > nlenmax0*3 \|\| resultlen > N )
264	{
265	fprintf(stderr, "Error in main1\n");
266	exit( 1 );
267	}
268	FreeTmpSeqs( mseq2, mseq1 );
269	alloclen = strlen( result[0] ) * 2.0;
270	fprintf( stderr, "\n\ntrying to allocate TmpSeqs\n\n" );
271	AllocateTmpSeqs( &mseq2, &mseq1, alloclen );
272	}
273	for( i=0; i<locnjob; i++ ) strcpy( mseq2[i], result[i] );
274
275	if( checkC )
276	fprintf( stderr, "wm in iteration == %f\n", wm );
277
278	strcpy( mseq1, mseq2[locnjob-1] );
279	/*
280	Write( stdout, locnjob, name, nlen, mseq2 );
281	*/
282	for( i=locnjob-2; i>=s; i-- ) strcpy( mseq2[i+1], mseq2[i] );
283	strcpy( mseq2[s], mseq1 );
284	if( checkC )
285	{
286	tmpscore1= score_m_1_0( mseq2, locnjob, s, eff, effarr );
287	fprintf( stderr, "pick up %d, before ALIGNM1 score_m_1_0 = %f\n", s+1, tmpscore );
288	fprintf( stderr, "pick up %d, after ALIGNM1 score_m_1_0 = %f\n", s+1, tmpscore1 );
289	if( tmpscore1 < tmpscore )
290	{
291	fprintf( stderr, "\7" );
292	fprintf( trap_g, ">>>>>>>n\n" );
293	}
294	if( fabs( wm - tmpscore1 ) / wm > 0.001 )
295	{
296	fprintf( stderr, "\7sorry\n" );
297	exit( 1 );
298	}
299	}
300
301	identity = !strcmp( mseq2[s], aseq[s] );
302	if( s == locnjob - 1 ) ss = 0; else ss=s+1;
303
304	identity *= !strcmp( mseq2[ss], aseq[ss] );
305
306	if( !identity )
307	{
308	tmpscore = score_calc0( mseq2, locnjob, eff, s );
309	}
310	else tmpscore = tscore;
311
312	if( disp )
313	{
314	fprintf( stderr, "% 3d % 3d / the rest \n", l+1, s+1 );
315	display( mseq2, locnjob );
316	}
317	fprintf( stderr, "% 3d % 3d / the rest \n", l+1, s+1 );
318	fprintf( stderr, "score = %f mscore = %f ", tmpscore, mscore );
319
320	fprintf( trap_g, "%#4d %#4d / the rest ", l+1, s+1 );
321	fprintf( trap_g, "score = %f mscore = %f ", tmpscore, mscore );
322
323	if( identity )
324	{
325	fprintf( stderr, "( identical )\n" );
326	fprintf( trap_g, "( identical )\n" );
327	sai[s] = 2;
328	}
329
330	else if( tmpscore > mscore - cut )
331	{
332	fprintf( stderr, "accepted\n" );
333	fprintf( trap_g, "accepted\n" );
334	for( i=0; i<locnjob; i++ ) strcpy( aseq[i], mseq2[i] );
335	strcpy( sai, sai1 ); /* kokoka ? */
336	if( !tbitr && !tbweight )
337	{
338	writePre( locnjob, name, nlen, aseq, 0 );
339	}
340	strcpy( sai, sai1 );
341	tscore = tmpscore;
342	/*
343	tscore = tmpscore = score_calc0( aseq, locnjob, eff, s ); * ? *
344	*/
345	if( tmpscore > mscore )
346	{
347	for( i=0; i<locnjob; i++ ) strcpy( bseq[i], mseq2[i] );
348	treeconstruction( bseq, locnjob, topol, len, eff );
349	tscore = mscore = score_calc0( bseq, locnjob, eff, s );
350	fprintf( trap_g, " -> %f\n", mscore );
351	strcpy( sai, sai1 ); /* kokoka ? */
352	#if 0
353	if( !tbitr && !tbweight )
354	{ prep = fopen( "best", "w" );
355	Write( prep, locnjob, name, nlen, bseq );
356	fclose( prep );
357	}
358	#endif
359	}
360	}
361
362	else
363	{
364	if( tmpscore == tscore )
365	{
366	fprintf( stderr, "occational coincidence \n" );
367	fprintf( trap_g, "occational coincidence\n" );
368	}
369	else
370	{
371	fprintf( stderr, "rejected\n" );
372	fprintf( trap_g, "rejected\n" );
373	}
374	for( i=0; i<locnjob; i++ ) strcpy( aseq[i], bseq[i] );
375	tscore = mscore;
376	sai[s] = 2;
377	}
378
379	/*
380	prep = fopen( "cur", "w" );
381	Write( prep, locnjob, name, nlen, mseq2 );
382	fclose( prep );
383	*/
384
385	his[ou][s][cyc[ou]] = tmpscore;
386	if( !strcmp( sai, sai2 ) )
387	{
388	returnvalue = 0;
389	fprintf( trap_g, "converged\n" );
390	break;
391	}
392	for( i=cyc[ou]-1; i>0; i-- )
393	{
394	if( tmpscore == his[ou][s][i] )
395	{
396	shindou = 1;
397	break;
398	}
399	}
400	fprintf( stderr, "\n" );
401	if( shindou == 1 )
402	{
403	returnvalue = -1;
404	fprintf( trap_g, "oscillating\n" );
405	break;
406	}
407	}
408	if( l == MAXITERATION ) returnvalue = -2;
409	FreeDoubleCub( his );
410	return( returnvalue );
411	}
412

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

Context Navigation

source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/core/iteration.c

Download in other formats: