Context Navigation

tditeration.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: 77.7 KB

Line
1
2	/*
3	tree-dependent iteration
4	algorithm A+ when group-to-group, C when group-to-singleSeqence
5	OR
6	algorithm A+
7	*/
8
9	#include "mltaln.h"
10
11
12	#define DEBUG 0
13	#define RECORD 0
14
15	extern char **seq_g;
16	extern char **res_g;
17
18	static int nwa;
19
20	#ifdef enablemultithread
21	typedef struct _threadarg
22	{
23	int thread_no;
24	int *jobposintpt;
25	int *ndonept;
26	int *ntrypt;
27	int *collectingpt;
28	int njob;
29	int nbranch;
30	int maxiter;
31	int nkozo;
32	int *subgenerationpt;
33	float *basegainpt;
34	float *gainlist;
35	float *tscorelist;
36	int *generationofinput;
37	char *kozoarivec;
38	char **mastercopy;
39	char ***candidates;
40	int *generationofmastercopypt;
41	int *branchtable;
42	RNApair ***singlerna;
43	LocalHom **localhomtable;
44	int alloclen;
45	Node *stopol;
46	int ***topol;
47	// double **len;
48	float **tscorehistory_detail;
49	int *finishpt;
50	int **skipthisbranch;
51	pthread_mutex_t *mutex;
52	pthread_cond_t *collection_end;
53	pthread_cond_t *collection_start;
54	} threadarg_t;
55	#endif
56
57	#if 1
58	static void shuffle( int *arr, int n )
59	{
60	int i;
61	int x;
62	int b;
63
64	for( i=1; i<n; i++ )
65	{
66	x = rand() % (i+1);
67	if( x != i )
68	{
69	b = arr[i];
70	arr[i] = arr[x];
71	arr[x] = b;
72	}
73	}
74	}
75	#endif
76
77	static void Writeoption2( FILE *fp, int cycle, double cut )
78	{
79	fprintf( fp, "%dth cycle\n", cycle );
80	fprintf( fp, "marginal score to search : current score * (100-%d) / 100\n", (int)cut );
81	}
82
83	static void Writeoptions( FILE *fp )
84	{
85	fprintf( fp, "Tree-dependent-iteration\n" );
86	if( scoremtx == 1 )
87	fprintf( fp, "Blosum %d\n", nblosum );
88	else if( scoremtx == -1 )
89	fprintf( fp, "DNA\n" );
90	else if( scoremtx == 2 )
91	fprintf( fp, "Miyata-Yasunaga\n" );
92	else
93	fprintf( fp, "JTT %dPAM\n", pamN );
94
95	if( scoremtx == 0 \|\| scoremtx == 1 )
96	fprintf( fp, "Gap Penalty = %+5.3f, %5.2f, %+5.3f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
97	else
98	fprintf( fp, "Gap Penalty = %+5.3f\n", (double)penalty/1000 );
99
100
101	if( scmtd == 3 )
102	fprintf( fp, "score of rnd or sco\n" );
103	else if( scmtd == 4 )
104	fprintf( fp, "score = sigma( score for a pair of homologous amino acids ) / ( number of amino acids pairs )\n" );
105	else if( scmtd == 5 )
106	fprintf( fp, "score : SP\n" );
107	if( mix )
108	fprintf( fp, "?\n" );
109	else
110	{
111	if( weight == 2 )
112	fprintf( fp, "weighted rationale-1, geta2 = %f\n", geta2 );
113	else if( weight == 3 )
114	fprintf( fp, "weighted like ClustalW," );
115	else if( weight == 4 )
116	fprintf( fp, "weighted rationale-2, geta2 = %f\n", geta2 );
117	else
118	fprintf( fp, "unweighted\n" );
119	}
120	if( weight && utree )
121	fprintf( fp, "using tree defined by the file hat2.\n" );
122	if( weight && !utree )
123	fprintf( fp, "using temporary tree.\n" );
124
125	if( treemethod == 'n' )
126	fprintf( fp, "Tree is calculated with Neighbor-Joining Method.\n" );
127	else if( treemethod == 'q' )
128	fprintf( fp, "Tree is calculated with Minimum linkage.\n" );
129	else if( treemethod == 'X' )
130	fprintf( fp, "Tree is calculated with simplified UPG Method and UPG Method.\n" );
131	else if( treemethod == 'E' )
132	fprintf( fp, "Tree is calculated with UPG Method.\n" );
133	else
134	fprintf( fp, "Tree is calculated with unknown Method.\n" );
135
136	if( alg == 'C' )
137	fprintf( fp, "Algorithm A+ / C\n" );
138	else if( alg == 'A' )
139	fprintf( fp, "Algorithm A+ \n" );
140	else if( alg == 'a' )
141	fprintf( fp, "Algorithm A \n" );
142	else
143	fprintf( fp, "Algorithm ? \n" );
144
145	if( use_fft )
146	{
147	if( scoremtx == -1 )
148	{
149	fprintf( fp, "Basis : 4 nucleotides\n" );
150	}
151	else
152	{
153	if( fftscore )
154	fprintf( fp, "Basis : Polarity and Volume\n" );
155	else
156	fprintf( fp, "Basis : 20 amino acids\n" );
157	}
158	fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
159	fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
160	}
161	}
162
163	#ifdef enablemultithread
164
165	static void freelocalarrays(
166	float *tscorehistory,
167	RNApair *grouprna1, RNApair *grouprna2,
168	RNApair *rnapairboth,
169	char indication1, char indication2,
170	double effarr, double effarrforlocalhom, double effarr1, double effarr2,
171	char mseq1, char mseq2,
172	char **localcopy,
173	int gapmap1, int gapmap2,
174	double effarr1_kozo, double effarr2_kozo, double *effarr_kozo,
175	int **memlist,
176	char *pairbuf,
177	LocalHom *** localhomshrink
178	)
179	{
180	// fprintf( stderr, "Skipping freelocalarrays\n" );
181	// return;
182	int i;
183	if( commonIP ) FreeIntMtx( commonIP );
184	commonIP = NULL;
185	Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
186	Falign_localhom( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL,NULL, 0, NULL );
187	if( rnakozo && rnaprediction == 'm' )
188	{
189	free( grouprna1 ); // nakamimo?
190	free( grouprna2 ); // nakamimo?
191	}
192
193	free( tscorehistory );
194	free( indication1 );
195	free( indication2 );
196	free( effarr );
197	free( effarrforlocalhom );
198	free( effarr1 );
199	free( effarr2 );
200	free( mseq1 );
201	free( mseq2 );
202	FreeCharMtx( localcopy );
203	free( gapmap1 );
204	free( gapmap2 );
205
206	free( effarr1_kozo );
207	free( effarr2_kozo );
208	free( effarr_kozo );
209
210	FreeIntMtx( memlist );
211	free( pairbuf );
212
213	if( rnakozo ) free( rnapairboth );
214
215	if( constraint )
216	{
217	for( i=0; i<njob; i++)
218	{
219	free( localhomshrink[i] ); // nakamimo??
220	}
221	free( localhomshrink );
222	}
223	}
224
225
226	static void athread( void arg )
227	{
228
229	threadarg_t targ = (threadarg_t )arg;
230	int thread_no = targ->thread_no;
231	int njob = targ->njob;
232	int nbranch = targ->nbranch;
233	int maxiter = targ->maxiter;
234	int *ndonept = targ->ndonept;
235	int *ntrypt = targ->ntrypt;
236	int *collectingpt = targ->collectingpt;
237	int *jobposintpt = targ->jobposintpt;
238	int nkozo = targ->nkozo;
239	float *gainlist = targ->gainlist;
240	float *tscorelist = targ->tscorelist;
241	int *generationofinput = targ->generationofinput;
242	int *subgenerationpt = targ->subgenerationpt;
243	float *basegainpt = targ->basegainpt;
244	char *kozoarivec = targ->kozoarivec;
245	char **mastercopy = targ->mastercopy;
246	char ***candidates = targ->candidates;
247	int *generationofmastercopypt = targ->generationofmastercopypt;
248	int *branchtable = targ->branchtable;
249	RNApair ***singlerna = targ->singlerna;
250	LocalHom **localhomtable = targ->localhomtable;
251	int alloclen = targ->alloclen;
252	Node * stopol = targ->stopol;
253	int ***topol = targ->topol;
254	// double **len = targ->len;
255	float **tscorehistory_detail = targ->tscorehistory_detail;
256	int *finishpt = targ->finishpt;
257	int **skipthisbranch = targ->skipthisbranch;
258
259	int i, k, l, ii;
260	float gain;
261	int iterate;
262	int **memlist;
263	char *pairbuf;
264	int locnjob;
265	int s1, s2;
266	int clus1, clus2;
267	char **localcopy;
268	char mseq1, mseq2;
269	double effarr, effarr_kozo; // re-calc
270	double effarr1, effarr2, effarr1_kozo, effarr2_kozo;
271	char indication1, indication2;
272	int length;
273	RNApair *grouprna1, *grouprna2;
274	RNApair *rnapairboth;
275	LocalHom ***localhomshrink;
276	int gapmap1, gapmap2;
277	float tscore, mscore, oimpmatch, impmatch;
278	int identity;
279	double tmpdouble;
280	float naivescore0 = 0, naivescore1;
281	double *effarrforlocalhom;
282	float *tscorehistory;
283	int intdum;
284	#if 0
285	int oscillating;
286	int lin, ldf;
287	#endif
288	float maxgain;
289	int bestthread;
290	int branchpos;
291	int subgenerationatfirst;
292	double unweightedspscore;
293	int myjob;
294	int converged2 = 0;
295	int chudanres;
296
297
298	locnjob = njob;
299
300	if( utree == 0 )
301	{
302	fprintf( stderr, "Dynamic tree is not supported in the multithread version.\n" );
303	exit( 1 );
304	}
305	if( score_check == 2 )
306	{
307	fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
308	exit( 1 );
309	}
310
311	if( weight == 2 )
312	{
313	fprintf( stderr, "Weight 2 is not supported in the multithread version.\n" );
314	exit( 1 );
315	}
316	if( cooling && cut > 0.0 )
317	{
318	fprintf( stderr, "Cooling is not supported in the multithread version.\n" );
319	exit( 1 );
320	}
321
322	tscorehistory = calloc( maxiter, sizeof( float ) );
323
324	if( rnakozo && rnaprediction == 'm' )
325	{
326	grouprna1 = (RNApair *)calloc( njob, sizeof( RNApair ) );
327	grouprna2 = (RNApair *)calloc( njob, sizeof( RNApair ) );
328	}
329	else
330	{
331	grouprna1 = grouprna2 = NULL;
332	}
333
334	indication1 = AllocateCharVec( 150 );
335	indication2 = AllocateCharVec( 150 );
336	effarr = AllocateDoubleVec( locnjob );
337	effarrforlocalhom = AllocateDoubleVec( locnjob );
338	effarr1 = AllocateDoubleVec( locnjob );
339	effarr2 = AllocateDoubleVec( locnjob );
340	mseq1 = AllocateCharMtx( locnjob, 0 );
341	mseq2 = AllocateCharMtx( locnjob, 0 );
342	localcopy = AllocateCharMtx( locnjob, alloclen );
343	gapmap1 = AllocateIntVec( alloclen );
344	gapmap2 = AllocateIntVec( alloclen );
345
346	effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
347	effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
348	effarr_kozo = AllocateDoubleVec( locnjob );
349	for( i=0; i<locnjob; i++ )
350	effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
351
352	memlist = AllocateIntMtx( 2, locnjob );
353	pairbuf = AllocateCharVec( locnjob );
354
355
356	if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
357
358	if( constraint )
359	{
360	localhomshrink = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
361	for( i=0; i<njob; i++)
362	{
363	localhomshrink[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
364	}
365	}
366
367
368	if( thread_no == 0 )
369	{
370	*ntrypt = 0;
371	srand( randomseed );
372	*finishpt = 0;
373	for( iterate=0; iterate<maxiter; iterate++ )
374	{
375	pthread_mutex_lock( targ->mutex );
376
377	if( *collectingpt == 1 )
378	{
379	*collectingpt = 0;
380	*generationofmastercopypt = iterate;
381	*subgenerationpt = 0;
382	*basegainpt = 0.0;
383	*ndonept = 0;
384	*jobposintpt = 0;
385	for( i=0; i<nwa; i++ ) gainlist[i] = 0;
386	for( i=0; i<nwa; i++ ) tscorelist[i] = 0.0;
387	for( i=0; i<nbranch; i++ ) generationofinput[i] = -1;
388	if( parallelizationstrategy != BESTFIRST && randomseed != 0 ) shuffle( branchtable, nbranch );
389	pthread_cond_broadcast( targ->collection_end );
390	pthread_mutex_unlock( targ->mutex );
391	}
392	else
393	{
394	pthread_cond_broadcast( targ->collection_end );
395	pthread_mutex_unlock( targ->mutex );
396	freelocalarrays
397	(
398	tscorehistory,
399	grouprna1, grouprna2,
400	rnapairboth,
401	indication1, indication2,
402	effarr, effarrforlocalhom, effarr1, effarr2,
403	mseq1, mseq2,
404	localcopy,
405	gapmap1, gapmap2,
406	effarr1_kozo, effarr2_kozo, effarr_kozo,
407	memlist, pairbuf,
408	localhomshrink
409	);
410	// return( NULL );
411	pthread_exit( NULL );
412	}
413
414	pthread_mutex_lock( targ->mutex );
415	while( *ndonept < nbranch )
416	pthread_cond_wait( targ->collection_start, targ->mutex );
417	pthread_mutex_unlock( targ->mutex );
418	// fprintf( stderr, "Thread 0 got a signal, collectionpt = %d\n", collectingpt );
419
420	/*
421	Hoka no thread ga keisan
422	*/
423
424	pthread_mutex_lock( targ->mutex );
425	*collectingpt = 1; // chofuku
426
427	#if 0
428	for( i=0; i<nbranch; i++ )
429	{
430	if( generationofinput[i] != iterate )
431	{
432	fprintf( stderr, "Error! generationofinput[%d] = %d, but iterate=%d\n", i, generationofinput[i], iterate );
433	exit( 1 );
434
435	}
436	}
437	#endif
438
439	maxgain = gainlist[1];
440	bestthread = 1;
441	for( i=2; i<nwa; i++ )
442	{
443	if( gainlist[i] > maxgain )
444	{
445	maxgain = gainlist[i];
446	bestthread = i;
447	}
448	}
449
450	if( maxgain > 0.0 )
451	{
452	// fprintf( stderr, "\nGain = %f\n", maxgain );
453	// fprintf( stderr, "best gain = %f by thread %d\n", gainlist[bestthread], bestthread );
454	// fprintf( stderr, "tscorelist[best] = %f by thread %d\n", tscorelist[bestthread], bestthread );
455	if( parallelizationstrategy == BESTFIRST )
456	{
457	for( i=0; i<locnjob; i++ ) strcpy( mastercopy[i], candidates[bestthread][i] );
458	if( scoreout )
459	{
460	unweightedspscore = plainscore( locnjob, mastercopy );
461	fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * nbranch, unweightedspscore );
462	fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( locnjob * strlen( mastercopy[0] ) ) );
463	if( weight \|\| constraint ) fprintf( stderr, " (differs from the objective score)" );
464	fprintf( stderr, "\n" );
465	}
466	}
467	#if 1
468	// fprintf( stderr, "gain(%d, by %d) = %f\n", iterate, bestthread, maxgain );
469	for( i=iterate-1; i>0; i-- )
470	{
471	// if( iterate-i < 15 ) fprintf( stderr, "hist[%d] = %f\n", i, tscorehistory[i] );
472	if( tscorehistory[i] == tscorelist[bestthread] )
473	{
474	fprintf( stderr, "\nOscillating? %f == %f\n", tscorehistory[i], tscorelist[bestthread] );
475	*collectingpt = -1;
476	break;
477	}
478	}
479	tscorehistory[iterate] = tscorelist[bestthread];
480	#endif
481	}
482	else
483	{
484	fprintf( stderr, "\nConverged.\n" );
485	*collectingpt = -1;
486	// pthread_cond_broadcast( targ->collection_end );
487	// pthread_mutex_unlock( targ->mutex );
488	// freelocalarrays();
489	// return( NULL );
490	// pthread_exit( NULL );
491	}
492
493	#if 1
494	if( *finishpt )
495	{
496	fprintf( stderr, "\nConverged2.\n" );
497	*collectingpt = -1;
498	}
499	#endif
500
501	pthread_mutex_unlock( targ->mutex );
502	}
503	pthread_mutex_lock( targ->mutex );
504	fprintf( stderr, "\nReached %d\n", maxiter );
505	*collectingpt = -1;
506	pthread_cond_broadcast( targ->collection_end );
507	pthread_mutex_unlock( targ->mutex );
508	freelocalarrays
509	(
510	tscorehistory,
511	grouprna1, grouprna2,
512	rnapairboth,
513	indication1, indication2,
514	effarr, effarrforlocalhom, effarr1, effarr2,
515	mseq1, mseq2,
516	localcopy,
517	gapmap1, gapmap2,
518	effarr1_kozo, effarr2_kozo, effarr_kozo,
519	memlist, pairbuf,
520	localhomshrink
521	);
522	return( NULL );
523	pthread_exit( NULL );
524	}
525	else
526	{
527	while( 1 )
528	{
529	#if 0
530	if( iterate % 2 == 0 )
531	{
532	lin = 0; ldf = +1;
533	}
534	else
535	{
536	lin = locnjob - 2; ldf = -1;
537	}
538	for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
539	for( k=0; k<2; k++ )
540	#endif
541
542	pthread_mutex_lock( targ->mutex );
543	while( *collectingpt > 0 )
544	pthread_cond_wait( targ->collection_end, targ->mutex );
545	if( *collectingpt == -1 )
546	{
547	pthread_mutex_unlock( targ->mutex );
548	freelocalarrays
549	(
550	tscorehistory,
551	grouprna1, grouprna2,
552	rnapairboth,
553	indication1, indication2,
554	effarr, effarrforlocalhom, effarr1, effarr2,
555	mseq1, mseq2,
556	localcopy,
557	gapmap1, gapmap2,
558	effarr1_kozo, effarr2_kozo, effarr_kozo,
559	memlist, pairbuf,
560	localhomshrink
561	);
562	return( NULL );
563	pthread_exit( NULL );
564	}
565	// pthread_mutex_unlock( targ->mutex );
566
567
568	// pthread_mutex_lock( targ->mutex );
569	if( *jobposintpt == nbranch )
570	{
571	if( collectingpt != -1 ) collectingpt = 1; // chofuku
572	pthread_mutex_unlock( targ->mutex );
573	continue;
574	}
575	// fprintf( stderr, "JOB jobposintpt=%d\n", *jobposintpt );
576	myjob = branchtable[*jobposintpt];
577	l = myjob / 2;
578	if( l == locnjob-2 ) k = 1;
579	else k = myjob - l * 2;
580	// fprintf( stderr, "JOB l=%d, k=%d\n", l, k );
581
582
583	branchpos = myjob;
584	(*jobposintpt)++;
585	iterate = *generationofmastercopypt;
586	(*ntrypt)++;
587	pthread_mutex_unlock( targ->mutex );
588
589	// fprintf( stderr, "\n IRANAI IRANAI jobposintpt=%d, nbranch = %d\n", jobposintpt, nbranch );
590
591	// fprintf( stderr, "branchpos = %d (thread %d)\n", branchpos, thread_no );
592
593	// fprintf( stderr, "iterate=%d, l=%d, k=%d (thread %d)\n", iterate, l, k, thread_no );
594
595	#if 0
596	fprintf( stderr, "STEP %03d-%03d-%d (Thread %d) ", iterate+1, l+1, k, thread_no );
597	fprintf( stderr, "STEP %03d-%03d-%d (thread %d) %s ", iterate+1, l+1, k, thread_no, use_fft?"\n":"\n" );
598	#endif
599
600	// for( i=0; i<2; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
601	OneClusterAndTheOther_fast( locnjob, memlist[0], memlist[1], &s1, &s2, pairbuf, topol, l, k );
602
603
604	#if 0
605	fprintf( stderr, "STEP%d-%d\n", l, k );
606	for( i=0; i<locnjob; i++ )
607	{
608	for( j=0; j<locnjob; j++ )
609	{
610	fprintf( stderr, "%#3d", pair[i][j] );
611	}
612	fprintf( stderr, "\n" );
613	}
614	#endif
615	if( !weight )
616	{
617	for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
618	if( nkozo )
619	{
620	for( i=0; i<locnjob; i++ )
621	{
622	if( kozoarivec[i] )
623	effarr_kozo[i] = 1.0;
624	else
625	effarr_kozo[i] = 0.0;
626	}
627	}
628	}
629	else if( weight == 4 )
630	{
631	weightFromABranch( locnjob, effarr, stopol, topol, l, k );
632	if( nkozo ) // hitomadu single weight.
633	{
634	for( i=0; i<locnjob; i++ )
635	{
636	if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
637	else effarr_kozo[i] = 0.0;
638	}
639	}
640	}
641	else
642	{
643	fprintf( stderr, "weight error!\n" );
644	exit( 1 );
645	}
646
647	yarinaoshi:
648
649	pthread_mutex_lock( targ->mutex );
650	for( i=0; i<locnjob; i++ ) strcpy( localcopy[i], mastercopy[i] );
651	subgenerationatfirst = *subgenerationpt;
652	pthread_mutex_unlock( targ->mutex );
653	length = strlen( localcopy[0] );
654
655	if( nkozo )
656	{
657	// double tmptmptmp;
658	// tmptmptmp = 0.0;
659	// clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair[0], s1, localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
660	clus1 = fastconjuction_noname_kozo( memlist[0], localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
661	for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
662	// tmptmptmp = 0.0;
663	// clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair[1], s2, localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
664	clus2 = fastconjuction_noname_kozo( memlist[1], localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
665	for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
666
667	#if 0
668	fprintf( stderr, "\ngroup1 = %s\n", indication1 );
669	for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[] = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
670	fprintf( stderr, "\ngroup2 = %s\n", indication2 );
671	for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[] = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
672	#endif
673	}
674	else
675	{
676	// clus1 = conjuctionfortbfast( pair[0], s1, localcopy, mseq1, effarr1, effarr, indication1 );
677	// clus2 = conjuctionfortbfast( pair[1], s2, localcopy, mseq2, effarr2, effarr, indication2 );
678	clus1 = fastconjuction_noname( memlist[0], localcopy, mseq1, effarr1, effarr, indication1 );
679	clus2 = fastconjuction_noname( memlist[1], localcopy, mseq2, effarr2, effarr, indication2 );
680	}
681
682	if( rnakozo && rnaprediction == 'm' )
683	{
684	// makegrouprnait( grouprna1, singlerna, pair[0], s1 );
685	// makegrouprnait( grouprna2, singlerna, pair[1], s2 );
686	makegrouprna( grouprna1, singlerna, memlist[0] );
687	makegrouprna( grouprna2, singlerna, memlist[1] );
688	}
689
690	if( score_check == 2 )
691	{
692	fprintf( stderr, "Score_check 2 is not supported in the multithread version.\n" );
693	exit( 1 );
694	}
695	else if( score_check )
696	{
697	if( RNAscoremtx == 'r' )
698	intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
699	else
700	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
701
702	if( constraint )
703	{
704	// shrinklocalhom( pair, int s1, int s2, localhomtable, localhomshrink );
705	// msshrinklocalhom( pair[0], pair[1], s1, s2, localhomtable, localhomshrink );
706	msshrinklocalhom_fast( memlist[0], memlist[1], localhomtable, localhomshrink );
707	oimpmatch = 0.0;
708	if( use_fft )
709	{
710	if( alg == 'Q' )
711	{
712	part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
713	if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
714	for( i=length-1; i>=0; i-- )
715	{
716	oimpmatch += part_imp_match_out_scQ( i, i );
717	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
718	}
719	}
720	else
721	{
722	part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
723	if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
724	for( i=length-1; i>=0; i-- )
725	{
726	oimpmatch += part_imp_match_out_sc( i, i );
727	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
728	}
729	}
730	// fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
731	}
732	else
733	{
734	if( alg == 'Q' )
735	{
736	imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
737	if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
738
739	for( i=length-1; i>=0; i-- )
740	{
741	oimpmatch += imp_match_out_scQ( i, i );
742	// fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
743	}
744	}
745	else
746	{
747	imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
748
749	fprintf( stderr, "not supported\n" );
750	exit( 1 );
751
752	for( i=length-1; i>=0; i-- )
753	{
754	oimpmatch += imp_match_out_sc( i, i );
755	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
756	}
757	}
758	// fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
759	}
760	// fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
761	}
762	else
763	{
764	oimpmatch = 0.0;
765	}
766
767
768	// fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
769	mscore = (double)oimpmatch + tmpdouble;
770	}
771	else
772	{
773	fprintf( stderr, "score_check = %d\n", score_check );
774	fprintf( stderr, "Not supported\n" );
775	exit( 1 );
776	}
777
778
779	// if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
780
781	// if( !use_fft && !rnakozo )
782	if( !use_fft )
783	{
784	commongappick_record( clus1, mseq1, gapmap1 );
785	commongappick_record( clus2, mseq2, gapmap2 );
786	}
787
788	#if 0
789	fprintf( stderr, "##### mscore = %f\n", mscore );
790	#endif
791
792	#if DEBUG
793	if( !devide )
794	{
795	fprintf( trap_g, "\n%d-%d-%d\n", iterate+1, l+1, k );
796	fprintf( trap_g, "group1 = %s\n", indication1 );
797	fprintf( trap_g, "group2 = %s\n", indication2 );
798	fflush( trap_g );
799	}
800
801	#endif
802	#if 0
803	printf( "STEP %d-%d-%d\n", iterate, l, k );
804	for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
805	printf( "\n" );
806	#endif
807	if( !skipthisbranch[l][k] )
808	{
809
810	if( constraint == 2 )
811	{
812	if( use_fft )
813	{
814	// if( alg == 'Q' )
815	// part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
816	// else
817	// part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
818	chudanres = 0;
819	Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, subgenerationpt, subgenerationatfirst, &chudanres );
820	// fprintf( stderr, "##### impmatch = %f\n", impmatch );
821	if( chudanres && parallelizationstrategy == BAATARI2 )
822	{
823	// fprintf( stderr, "#### yarinaoshi!!! INS-i\n" );
824	goto yarinaoshi;
825	}
826	}
827	else
828	{
829	if( alg == 'Q' )
830	{
831	float wm;
832	// imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap wo tsukuttakara iranai.
833	// if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
834
835	wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
836	fprintf( stderr, "wm = %f\n", wm );
837	#if 0
838	fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
839	naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
840	fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
841
842	if( naivescore1 > naivescore0 )
843	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
844	else if( naivescore1 < naivescore0 )
845	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
846	else
847	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
848	#if 0 // chuui
849	if( abs( wm - naivescore1 ) > 100 )
850	{
851	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
852	// rewind( stdout );
853	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
854	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
855	// exit( 1 );
856	}
857	#endif
858	#endif
859	}
860	else if( alg == 'R' )
861	{
862	float wm;
863	imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
864	wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
865	// fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
866	naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
867	// fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
868
869	if( naivescore1 > naivescore0 )
870	fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
871	else if( naivescore1 < naivescore0 )
872	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
873	else
874	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
875	#if 0 // chuui
876	if( abs( wm - naivescore1 ) > 100 )
877	{
878	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
879	rewind( stdout );
880	for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
881	for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
882	exit( 1 );
883	}
884	#endif
885	}
886	else if( alg == 'H' )
887	{
888	float wm;
889	imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
890	wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
891	fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
892	naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
893	fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
894
895	if( naivescore1 > naivescore0 )
896	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
897	else if( naivescore1 < naivescore0 )
898	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
899	else
900	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
901	#if 0 // chuui
902	if( abs( wm - naivescore1 ) > 100 )
903	{
904	// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
905	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
906	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
907	// exit( 1 );
908	}
909	#endif
910	}
911	else
912	{
913	// imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
914	A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
915	fprintf( stderr, "A__align_gapmap\n" );
916	// fprintf( stderr, "##### impmatch = %f\n", impmatch );
917	}
918	}
919	}
920	else if( use_fft )
921	{
922	float totalwm;
923	chudanres = 0;
924	totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, subgenerationpt, subgenerationatfirst, &chudanres );
925	if( chudanres && parallelizationstrategy == BAATARI2 )
926	{
927	// fprintf( stderr, "#### yarinaoshi!!! FFT-NS-i\n" );
928	goto yarinaoshi;
929	}
930
931	// fprintf( stderr, "totalwm = %f\n", totalwm );
932	#if 0
933	if( alg == 'Q' )
934	{
935	fprintf( stderr, "totalwm = %f\n", totalwm );
936	naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
937
938	if( naivescore1 > naivescore0 )
939	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
940	else if( naivescore1 < naivescore0 )
941	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
942	else
943	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
944	#if 1 // chuui
945	if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
946	{
947	// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
948	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
949	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
950	// exit( 1 );
951	}
952	#endif
953	}
954	#endif
955	if( alg == 'R' )
956	{
957	fprintf( stderr, "totalwm = %f\n", totalwm );
958	naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
959
960	if( naivescore1 > naivescore0 )
961	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
962	else if( naivescore1 < naivescore0 )
963	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
964	else
965	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
966	#if 1 // chuui
967	if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
968	{
969	// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
970	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
971	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
972	// exit( 1 );
973	}
974	}
975	#endif
976	}
977	else
978	{
979	if( alg == 'M' )
980	{
981	chudanres = 0;
982	MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, subgenerationpt, subgenerationatfirst, &chudanres, outgap, outgap );
983	if( chudanres && parallelizationstrategy == BAATARI2 )
984	{
985	// fprintf( stderr, "#### yarinaoshi!!! NW-NS-i\n" );
986	goto yarinaoshi;
987	}
988	}
989	else if( alg == 'A' )
990	{
991	A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); //outgap==1
992	}
993	else if( alg == 'Q' )
994	{
995	float wm;
996	wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
997	fprintf( stderr, "wm = %f\n", wm );
998	fprintf( stderr, "impmatch = %f\n", impmatch );
999	naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1000
1001	if( naivescore1 > naivescore0 )
1002	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1003	else if( naivescore1 < naivescore0 )
1004	fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
1005	else
1006	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1007	#if 1 // chuui
1008	if( abs( wm - naivescore1 ) > 100 )
1009	{
1010	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
1011	// rewind( stderr );
1012	// rewind( stdout );
1013	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1014	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1015	// exit( 1 );
1016	}
1017	#endif
1018	}
1019	else if( alg == 'R' )
1020	{
1021	float wm;
1022	wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
1023	naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1024
1025	if( naivescore1 > naivescore0 )
1026	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1027	else if( naivescore1 < naivescore0 )
1028	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
1029	else
1030	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1031	#if 1 // chuui
1032	if( abs( wm - naivescore1 ) > 100 )
1033	{
1034	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
1035	// rewind( stderr );
1036	// rewind( stdout );
1037	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1038	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1039	// exit( 1 );
1040	}
1041	#endif
1042	}
1043	else if( alg == 'H' )
1044	{
1045	float wm;
1046	wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
1047	naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
1048
1049	if( naivescore1 > naivescore0 )
1050	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
1051	else if( naivescore1 < naivescore0 )
1052	{
1053	fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
1054	}
1055	else
1056	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
1057
1058	#if 0 // chuui
1059	if( abs( wm - naivescore1 ) > 100 )
1060	{
1061	rewind( stdout );
1062	for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
1063	for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
1064	exit( 1 );
1065	}
1066	#endif
1067	}
1068	else if( alg == 'a' )
1069	{
1070	Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
1071	}
1072	else ErrorExit( "Sorry!" );
1073	}
1074	// fprintf( stderr, "## impmatch = %f\n", impmatch );
1075
1076	#if 1
1077	if( parallelizationstrategy == BAATARI2 && *subgenerationpt != subgenerationatfirst )
1078	{
1079	// fprintf( stderr, "\nYarinaoshi2!! (Thread %d)\n", thread_no );
1080	goto yarinaoshi;
1081	}
1082	#endif
1083
1084	identity = !strcmp( localcopy[s1], mastercopy[s1] );
1085	identity *= !strcmp( localcopy[s2], mastercopy[s2] );
1086	fprintf( stderr, "%03d-%04d-%d (thread %4d) identical \r", iterate+1, *ndonept, k, thread_no );
1087
1088	}
1089	else
1090	{
1091	identity = 1;
1092	fprintf( stderr, "%03d-%04d-%d (thread %4d) skip \r", iterate+1, *ndonept, k, thread_no );
1093	}
1094
1095
1096	/* Bug? : idnetitcal but score change when scoreing mtx != JTT */
1097
1098	length = strlen( mseq1[0] );
1099
1100	if( identity )
1101	{
1102	tscore = mscore;
1103	}
1104	else
1105	{
1106	if( score_check )
1107	{
1108	if( constraint == 2 )
1109	{
1110	#if 1
1111	if( RNAscoremtx == 'r' )
1112	intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1113	else
1114	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1115	#else
1116	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1117	#endif
1118
1119	tscore = impmatch + tmpdouble;
1120
1121	// fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
1122	}
1123	else
1124	{
1125	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1126	tscore = tmpdouble;
1127	}
1128	// fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
1129	#if 0
1130	for( i=0; i<1; i++ )
1131	fprintf( stderr, "%s\n", mseq1[i] );
1132	fprintf( stderr, "+++++++\n" );
1133	for( i=0; i<1; i++ )
1134	fprintf( stderr, "%s\n", mseq2[i] );
1135	#endif
1136
1137	}
1138	else
1139	{
1140	tscore = mscore + 1.0;
1141	// tscore = 0.0;
1142	// fprintf( stderr, "in line 705, tscore=%f\n", tscore );
1143	// for( i=0; i<length; i++ )
1144	// tscore = tscore + (double)mseq1[0][i];
1145	// mscore = tscore - 1.0;
1146	}
1147
1148	if( isnan( mscore ) )
1149	{
1150	fprintf( stderr, "\n\nmscore became NaN\n" );
1151	exit( 1 );
1152	}
1153	if( isnan( tscore ) )
1154	{
1155	fprintf( stderr, "\n\ntscore became NaN\n" );
1156	exit( 1 );
1157	}
1158
1159
1160
1161	// fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
1162
1163	#if 1
1164	if( parallelizationstrategy == BAATARI1 && *subgenerationpt != subgenerationatfirst )
1165	{
1166	// fprintf( stderr, "\nYarinaoshi1!! (Thread %d)\n", thread_no );
1167	goto yarinaoshi;
1168	}
1169	#endif
1170	gain = tscore - ( mscore - cut/100.0*mscore );
1171	if( gain > 0 )
1172	{
1173	if( parallelizationstrategy == BESTFIRST )
1174	{
1175	if( gain > gainlist[thread_no] )
1176	{
1177	gainlist[thread_no] = gain;
1178	for( i=0; i<locnjob; i++ ) strcpy( candidates[thread_no][i], localcopy[i] );
1179	tscorelist[thread_no] = tscore;
1180	// if( iterate == 0 ) fprintf( stderr, "hist %d-%d-%d, gain=%f (Thread %d)\n", iterate, l, k, gain, thread_no );
1181	}
1182	}
1183	else
1184	{
1185	pthread_mutex_lock( targ->mutex );
1186	for( i=0; i<locnjob; i++ ) strcpy( mastercopy[i], localcopy[i] );
1187	*subgenerationpt += 1;
1188	gainlist[thread_no] = *basegainpt + gain;
1189	*basegainpt += gain;
1190
1191	if( scoreout )
1192	{
1193	unweightedspscore = plainscore( locnjob, localcopy );
1194	fprintf( stderr, "\nSCORE %d = %.0f, ", *ntrypt, unweightedspscore );
1195	fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( locnjob * strlen( mastercopy[0] ) ) );
1196	if( weight \|\| constraint ) fprintf( stderr, " (differs from the objective score)" );
1197	fprintf( stderr, "\n" );
1198	}
1199
1200	pthread_mutex_unlock( targ->mutex );
1201	tscorelist[thread_no] = tscore;
1202	}
1203	#if 0
1204	fprintf( stderr, "tscore = %f mscore = %f accepted.\n", tscore, mscore );
1205	fprintf( stderr, "\nbetter! gain = %f (thread %d)\r", gain, thread_no );
1206	#else
1207	fprintf( stderr, "%03d-%04d-%d (thread %4d) better \r", iterate+1, *ndonept, k, thread_no );
1208	#endif
1209
1210	}
1211	else
1212	{
1213	#if 0
1214	fprintf( stderr, "tscore = %f mscore = %f rejected.\r", tscore, mscore );
1215	fprintf( stderr, "worse! gain = %f", gain );
1216	#else
1217	fprintf( stderr, "%03d-%04d-%d (thread %4d) worse \r", iterate+1, *ndonept, k, thread_no );
1218	#endif
1219	tscore = mscore;
1220	}
1221	}
1222	converged2 = 0;
1223	for( ii=iterate-2; ii>=0; ii-=1 )
1224	{
1225	// fprintf( stderr, "Checking tscorehistory %f ?= %f\n", tscore, tscorehistory_detail[ii][branchpos] );
1226	if( tscore == tscorehistory_detail[ii][branchpos] )
1227	{
1228	converged2 = 1;
1229	break;
1230	}
1231	}
1232	if( parallelizationstrategy != BESTFIRST && converged2 )
1233	{
1234	// fprintf( stderr, "\nFINISH!\n" );
1235	pthread_mutex_lock( targ->mutex );
1236	*finishpt = 1;
1237	pthread_mutex_unlock( targ->mutex );
1238	}
1239
1240	tscorehistory_detail[iterate][branchpos] = tscore;
1241	fprintf( stderr, "\r" );
1242
1243	pthread_mutex_lock( targ->mutex );
1244	(*ndonept)++;
1245	// fprintf( stderr, "ndonept = %d, nbranch = %d (thread %d) iterate=%d\n", ndonept, nbranch, thread_no, iterate );
1246	generationofinput[branchpos] = iterate;
1247	if( *ndonept == nbranch )
1248	{
1249	if( collectingpt != -1 ) collectingpt = 1; // chofuku
1250	// fprintf( stderr, "Thread %d sends a signal, ndonept = %d\n", thread_no, ndonept );
1251	pthread_cond_signal( targ->collection_start );
1252	}
1253	pthread_mutex_unlock( targ->mutex );
1254	} /* while( 1 ) */
1255
1256	} /* for( iterate ) */
1257	// return( NULL );
1258	}
1259	#endif
1260
1261
1262	int TreeDependentIteration( int locnjob, char **name, int nlen[M],
1263	char aseq, char bseq, int *topol, double len,
1264	int **skipthisbranch,
1265	int alloclen, LocalHom **localhomtable,
1266	RNApair ***singlerna,
1267	int nkozo, char *kozoarivec )
1268	{
1269	int i, j, k, l, iterate, ii, iu, ju;
1270	int lin, ldf, length;
1271	int clus1, clus2;
1272	int s1, s2;
1273	static double **imanoten;
1274	static Node *stopol;
1275	static double *effarrforlocalhom = NULL;
1276	static double *effarr = NULL;
1277	static double *effarr1 = NULL;
1278	static double *effarr2 = NULL;
1279	static double *effarr_kozo = NULL;
1280	static double *effarr1_kozo = NULL;
1281	static double *effarr2_kozo = NULL;
1282	static double **mtx = NULL;
1283	static int **node = NULL;
1284	static int *branchnode = NULL;
1285	static double **branchWeight = NULL;
1286	static char mseq1, mseq2;
1287	static float ***history;
1288	FILE *trap;
1289	double tscore, mscore;
1290	int identity;
1291	int converged;
1292	int oscillating;
1293	float naivescore0 = 0.0; // by D.Mathog, a guess
1294	float naivescore1;
1295	#if 0
1296	char pair[njob][njob];
1297	#else
1298	static int **memlist;
1299	static char *pairbuf;
1300	#endif
1301	#if DEBUG + RECORD
1302	double score_for_check0, score_for_check1;
1303	static double **effmtx = NULL;
1304	extern double score_calc0();
1305	#endif
1306	static char indication1, indication2;
1307	static LocalHom ***localhomshrink = NULL;
1308	float impmatch = 0.0, oimpmatch = 0.0;
1309	static int *gapmap1;
1310	static int *gapmap2;
1311	double tmpdouble;
1312	int intdum;
1313	static RNApair *rnapairboth;
1314	RNApair *grouprna1, *grouprna2;
1315	double unweightedspscore;
1316
1317	if( rnakozo && rnaprediction == 'm' )
1318	{
1319	grouprna1 = (RNApair *)calloc( njob, sizeof( RNApair ) );
1320	grouprna2 = (RNApair *)calloc( njob, sizeof( RNApair ) );
1321	}
1322	else
1323	{
1324	grouprna1 = grouprna2 = NULL;
1325	}
1326
1327	Writeoptions( trap_g );
1328	fflush( trap_g );
1329
1330	if( effarr == NULL ) /* locnjob == njob ni kagiru */
1331	{
1332	indication1 = AllocateCharVec( 150 );
1333	indication2 = AllocateCharVec( 150 );
1334	effarr = AllocateDoubleVec( locnjob );
1335	effarrforlocalhom = AllocateDoubleVec( locnjob );
1336	effarr1 = AllocateDoubleVec( locnjob );
1337	effarr2 = AllocateDoubleVec( locnjob );
1338	mseq1 = AllocateCharMtx( locnjob, 0 );
1339	mseq2 = AllocateCharMtx( locnjob, 0 );
1340	mtx = AllocateDoubleMtx( locnjob, locnjob );
1341	node = AllocateIntMtx( locnjob, locnjob );
1342	branchnode = AllocateIntVec( locnjob );
1343	branchWeight = AllocateDoubleMtx( locnjob, 2 );
1344	history = AllocateFloatCub( niter, locnjob, 2 );
1345	stopol = (Node )calloc( locnjob 2, sizeof( Node ) );
1346	gapmap1 = AllocateIntVec( alloclen );
1347	gapmap2 = AllocateIntVec( alloclen );
1348	if( score_check == 2 ) imanoten = AllocateDoubleMtx( njob, njob );
1349
1350	effarr1_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
1351	effarr2_kozo = AllocateDoubleVec( locnjob ); // tsuneni allocate suru.
1352	effarr_kozo = AllocateDoubleVec( locnjob );
1353	for( i=0; i<locnjob; i++ )
1354	effarr_kozo[i] = effarr1_kozo[i] = effarr2_kozo[i] = 0.0;
1355
1356	#if 0
1357	#else
1358	pairbuf = AllocateCharVec( locnjob );
1359	memlist = AllocateIntMtx( 2, locnjob );
1360	if( rnakozo ) rnapairboth = (RNApair *)calloc( alloclen, sizeof( RNApair ) );
1361
1362	if( constraint )
1363	{
1364	localhomshrink = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1365	for( i=0; i<njob; i++)
1366	{
1367	localhomshrink[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1368	}
1369	}
1370	#endif
1371	}
1372	#if DEBUG + RECORD
1373	if( !effmtx ) effmtx = AllocateDoubleMtx( locnjob, locnjob );
1374	for( i=0; i<locnjob; i++ ) for( j=0; j<locnjob; j++ ) effmtx[i][j] = 1.0;
1375	#endif
1376
1377	for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
1378
1379	writePre( locnjob, name, nlen, aseq, 0 );
1380
1381	if( utree )
1382	{
1383	if( constraint )
1384	{
1385	counteff_simple( locnjob, topol, len, effarrforlocalhom );
1386	calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
1387	}
1388
1389	if( weight == 2 )
1390	{
1391	countnode_int( locnjob, topol, node );
1392	if( nkozo )
1393	{
1394	fprintf( stderr, "Not supported, weight=%d nkozo=%d.\n", weight, nkozo );
1395	}
1396	}
1397	else if( weight == 4 )
1398	{
1399	treeCnv( stopol, locnjob, topol, len, branchWeight );
1400	calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
1401	}
1402	}
1403
1404	#ifdef enablemultithread
1405	if( nthread > 0 )
1406	{
1407	threadarg_t *targ;
1408	pthread_t *handle;
1409	pthread_mutex_t mutex;
1410	pthread_cond_t collection_end;
1411	pthread_cond_t collection_start;
1412	int jobposint;
1413	int generationofmastercopy;
1414	int subgeneration;
1415	float basegain;
1416	int *generationofinput;
1417	float *gainlist;
1418	float *tscorelist;
1419	int ndone;
1420	int ntry;
1421	int collecting;
1422	int nbranch;
1423	int maxiter;
1424	char ***candidates;
1425	int *branchtable;
1426	float **tscorehistory_detail;
1427	int finish;
1428
1429	nwa = nthread + 1;
1430	nbranch = (njob-1) * 2 - 1;
1431	maxiter = niter;
1432
1433	targ = calloc( nwa, sizeof( threadarg_t ) );
1434	handle = calloc( nwa, sizeof( pthread_t ) );
1435	pthread_mutex_init( &mutex, NULL );
1436	pthread_cond_init( &collection_end, NULL );
1437	pthread_cond_init( &collection_start, NULL );
1438
1439	gainlist = calloc( nwa, sizeof( float ) );
1440	tscorelist = calloc( nwa, sizeof( float ) );
1441	branchtable = calloc( nbranch, sizeof( int ) );
1442	generationofinput = calloc( nbranch, sizeof( int ) );
1443	if( parallelizationstrategy == BESTFIRST )
1444	candidates = AllocateCharCub( nwa, locnjob, alloclen );
1445	for( i=0; i<nbranch; i++ ) branchtable[i] = i;
1446	tscorehistory_detail = AllocateFloatMtx( maxiter, nbranch );
1447
1448	collecting = 1;
1449
1450	for( i=0; i<nwa; i++ )
1451	{
1452	targ[i].thread_no = i;
1453	targ[i].njob = njob;
1454	targ[i].nbranch = nbranch;
1455	targ[i].maxiter = maxiter;
1456	targ[i].ndonept = &ndone;
1457	targ[i].ntrypt = &ntry;
1458	targ[i].collectingpt = &collecting;
1459	targ[i].jobposintpt = &jobposint;
1460	targ[i].gainlist = gainlist;
1461	targ[i].tscorelist = tscorelist;
1462	targ[i].nkozo = nkozo;
1463	targ[i].kozoarivec = kozoarivec;
1464	targ[i].mastercopy = bseq;
1465	targ[i].candidates = candidates;
1466	targ[i].subgenerationpt = &subgeneration;
1467	targ[i].basegainpt = &basegain;
1468	targ[i].generationofmastercopypt = &generationofmastercopy;
1469	targ[i].generationofinput = generationofinput;
1470	targ[i].branchtable = branchtable;
1471	targ[i].singlerna = singlerna;
1472	targ[i].localhomtable = localhomtable;
1473	targ[i].alloclen = alloclen;
1474	targ[i].stopol = stopol;
1475	targ[i].topol = topol;
1476	targ[i].skipthisbranch = skipthisbranch;
1477	// targ[i].len = len;
1478	targ[i].mutex = &mutex;
1479	targ[i].collection_end = &collection_end;
1480	targ[i].collection_start = &collection_start;
1481	targ[i].tscorehistory_detail = tscorehistory_detail;
1482	targ[i].finishpt = &finish;
1483
1484	pthread_create( handle+i, NULL, athread, (void *)(targ+i) );
1485	}
1486
1487	for( i=0; i<nwa; i++ )
1488	{
1489	pthread_join( handle[i], NULL );
1490	}
1491
1492	pthread_mutex_destroy( &mutex );
1493	pthread_cond_destroy( &collection_end );
1494	pthread_cond_destroy( &collection_start );
1495
1496	free( targ );
1497	free( handle );
1498	free( gainlist );
1499	free( tscorelist );
1500	free( branchtable );
1501	free( generationofinput );
1502	if( parallelizationstrategy == BESTFIRST )
1503	FreeCharCub( candidates );
1504	FreeFloatMtx( tscorehistory_detail );
1505	}
1506	else
1507	#endif
1508	{
1509	#if 0
1510	int *branchtable;
1511	int jobpos;
1512	int myjob;
1513
1514	int nbranch;
1515	nbranch = (njob-1) * 2 - 1;
1516
1517	branchtable = calloc( nbranch, sizeof( int ) );
1518	for( i=0; i<nbranch; i++ ) branchtable[i] = i;
1519
1520	srand( randomseed );
1521	#endif
1522
1523	if( parallelizationstrategy == BESTFIRST )
1524	{
1525	fprintf( stderr, "Not implemented. Try --thread 1 --bestfirst\n" );
1526	exit( 1 );
1527	}
1528	converged = 0;
1529	if( cooling ) cut *= 2.0;
1530	for( iterate = 0; iterate<niter; iterate++ )
1531	{
1532	if( cooling ) cut = 0.5; / ... */
1533
1534	#if 0
1535	if( randomseed != 0 ) shuffle( branchtable, nbranch );
1536	#endif
1537
1538	fprintf( trap_g, "\n" );
1539	Writeoption2( trap_g, iterate, cut );
1540	fprintf( trap_g, "\n" );
1541
1542
1543	if( utree == 0 )
1544	{
1545	if( nkozo )
1546	{
1547	fprintf( stderr, "The combination of dynamic tree and kozo is not supported.\n" );
1548	exit( 1 );
1549	}
1550	if( devide )
1551	{
1552	static char *buff1 = NULL;
1553	static char *buff2 = NULL;
1554	if( !buff1 )
1555	{
1556	buff1 = AllocateCharVec( alloclen );
1557	buff2 = AllocateCharVec( alloclen );
1558	}
1559
1560	for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
1561	{
1562	buff1[0] = buff2[0] = 0;
1563	strcat( buff1, res_g[i] ); strcat( buff2, res_g[j] );
1564	strcat( buff1, bseq[i] ); strcat( buff2, bseq[j] );
1565	strcat( buff1, seq_g[i] ); strcat( buff2, seq_g[j] );
1566
1567	mtx[i][j] = (double)substitution_hosei( buff1, buff2 );
1568	}
1569	}
1570	else
1571	{
1572	for( i=0; i<locnjob-1; i++ ) for( j=i+1; j<locnjob; j++ )
1573	mtx[i][j] = (double)substitution_hosei( bseq[i], bseq[j] );
1574	}
1575
1576	if ( treemethod == 'n' )
1577	nj( locnjob, mtx, topol, len );
1578	else if( treemethod == 's' )
1579	spg( locnjob, mtx, topol, len );
1580	else if( treemethod == 'X' )
1581	supg( locnjob, mtx, topol, len );
1582	else if( treemethod == 'p' )
1583	upg2( locnjob, mtx, topol, len );
1584	/* veryfastsupg$B$O!":#$N$H$3$m;H$($^$;$s!#(B*/
1585	/* $B=gHV$NLdBj$,$"$k$N$G(B */
1586
1587	if( weight == 2 )
1588	countnode_int( locnjob, topol, node );
1589	else if( weight == 4 )
1590	{
1591	treeCnv( stopol, locnjob, topol, len, branchWeight );
1592	calcBranchWeight( branchWeight, locnjob, stopol, topol, len );
1593	}
1594	trap = fopen( "hat2", "w" );
1595	if( !trap ) ErrorExit( "Cannot open hat2." );
1596	WriteHat2_pointer( trap, locnjob, name, mtx );
1597	fclose( trap );
1598	if( constraint )
1599	{
1600	counteff_simple( locnjob, topol, len, effarrforlocalhom );
1601	calcimportance( locnjob, effarrforlocalhom, aseq, localhomtable );
1602	}
1603	}
1604
1605	if( iterate % 2 == 0 )
1606	{
1607	lin = 0; ldf = +1;
1608	}
1609	else
1610	{
1611	lin = locnjob - 2; ldf = -1;
1612	}
1613
1614	if( score_check == 2 )
1615	{
1616	effarr1[0] = 1.0;
1617	effarr2[0] = 1.0;
1618	length = strlen( bseq[0] );
1619	for( i=0; i<locnjob-1; i++ )
1620	for( j=i+1; j<locnjob; j++ )
1621	intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
1622	}
1623
1624	#if 1
1625	for( l=lin; l < locnjob-1 && l >= 0 ; l+=ldf )
1626	{
1627
1628
1629	for( k=0; k<2; k++ )
1630	{
1631	if( l == locnjob-2 ) k = 1;
1632	#else
1633
1634	for( jobpos=0; jobpos<nbranch; jobpos++)
1635	{
1636	{
1637	myjob = branchtable[jobpos];
1638	l = myjob / 2;
1639	if( l == locnjob-2 ) k = 1;
1640	else k = myjob - l * 2;
1641	#endif
1642	#if 0 // IRANAI!!!!
1643	fprintf( stderr, "\nSTEP %d-%d\n", l, k );
1644	for( i=0; topol[l][k][i]!=-1; i++ ) fprintf( stderr, " %d ", topol[l][k][i]+1 );
1645	fprintf( stderr, "\n" );
1646	fprintf( stderr, "SKIP %d\n", skipthisbranch[l][k] );
1647	#endif
1648	#if 1
1649	fprintf( stderr, "STEP %03d-%03d-%d ", iterate+1, l+1, k );
1650	fflush( stderr );
1651	#else
1652	fprintf( stderr, "STEP %03d-%03d-%d %s", iterate+1, l+1, k, use_fft?"\n":"\n" );
1653	#endif
1654	if( skipthisbranch[l][k] )
1655	{
1656	fprintf( stderr, " skip. \r" );
1657	continue;
1658	}
1659
1660	// for( i=0; i<2; i++ ) for( j=0; j<locnjob; j++ ) pair[i][j] = 0;
1661	OneClusterAndTheOther_fast( locnjob, memlist[0], memlist[1], &s1, &s2, pairbuf, topol, l, k );
1662	#if 0 // IRANAI!!!!
1663	fprintf( stderr, "STEP%d-%d\n", l, k );
1664	for( i=0; i<2; i++ )
1665	{
1666	for( j=0; j<locnjob; j++ )
1667	{
1668	fprintf( stderr, "%#3d", pair[i][j] );
1669	}
1670	fprintf( stderr, "\n" );
1671	}
1672	#endif
1673	if( !weight )
1674	{
1675	for( i=0; i<locnjob; i++ ) effarr[i] = 1.0;
1676	if( nkozo )
1677	{
1678	for( i=0; i<locnjob; i++ )
1679	{
1680	if( kozoarivec[i] )
1681	effarr_kozo[i] = 1.0;
1682	else
1683	effarr_kozo[i] = 0.0;
1684	}
1685	}
1686	}
1687	else if( weight == 2 )
1688	{
1689	nodeFromABranch( locnjob, branchnode, node, topol, len, l, k );
1690	node_eff( locnjob, effarr, branchnode );
1691	}
1692	else if( weight == 4 )
1693	{
1694	weightFromABranch( locnjob, effarr, stopol, topol, l, k );
1695	#if 0
1696	if( nkozo )
1697	{
1698	assignstrweight( locnjob, effarr_kozo, stopol, topol, l, k, kozoarivec, effarr );
1699	}
1700
1701	#else
1702	if( nkozo ) // hitomadu single weight.
1703	for( i=0; i<locnjob; i++ )
1704	{
1705	if( kozoarivec[i] ) effarr_kozo[i] = effarr[i];
1706	else effarr_kozo[i] = 0.0;
1707	}
1708	#endif
1709	#if 0
1710	fprintf( stderr, "\n" );
1711	fprintf( stderr, "effarr_kozo = \n" );
1712	for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr_kozo[i] );
1713	fprintf( stderr, "\n" );
1714	fprintf( stderr, "effarr = \n" );
1715	for( i=0; i<locnjob; i++ ) fprintf( stderr, "%5.3f ", effarr[i] );
1716	fprintf( stderr, "\n\n" );
1717	#endif
1718	}
1719
1720	for( i=0; i<locnjob; i++ ) strcpy( aseq[i], bseq[i] );
1721	length = strlen( aseq[0] );
1722
1723	if( nkozo )
1724	{
1725	// double tmptmptmp;
1726	// tmptmptmp = 0.0;
1727	// clus1 = conjuctionfortbfast_kozo( &tmptmptmp, pair[0], s1, aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1728	clus1 = fastconjuction_noname_kozo( memlist[0], aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
1729	for( i=0; i<clus1; i++ ) effarr1_kozo[i] *= 1.0; // 0.5 ga sairyo ?
1730	// tmptmptmp = 0.0;
1731	// clus2 = conjuctionfortbfast_kozo( &tmptmptmp, pair[1], s2, aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1732	clus2 = fastconjuction_noname_kozo( memlist[1], aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
1733	for( i=0; i<clus2; i++ ) effarr2_kozo[i] *= 1.0; // 0.5 ga sairyo ?
1734
1735	#if 0
1736	fprintf( stderr, "\ngroup1 = %s\n", indication1 );
1737	for( i=0; i<clus1; i++ ) fprintf( stderr, "effarr1_kozo[%d], effarr1[] = %f, %f\n", i, effarr1_kozo[i], effarr1[i] );
1738	fprintf( stderr, "\ngroup2 = %s\n", indication2 );
1739	for( i=0; i<clus2; i++ ) fprintf( stderr, "effarr2_kozo[%d], effarr2[] = %f, %f\n", i, effarr2_kozo[i], effarr2[i] );
1740	#endif
1741
1742	}
1743	else
1744	{
1745	// clus1 = conjuctionfortbfast( pair[0], s1, aseq, mseq1, effarr1, effarr, indication1 );
1746	// clus2 = conjuctionfortbfast( pair[1], s2, aseq, mseq2, effarr2, effarr, indication2 );
1747	clus1 = fastconjuction_noname( memlist[0], aseq, mseq1, effarr1, effarr, indication1 );
1748	clus2 = fastconjuction_noname( memlist[1], aseq, mseq2, effarr2, effarr, indication2 );
1749	}
1750
1751
1752	if( rnakozo && rnaprediction == 'm' )
1753	{
1754	// makegrouprnait( grouprna1, singlerna, pair[0], s1 );
1755	// makegrouprnait( grouprna2, singlerna, pair[1], s2 );
1756	makegrouprna( grouprna1, singlerna, memlist[0] );
1757	makegrouprna( grouprna2, singlerna, memlist[1] );
1758	}
1759
1760
1761	if( score_check == 2 )
1762	{
1763	fprintf( stderr, "Not supported\n" );
1764	exit( 1 );
1765	if( constraint )
1766	{
1767	// msshrinklocalhom( pair[0], pair[1], s1, s2, localhomtable, localhomshrink );
1768	msshrinklocalhom_fast( memlist[0], memlist[1], localhomtable, localhomshrink );
1769	oimpmatch = 0.0;
1770	if( use_fft )
1771	{
1772	if( alg == 'Q' )
1773	{
1774	part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1775	if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1776	for( i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_scQ( i, i );
1777	}
1778	else
1779	{
1780	part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1781	if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1782	for( i=length-1; i>=0; i-- ) oimpmatch += part_imp_match_out_sc( i, i );
1783	}
1784	}
1785	else
1786	{
1787	if( alg == 'Q' )
1788	{
1789	imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1790	if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1791	for( i=length-1; i>=0; i-- ) oimpmatch += imp_match_out_scQ( i, i );
1792	}
1793	else
1794	{
1795	imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1796	fprintf( stderr, "not supported\n" );
1797	exit( 1 );
1798	}
1799	}
1800	// fprintf( stderr, "### oimpmatch = %f\n", oimpmatch );
1801	}
1802	else
1803	{
1804	oimpmatch = 0.0;
1805	}
1806	#if 0
1807	tmpdouble = 0.0;
1808	iu=0;
1809	for( i=s1; i<locnjob; i++ )
1810	{
1811	if( !pair[0][i] ) continue;
1812	ju=0;
1813	for( j=s2; j<locnjob; j++ )
1814	{
1815	if( !pair[1][j] ) continue;
1816	// fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
1817	tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(i,j)][MAX(i,j)];
1818	ju++;
1819	}
1820	iu++;
1821	}
1822	#else // not yet checked
1823	fprintf( stderr, "##### NOT YET CHECKED!!!!\n" );
1824	exit( 1 );
1825	tmpdouble = 0.0;
1826	iu=0;
1827	for( i=0; (s1=memlist[0][i])!=-1; i++ )
1828	{
1829	ju=0;
1830	for( j=0; (s2=memlist[1][j])!=-1; j++ )
1831	{
1832	// fprintf( stderr, "i = %d, j = %d, effarr1=%f, effarr2=%f\n", i, j, effarr1[iu], effarr2[ju] );
1833	tmpdouble += effarr1[iu] * effarr2[ju] * imanoten[MIN(s1,s2)][MAX(s1,s2)];
1834	ju++;
1835	}
1836	iu++;
1837	}
1838	#endif
1839	mscore = oimpmatch + tmpdouble;
1840	}
1841	else if( score_check )
1842	{
1843	#if 1
1844	if( RNAscoremtx == 'r' )
1845	intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1846	else
1847	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1848	#else
1849	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble ); // gappick mae denaito dame
1850	#endif
1851
1852	if( constraint )
1853	{
1854	// shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1855	msshrinklocalhom_fast( memlist[0], memlist[1], localhomtable, localhomshrink );
1856	// weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink ); // >>>
1857	oimpmatch = 0.0;
1858	if( use_fft )
1859	{
1860	if( alg == 'Q' )
1861	{
1862	part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1863	if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1864	for( i=length-1; i>=0; i-- )
1865	{
1866	oimpmatch += part_imp_match_out_scQ( i, i );
1867	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1868	}
1869	}
1870	else
1871	{
1872	part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1873	if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1874	for( i=length-1; i>=0; i-- )
1875	{
1876	oimpmatch += part_imp_match_out_sc( i, i );
1877	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1878	}
1879	}
1880	// fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
1881	}
1882	else
1883	{
1884	if( alg == 'Q' )
1885	{
1886	imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1887	if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1888
1889	for( i=length-1; i>=0; i-- )
1890	{
1891	oimpmatch += imp_match_out_scQ( i, i );
1892	// fprintf( stderr, "#### i=%d, initial impmatch = %f\n", i, oimpmatch );
1893	}
1894	}
1895	else
1896	{
1897	imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1898
1899	fprintf( stderr, "not supported\n" );
1900	exit( 1 );
1901
1902	for( i=length-1; i>=0; i-- )
1903	{
1904	oimpmatch += imp_match_out_sc( i, i );
1905	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1906	}
1907	}
1908	// fprintf( stderr, "otmpmatch = %f\n", oimpmatch );
1909	}
1910	// fprintf( stderr, "#### initial impmatch = %f\n", oimpmatch );
1911	}
1912	else
1913	{
1914	oimpmatch = 0.0;
1915	}
1916
1917
1918	// fprintf( stderr, "#### tmpdouble = %f\n", tmpdouble );
1919	mscore = (double)oimpmatch + tmpdouble;
1920	}
1921	else
1922	{
1923	// fprintf( stderr, "score_check = %d\n", score_check );
1924	/* atode kousokuka */
1925	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
1926	mscore = tmpdouble;
1927	/* atode kousokuka */
1928
1929	if( constraint )
1930	{
1931	oimpmatch = 0.0;
1932	// shrinklocalhom( pair, s1, s2, localhomtable, localhomshrink );
1933	msshrinklocalhom_fast( memlist[0], memlist[1], localhomtable, localhomshrink );
1934	if( use_fft )
1935	{
1936	if( alg == 'Q' )
1937	{
1938	part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1939	if( rnakozo ) part_imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1940	}
1941	else
1942	{
1943	part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1944	if( rnakozo ) part_imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1945	}
1946	}
1947	else
1948	{
1949	if( alg == 'Q' )
1950	{
1951	imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1952	if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, gapmap1, gapmap2, NULL );
1953	}
1954	else
1955	{
1956	imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1957	fprintf( stderr, "Not supported\n" );
1958	exit( 1 );
1959	}
1960	}
1961	}
1962	}
1963
1964	// oimpmatch = 0.0;
1965	if( constraint )
1966	{
1967	#if 0 // iranai
1968	if( alg == 'Q' )
1969	{
1970	imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1971	for( i=length-1; i>=0; i-- )
1972	{
1973	oimpmatch += imp_match_out_scQ( i, i );
1974	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1975	}
1976	}
1977	else
1978	{
1979	imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1980	for( i=length-1; i>=0; i-- )
1981	{
1982	oimpmatch += imp_match_out_sc( i, i );
1983	// fprintf( stderr, "#### i=%d, initial impmatch = %f seq1 = %c, seq2 = %c\n", i, oimpmatch, mseq1[0][i], mseq2[0][i] );
1984	}
1985	}
1986	#endif
1987	}
1988	#if 0
1989	if( alg == 'H' )
1990	naivescore0 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1991	else if( alg == 'Q' )
1992	naivescore0 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1993	else if( alg == 'R' )
1994	naivescore0 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + oimpmatch;
1995	#endif
1996
1997	// if( rnakozo ) foldalignedrna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, rnapairboth );
1998
1999	// if( !use_fft && !rnakozo )
2000	if( !use_fft )
2001	{
2002	commongappick_record( clus1, mseq1, gapmap1 );
2003	commongappick_record( clus2, mseq2, gapmap2 );
2004	}
2005
2006	#if 0
2007	fprintf( stderr, "##### mscore = %f\n", mscore );
2008	#endif
2009
2010	#if DEBUG
2011	if( !devide )
2012	{
2013	fprintf( trap_g, "\nSTEP%d-%d-%d\n", iterate+1, l+1, k );
2014	fprintf( trap_g, "group1 = %s\n", indication1 );
2015	fprintf( trap_g, "group2 = %s\n", indication2 );
2016	fflush( trap_g );
2017	}
2018
2019	#endif
2020	#if 0
2021	printf( "STEP %d-%d-%d\n", iterate, l, k );
2022	for( i=0; i<clus2; i++ ) printf( "%f ", effarr2[i] );
2023	printf( "\n" );
2024	#endif
2025	if( constraint == 2 )
2026	{
2027	if( use_fft )
2028	{
2029	// if( alg == 'Q' )
2030	// part_imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
2031	// else
2032	// part_imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
2033	Falign_localhom( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2, NULL, 0, NULL );
2034	// fprintf( stderr, "##### impmatch = %f\n", impmatch );
2035	}
2036	else
2037	{
2038	if( alg == 'Q' )
2039	{
2040	float wm;
2041	// imp_match_init_strictQ( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap wo tsukuttakara iranai.
2042	// if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, gapmap1, gapmap2, rnapairboth );
2043
2044	wm = Q__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL, gapmap1, gapmap2 );
2045	fprintf( stderr, "wm = %f\n", wm );
2046	#if 0
2047	fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2048	naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2049	fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2050
2051	if( naivescore1 > naivescore0 )
2052	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2053	else if( naivescore1 < naivescore0 )
2054	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2055	else
2056	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2057	#if 0 // chuui
2058	if( abs( wm - naivescore1 ) > 100 )
2059	{
2060	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2061	// rewind( stdout );
2062	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2063	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2064	// exit( 1 );
2065	}
2066	#endif
2067	#endif
2068	}
2069	else if( alg == 'R' )
2070	{
2071	float wm;
2072	imp_match_init_strictR( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
2073	wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
2074	// fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2075	naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2076	// fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2077
2078	if( naivescore1 > naivescore0 )
2079	fprintf( stderr, "%d-%d, ns: %f->%f UP!\n", clus1, clus2, naivescore0, naivescore1 );
2080	else if( naivescore1 < naivescore0 )
2081	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2082	else
2083	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2084	#if 0 // chuui
2085	if( abs( wm - naivescore1 ) > 100 )
2086	{
2087	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2088	rewind( stdout );
2089	for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2090	for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2091	exit( 1 );
2092	}
2093	#endif
2094	}
2095	else if( alg == 'H' )
2096	{
2097	float wm;
2098	imp_match_init_strictH( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 ); // Ichijiteki, gapmap ha mada
2099	wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, NULL, NULL, NULL, NULL );
2100	fprintf( stderr, "##### impmatch = %f->%f\n", oimpmatch, impmatch );
2101	naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty ) + impmatch;
2102	fprintf( stderr, "##### naivscore1 = %f\n", naivescore1 );
2103
2104	if( naivescore1 > naivescore0 )
2105	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2106	else if( naivescore1 < naivescore0 )
2107	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2108	else
2109	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2110	#if 0 // chuui
2111	if( abs( wm - naivescore1 ) > 100 )
2112	{
2113	// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2114	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2115	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2116	// exit( 1 );
2117	}
2118	#endif
2119	}
2120	else
2121	{
2122	// imp_match_init_strict( NULL, clus1, clus2, length, length, mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
2123	A__align_gapmap( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, localhomshrink, &impmatch, gapmap1, gapmap2 );
2124	// fprintf( stderr, "##### impmatch = %f\n", impmatch );
2125	}
2126	}
2127	}
2128	else if( use_fft )
2129	{
2130	float totalwm;
2131	totalwm = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, &intdum, NULL, 0, NULL );
2132
2133	// fprintf( stderr, "totalwm = %f\n", totalwm );
2134	#if 0
2135	if( alg == 'Q' )
2136	{
2137	fprintf( stderr, "totalwm = %f\n", totalwm );
2138	naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2139
2140	if( naivescore1 > naivescore0 )
2141	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2142	else if( naivescore1 < naivescore0 )
2143	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2144	else
2145	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2146	#if 1 // chuui
2147	if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
2148	{
2149	// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2150	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2151	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2152	// exit( 1 );
2153	}
2154	#endif
2155	}
2156	#endif
2157	if( alg == 'R' )
2158	{
2159	fprintf( stderr, "totalwm = %f\n", totalwm );
2160	naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2161
2162	if( naivescore1 > naivescore0 )
2163	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2164	else if( naivescore1 < naivescore0 )
2165	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2166	else
2167	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2168	#if 1 // chuui
2169	if( totalwm != 0.0 && abs( totalwm - naivescore1 ) > 100 )
2170	{
2171	// fprintf( stderr, "WARNING, totalwm=%f but naivescore=%f\n", totalwm, naivescore1 );
2172	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2173	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2174	// exit( 1 );
2175	}
2176	}
2177	#endif
2178	}
2179	else
2180	{
2181	if( alg == 'M' )
2182	{
2183	MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
2184	}
2185	else if( alg == 'A' )
2186	{
2187	A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL, NULL, 0, NULL, 1, 1 ); // outgap==1
2188	}
2189	else if( alg == 'Q' )
2190	{
2191	float wm;
2192	wm = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2193	fprintf( stderr, "wm = %f\n", wm );
2194	fprintf( stderr, "impmatch = %f\n", impmatch );
2195	naivescore1 = naiveQpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2196
2197	if( naivescore1 > naivescore0 )
2198	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2199	else if( naivescore1 < naivescore0 )
2200	fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
2201	else
2202	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2203	#if 1 // chuui
2204	if( abs( wm - naivescore1 ) > 100 )
2205	{
2206	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2207	// rewind( stderr );
2208	// rewind( stdout );
2209	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2210	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2211	// exit( 1 );
2212	}
2213	#endif
2214	}
2215	else if( alg == 'R' )
2216	{
2217	float wm;
2218	wm = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2219	naivescore1 = naiveRpairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2220
2221	if( naivescore1 > naivescore0 )
2222	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2223	else if( naivescore1 < naivescore0 )
2224	fprintf( stderr, "%d-%d, ns: DOWN! %f->%f\n", clus1, clus2, naivescore0, naivescore1 );
2225	else
2226	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2227	#if 1 // chuui
2228	if( abs( wm - naivescore1 ) > 100 )
2229	{
2230	// fprintf( stderr, "WARNING, wm=%f but naivescore=%f\n", wm, naivescore1 );
2231	// rewind( stderr );
2232	// rewind( stdout );
2233	// for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2234	// for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2235	// exit( 1 );
2236	}
2237	#endif
2238	}
2239	else if( alg == 'H' )
2240	{
2241	float wm;
2242	wm = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen, NULL, &impmatch, NULL, NULL, NULL, NULL );
2243	naivescore1 = naivepairscore( clus1, clus2, mseq1, mseq2, effarr1, effarr2, penalty );
2244
2245	if( naivescore1 > naivescore0 )
2246	fprintf( stderr, "%d-%d, ns: UP!\n", clus1, clus2 );
2247	else if( naivescore1 < naivescore0 )
2248	{
2249	fprintf( stderr, "%d-%d, ns: DOWN!\n", clus1, clus2 );
2250	}
2251	else
2252	fprintf( stderr, "%d-%d, ns: IDENTICAL\n", clus1, clus2 );
2253
2254	#if 0 // chuui
2255	if( abs( wm - naivescore1 ) > 100 )
2256	{
2257	rewind( stdout );
2258	for( i=0; i<clus1; i++ ) printf( ">\n%s\n", mseq1[i] );
2259	for( i=0; i<clus2; i++ ) printf( ">\n%s\n", mseq2[i] );
2260	exit( 1 );
2261	}
2262	#endif
2263	}
2264	else if( alg == 'a' )
2265	{
2266	Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, alloclen );
2267	}
2268	else ErrorExit( "Sorry!" );
2269	}
2270	// fprintf( stderr, "## impmatch = %f\n", impmatch );
2271
2272	if( checkC )
2273	{
2274	extern double DSPscore();
2275	extern double SSPscore();
2276	static double cur;
2277	static double pre;
2278
2279	/*
2280	pre = SSPscore( locnjob, bseq );
2281	cur = SSPscore( locnjob, aseq );
2282	*/
2283	pre = DSPscore( locnjob, bseq );
2284	cur = DSPscore( locnjob, aseq );
2285
2286	fprintf( stderr, "Previous Sscore = %f\n", pre );
2287	fprintf( stderr, "Currnet Sscore = %f\n\n", cur );
2288	}
2289
2290	// fprintf( stderr, "## impmatch = %f\n", impmatch );
2291	identity = !strcmp( aseq[s1], bseq[s1] );
2292	identity *= !strcmp( aseq[s2], bseq[s2] );
2293
2294
2295	/* Bug? : idnetitcal but score change when scoreing mtx != JTT */
2296
2297	length = strlen( mseq1[0] );
2298
2299	if( identity )
2300	{
2301	tscore = mscore;
2302	if( !devide ) fprintf( trap_g, "tscore = %f identical.\n", tscore );
2303	fprintf( stderr, " identical. " );
2304	converged++;
2305	}
2306	else
2307	{
2308	if( score_check )
2309	{
2310	if( constraint == 2 )
2311	{
2312	#if 1
2313	if( RNAscoremtx == 'r' )
2314	intergroup_score_gapnomi( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2315	else
2316	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2317	#else
2318	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2319	#endif
2320
2321	tscore = impmatch + tmpdouble;
2322
2323	// fprintf( stderr, "tmpdouble=%f, impmatch = %f -> %f, tscore = %f\n", tmpdouble, oimpmatch, impmatch, tscore );
2324	}
2325	else
2326	{
2327	intergroup_score( mseq1, mseq2, effarr1, effarr2, clus1, clus2, length, &tmpdouble );
2328	tscore = tmpdouble;
2329	}
2330	// fprintf( stderr, "#######ii=%d, iterate=%d score = %f -> %f \n", ii, iterate , mscore, tscore );
2331	#if 0
2332	for( i=0; i<1; i++ )
2333	fprintf( stderr, "%s\n", mseq1[i] );
2334	fprintf( stderr, "+++++++\n" );
2335	for( i=0; i<1; i++ )
2336	fprintf( stderr, "%s\n", mseq2[i] );
2337	#endif
2338
2339	}
2340	else
2341	{
2342	tscore = mscore + 1.0;
2343	// tscore = 0.0;
2344	// fprintf( stderr, "in line 705, tscore=%f\n", tscore );
2345	// for( i=0; i<length; i++ )
2346	// tscore = tscore + (double)mseq1[0][i];
2347	// mscore = tscore - 1.0;
2348	}
2349
2350	if( isnan( mscore ) )
2351	{
2352	fprintf( stderr, "\n\nmscore became NaN\n" );
2353	exit( 1 );
2354	}
2355	if( isnan( tscore ) )
2356	{
2357	fprintf( stderr, "\n\ntscore became NaN\n" );
2358	exit( 1 );
2359	}
2360
2361
2362
2363	// fprintf( stderr, "@@@@@ mscore,tscore = %f,%f\n", mscore, tscore );
2364
2365	if( tscore > mscore - cut/100.0*mscore )
2366	{
2367	writePre( locnjob, name, nlen, aseq, 0 );
2368	for( i=0; i<locnjob; i++ ) strcpy( bseq[i], aseq[i] );
2369	if( score_check == 2 )
2370	{
2371	effarr1[0] = 1.0;
2372	effarr2[0] = 1.0;
2373	for( i=0; i<locnjob-1; i++ )
2374	for( j=i+1; j<locnjob; j++ )
2375	intergroup_score( bseq+i, bseq+j, effarr1, effarr2, 1, 1, length, imanoten[i]+j );
2376	}
2377
2378	#if 0
2379	fprintf( stderr, "tscore = %f mscore = %f accepted.\n", tscore, mscore );
2380	#endif
2381	fprintf( stderr, " accepted." );
2382	converged = 0;
2383
2384	}
2385	else
2386	{
2387	#if 0
2388	fprintf( stderr, "tscore = %f mscore = %f rejected.\n", tscore, mscore );
2389	#endif
2390	fprintf( stderr, " rejected." );
2391	tscore = mscore;
2392	converged++;
2393	}
2394	}
2395	fprintf( stderr, "\r" );
2396
2397
2398	history[iterate][l][k] = (float)tscore;
2399
2400	// fprintf( stderr, "tscore = %f\n", tscore );
2401
2402	if( converged >= locnjob * 2 )
2403	{
2404	fprintf( trap_g, "Converged.\n\n" );
2405	fprintf( stderr, "\nConverged.\n\n" );
2406	if( scoreout )
2407	{
2408	unweightedspscore = plainscore( njob, bseq );
2409	fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2410	fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2411	if( weight \|\| constraint ) fprintf( stderr, " (differs from the objective score)" );
2412	fprintf( stderr, "\n\n" );
2413	}
2414	if( grouprna1 ) free( grouprna1 );
2415	if( grouprna2 ) free( grouprna2 );
2416	return( 0 );
2417	}
2418	if( iterate >= 1 )
2419	{
2420	/* oscillation check */
2421	oscillating = 0;
2422	for( ii=iterate-2; ii>=0; ii-=2 )
2423	{
2424	if( (float)tscore == history[ii][l][k] )
2425	{
2426	oscillating = 1;
2427	break;
2428	}
2429	}
2430	if( ( oscillating && !cooling ) \|\| ( oscillating && cut < 0.001 && cooling ) )
2431	{
2432	fprintf( trap_g, "Oscillating.\n" );
2433	fprintf( stderr, "\nOscillating.\n\n" );
2434	if( scoreout )
2435	{
2436	unweightedspscore = plainscore( njob, bseq );
2437	fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2438	fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2439	if( weight \|\| constraint ) fprintf( stderr, " (differs from the objective score)" );
2440	fprintf( stderr, "\n\n" );
2441	}
2442	#if 1 /* hujuubun */
2443	if( grouprna1 ) free( grouprna1 );
2444	if( grouprna2 ) free( grouprna2 );
2445	return( -1 );
2446	#endif
2447	}
2448	} /* if( iterate ) */
2449	} /* for( k ) */
2450	} /* for( l ) */
2451	if( scoreout )
2452	{
2453	unweightedspscore = plainscore( njob, bseq );
2454	fprintf( stderr, "\nSCORE %d = %.0f, ", iterate * ( (njob-1)*2-1 ), unweightedspscore );
2455	fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2456	if( weight \|\| constraint ) fprintf( stderr, " (differs from the objective score)" );
2457	fprintf( stderr, "\n\n" );
2458	}
2459	} /* for( iterate ) */
2460	}
2461	if( grouprna1 ) free( grouprna1 );
2462	if( grouprna2 ) free( grouprna2 );
2463	return( 2 );
2464	} /* int Tree... */

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

Context Navigation

source: branches/stable/GDE/MAFFT/mafft-7.055-with-extensions/core/tditeration.c

Download in other formats: