Context Navigation

tbfast.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: 57.9 KB

Line
1	#include "mltaln.h"
2
3	#define DEBUG 0
4	#define IODEBUG 0
5	#define SCOREOUT 0
6
7	static int nadd;
8	static int treein;
9	static int topin;
10	static int treeout;
11	static int distout;
12	static int noalign;
13	static int multidist;
14	static int subalignment;
15	static int subalignmentoffset;
16
17	#ifdef enablemultithread
18	typedef struct _jobtable
19	{
20	int i;
21	int j;
22	} Jobtable;
23
24	typedef struct _distancematrixthread_arg
25	{
26	int njob;
27	int thread_no;
28	float *selfscore;
29	float **iscore;
30	char **seq;
31	Jobtable *jobpospt;
32	pthread_mutex_t *mutex;
33	} distancematrixthread_arg_t;
34
35	typedef struct _treebasethread_arg
36	{
37	int thread_no;
38	int *nrunpt;
39	int njob;
40	int *nlen;
41	int *jobpospt;
42	int ***topol;
43	Treedep *dep;
44	char **aseq;
45	double *effarr;
46	int *alloclenpt;
47	LocalHom **localhomtable;
48	RNApair ***singlerna;
49	double *effarr_kozo;
50	int *fftlog;
51	char *mergeoralign;
52	pthread_mutex_t *mutex;
53	pthread_cond_t *treecond;
54	} treebasethread_arg_t;
55	#endif
56
57	void arguments( int argc, char *argv[] )
58	{
59	int c;
60
61	nthread = 1;
62	outnumber = 0;
63	scoreout = 0;
64	treein = 0;
65	topin = 0;
66	rnaprediction = 'm';
67	rnakozo = 0;
68	nevermemsave = 0;
69	inputfile = NULL;
70	addfile = NULL;
71	addprofile = 1;
72	fftkeika = 0;
73	constraint = 0;
74	nblosum = 62;
75	fmodel = 0;
76	calledByXced = 0;
77	devide = 0;
78	use_fft = 0; // chuui
79	force_fft = 0;
80	fftscore = 1;
81	fftRepeatStop = 0;
82	fftNoAnchStop = 0;
83	weight = 3;
84	utree = 1;
85	tbutree = 1;
86	refine = 0;
87	check = 1;
88	cut = 0.0;
89	disp = 0;
90	outgap = 1;
91	alg = 'A';
92	mix = 0;
93	tbitr = 0;
94	scmtd = 5;
95	tbweight = 0;
96	tbrweight = 3;
97	checkC = 0;
98	treemethod = 'X';
99	contin = 0;
100	scoremtx = 1;
101	kobetsubunkatsu = 0;
102	dorp = NOTSPECIFIED;
103	ppenalty = NOTSPECIFIED;
104	ppenalty_ex = NOTSPECIFIED;
105	poffset = NOTSPECIFIED;
106	kimuraR = NOTSPECIFIED;
107	pamN = NOTSPECIFIED;
108	geta2 = GETA2;
109	fftWinSize = NOTSPECIFIED;
110	fftThreshold = NOTSPECIFIED;
111	RNAppenalty = NOTSPECIFIED;
112	RNAppenalty_ex = NOTSPECIFIED;
113	RNApthr = NOTSPECIFIED;
114	TMorJTT = JTT;
115	consweight_multi = 1.0;
116	consweight_rna = 0.0;
117	multidist = 0;
118	subalignment = 0;
119	subalignmentoffset = 0;
120
121	while( --argc > 0 && (*++argv)[0] == '-' )
122	{
123	while ( ( c = *++argv[0] ) )
124	{
125	switch( c )
126	{
127	case 'i':
128	inputfile = *++argv;
129	fprintf( stderr, "inputfile = %s\n", inputfile );
130	--argc;
131	goto nextoption;
132	case 'I':
133	nadd = atoi( *++argv );
134	fprintf( stderr, "nadd = %d\n", nadd );
135	--argc;
136	goto nextoption;
137	case 'e':
138	RNApthr = (int)( atof( ++argv ) 1000 - 0.5 );
139	--argc;
140	goto nextoption;
141	case 'o':
142	RNAppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
143	--argc;
144	goto nextoption;
145	case 'f':
146	ppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
147	// fprintf( stderr, "ppenalty = %d\n", ppenalty );
148	--argc;
149	goto nextoption;
150	case 'g':
151	ppenalty_ex = (int)( atof( ++argv ) 1000 - 0.5 );
152	fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
153	--argc;
154	goto nextoption;
155	case 'h':
156	poffset = (int)( atof( ++argv ) 1000 - 0.5 );
157	// fprintf( stderr, "poffset = %d\n", poffset );
158	--argc;
159	goto nextoption;
160	case 'k':
161	kimuraR = atoi( *++argv );
162	fprintf( stderr, "kappa = %d\n", kimuraR );
163	--argc;
164	goto nextoption;
165	case 'b':
166	nblosum = atoi( *++argv );
167	scoremtx = 1;
168	fprintf( stderr, "blosum %d / kimura 200\n", nblosum );
169	--argc;
170	goto nextoption;
171	case 'j':
172	pamN = atoi( *++argv );
173	scoremtx = 0;
174	TMorJTT = JTT;
175	fprintf( stderr, "jtt/kimura %d\n", pamN );
176	--argc;
177	goto nextoption;
178	case 'm':
179	pamN = atoi( *++argv );
180	scoremtx = 0;
181	TMorJTT = TM;
182	fprintf( stderr, "tm %d\n", pamN );
183	--argc;
184	goto nextoption;
185	case 'l':
186	fastathreshold = atof( *++argv );
187	constraint = 2;
188	--argc;
189	goto nextoption;
190	case 'r':
191	consweight_rna = atof( *++argv );
192	rnakozo = 1;
193	--argc;
194	goto nextoption;
195	case 'c':
196	consweight_multi = atof( *++argv );
197	--argc;
198	goto nextoption;
199	case 'C':
200	nthread = atoi( *++argv );
201	fprintf( stderr, "nthread = %d\n", nthread );
202	--argc;
203	goto nextoption;
204	case 'H':
205	subalignment = 1;
206	subalignmentoffset = atoi( *++argv );
207	--argc;
208	goto nextoption;
209	case 'R':
210	rnaprediction = 'r';
211	break;
212	case 's':
213	RNAscoremtx = 'r';
214	break;
215	#if 1
216	case 'a':
217	fmodel = 1;
218	break;
219	#endif
220	case 'K':
221	addprofile = 0;
222	break;
223	case 'y':
224	distout = 1;
225	break;
226	case 't':
227	treeout = 1;
228	break;
229	case 'T':
230	noalign = 1;
231	break;
232	case 'D':
233	dorp = 'd';
234	break;
235	case 'P':
236	dorp = 'p';
237	break;
238	#if 1
239	case 'O':
240	outgap = 0;
241	break;
242	#else
243	case 'O':
244	fftNoAnchStop = 1;
245	break;
246	#endif
247	case 'S':
248	scoreout = 1;
249	break;
250	#if 0
251	case 'e':
252	fftscore = 0;
253	break;
254	case 'r':
255	fmodel = -1;
256	break;
257	case 'R':
258	fftRepeatStop = 1;
259	break;
260	case 's':
261	treemethod = 's';
262	break;
263	#endif
264	case 'X':
265	treemethod = 'X';
266	break;
267	case 'E':
268	treemethod = 'E';
269	break;
270	case 'q':
271	treemethod = 'q';
272	break;
273	case 'n' :
274	outnumber = 1;
275	break;
276	#if 0
277	case 'a':
278	alg = 'a';
279	break;
280	case 'H':
281	alg = 'H';
282	break;
283	#endif
284	case 'Q':
285	alg = 'Q';
286	break;
287	case 'A':
288	alg = 'A';
289	break;
290	case 'M':
291	alg = 'M';
292	break;
293	case 'N':
294	nevermemsave = 1;
295	break;
296	case 'B':
297	break;
298	case 'F':
299	use_fft = 1;
300	break;
301	case 'G':
302	force_fft = 1;
303	use_fft = 1;
304	break;
305	case 'U':
306	treein = 1;
307	break;
308	case 'V':
309	topin = 1;
310	break;
311	case 'u':
312	tbrweight = 0;
313	weight = 0;
314	break;
315	case 'v':
316	tbrweight = 3;
317	break;
318	case 'd':
319	multidist = 1;
320	break;
321	#if 0
322	case 'd':
323	disp = 1;
324	break;
325	#endif
326	/* Modified 01/08/27, default: user tree */
327	case 'J':
328	tbutree = 0;
329	break;
330	/* modification end. */
331	case 'z':
332	fftThreshold = atoi( *++argv );
333	--argc;
334	goto nextoption;
335	case 'w':
336	fftWinSize = atoi( *++argv );
337	--argc;
338	goto nextoption;
339	case 'Z':
340	checkC = 1;
341	break;
342	default:
343	fprintf( stderr, "illegal option %c\n", c );
344	argc = 0;
345	break;
346	}
347	}
348	nextoption:
349	;
350	}
351	if( argc == 1 )
352	{
353	cut = atof( (*argv) );
354	argc--;
355	}
356	if( argc != 0 )
357	{
358	fprintf( stderr, "options: Check source file !\n" );
359	exit( 1 );
360	}
361	if( tbitr == 1 && outgap == 0 )
362	{
363	fprintf( stderr, "conflicting options : o, m or u\n" );
364	exit( 1 );
365	}
366	if( alg == 'C' && outgap == 0 )
367	{
368	fprintf( stderr, "conflicting options : C, o\n" );
369	exit( 1 );
370	}
371	}
372
373	#if 0
374	static void distancematrixthread2( void arg )
375	{
376	distancematrixthread_arg_t targ = (distancematrixthread_arg_t )arg;
377	int njob = targ->njob;
378	int thread_no = targ->thread_no;
379	float *selfscore = targ->selfscore;
380	float **iscore = targ->iscore;
381	char **seq = targ->seq;
382	Jobtable *jobpospt = targ->jobpospt;
383
384	float ssi, ssj, bunbo;
385	int i, j;
386
387	while( 1 )
388	{
389	pthread_mutex_lock( targ->mutex );
390	i = jobpospt->i;
391	i++;
392	if( i == njob-1 )
393	{
394	pthread_mutex_unlock( targ->mutex );
395	return( NULL );
396	}
397	jobpospt->i = i;
398	pthread_mutex_unlock( targ->mutex );
399
400	ssi = selfscore[i];
401	if( i % 10 == 0 ) fprintf( stderr, "\r% 5d / %d (thread %4d)", i, njob, thread_no );
402	for( j=i+1; j<njob; j++ )
403	{
404	ssj = selfscore[j];
405	bunbo = MIN( ssi, ssj );
406	if( bunbo == 0.0 )
407	iscore[i][j-i] = 1.0;
408	else
409	iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
410	}
411	}
412	}
413	#endif
414
415	#ifdef enablemultithread
416	static void distancematrixthread( void arg )
417	{
418	distancematrixthread_arg_t targ = (distancematrixthread_arg_t )arg;
419	int njob = targ->njob;
420	int thread_no = targ->thread_no;
421	float *selfscore = targ->selfscore;
422	float **iscore = targ->iscore;
423	char **seq = targ->seq;
424	Jobtable *jobpospt = targ->jobpospt;
425
426	float ssi, ssj, bunbo;
427	int i, j;
428
429	while( 1 )
430	{
431	pthread_mutex_lock( targ->mutex );
432	j = jobpospt->j;
433	i = jobpospt->i;
434	j++;
435	if( j == njob )
436	{
437	i++;
438	j = i + 1;
439	if( i == njob-1 )
440	{
441	pthread_mutex_unlock( targ->mutex );
442	return( NULL );
443	}
444	}
445	jobpospt->j = j;
446	jobpospt->i = i;
447	pthread_mutex_unlock( targ->mutex );
448
449
450	if( j==i+1 && i % 10 == 0 ) fprintf( stderr, "\r% 5d / %d (thread %4d)", i, njob, thread_no );
451	ssi = selfscore[i];
452	ssj = selfscore[j];
453	bunbo = MIN( ssi, ssj );
454	if( bunbo == 0.0 )
455	iscore[i][j-i] = 1.0;
456	else
457	iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
458	}
459	}
460
461	static void treebasethread( void arg )
462	{
463	treebasethread_arg_t targ = (treebasethread_arg_t )arg;
464	int *nrunpt = targ->nrunpt;
465	int thread_no = targ->thread_no;
466	int njob = targ->njob;
467	int *nlen = targ->nlen;
468	int *jobpospt = targ->jobpospt;
469	int ***topol = targ->topol;
470	Treedep *dep = targ->dep;
471	char **aseq = targ->aseq;
472	double *effarr = targ->effarr;
473	int *alloclen = targ->alloclenpt;
474	LocalHom **localhomtable = targ->localhomtable;
475	RNApair ***singlerna = targ->singlerna;
476	double *effarr_kozo = targ->effarr_kozo;
477	int *fftlog = targ->fftlog;
478	char *mergeoralign = targ->mergeoralign;
479
480	char mseq1, mseq2;
481	char **localcopy;
482	int i, j, l;
483	int len1, len2;
484	int clus1, clus2;
485	float pscore;
486	char indication1, indication2;
487	double *effarr1 = NULL;
488	double *effarr2 = NULL;
489	double *effarr1_kozo = NULL;
490	double *effarr2_kozo = NULL;
491	LocalHom ***localhomshrink = NULL;
492	int m1, m2;
493	float dumfl = 0.0;
494	int ffttry;
495	RNApair *grouprna1 = NULL, *grouprna2 = NULL;
496
497
498	mseq1 = AllocateCharMtx( njob, 0 );
499	mseq2 = AllocateCharMtx( njob, 0 );
500	localcopy = calloc( njob, sizeof( char * ) );
501
502	if( rnakozo && rnaprediction == 'm' )
503	{
504	grouprna1 = (RNApair *)calloc( njob, sizeof( RNApair ) );
505	grouprna2 = (RNApair *)calloc( njob, sizeof( RNApair ) );
506	}
507	else
508	{
509	grouprna1 = grouprna2 = NULL;
510	}
511
512	effarr1 = AllocateDoubleVec( njob );
513	effarr2 = AllocateDoubleVec( njob );
514	indication1 = AllocateCharVec( 150 );
515	indication2 = AllocateCharVec( 150 );
516	#if 0
517	#else
518	if( constraint )
519	{
520	localhomshrink = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
521	for( i=0; i<njob; i++ )
522	localhomshrink[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
523	}
524	#endif
525	effarr1_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
526	effarr2_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
527	for( i=0; i<njob; i++ ) effarr1_kozo[i] = 0.0;
528	for( i=0; i<njob; i++ ) effarr2_kozo[i] = 0.0;
529
530
531	#if 0
532	#endif
533
534	#if 0
535	for( i=0; i<njob; i++ )
536	fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
537	#endif
538
539	#if 0 //-> main thread
540	if( constraint )
541	calcimportance( njob, effarr, aseq, localhomtable );
542	#endif
543
544
545	// writePre( njob, name, nlen, aseq, 0 );
546
547
548	// for( l=0; l<njob-1; l++ )
549	while( 1 )
550	{
551
552	pthread_mutex_lock( targ->mutex );
553	l = *jobpospt;
554	if( l == njob-1 )
555	{
556	pthread_mutex_unlock( targ->mutex );
557	if( commonIP ) FreeIntMtx( commonIP );
558	commonIP = NULL;
559	Falign( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, 0, NULL );
560	Falign_udpari_long( NULL, NULL, NULL, NULL, 0, 0, 0, NULL );
561	A__align( NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0, NULL, 0, 0 );
562	free( mseq1 );
563	free( mseq2 );
564	free( localcopy );
565	free( effarr1 );
566	free( effarr2 );
567	free( effarr1_kozo );
568	free( effarr2_kozo );
569	free( indication1 );
570	free( indication2 );
571	if( rnakozo && rnaprediction == 'm' )
572	{
573	if( grouprna1 ) free( grouprna1 ); // nakami ha?
574	if( grouprna2 ) free( grouprna2 ); // nakami ha?
575	grouprna1 = grouprna2 = NULL;
576	}
577	if( constraint )
578	{
579	if( localhomshrink ) // nen no tame
580	{
581	for( i=0; i<njob; i++ )
582	{
583	free( localhomshrink[i] );
584	localhomshrink[i] = NULL;
585	}
586	free( localhomshrink );
587	localhomshrink = NULL;
588	}
589	}
590	return( NULL );
591	}
592	*jobpospt = l+1;
593
594	if( dep[l].child0 != -1 )
595	{
596	while( dep[dep[l].child0].done == 0 )
597	pthread_cond_wait( targ->treecond, targ->mutex );
598	}
599	if( dep[l].child1 != -1 )
600	{
601	while( dep[dep[l].child1].done == 0 )
602	pthread_cond_wait( targ->treecond, targ->mutex );
603	}
604	// while( *nrunpt >= nthread )
605	// pthread_cond_wait( targ->treecond, targ->mutex );
606	(*nrunpt)++;
607
608	// pthread_mutex_unlock( targ->mutex );
609
610	if( mergeoralign[l] == 'n' )
611	{
612	// fprintf( stderr, "SKIP!\n" );
613	dep[l].done = 1;
614	(*nrunpt)--;
615	pthread_cond_broadcast( targ->treecond );
616	free( topol[l][0] );
617	free( topol[l][1] );
618	free( topol[l] );
619	pthread_mutex_unlock( targ->mutex );
620	continue;
621	}
622
623
624
625	m1 = topol[l][0][0];
626	m2 = topol[l][1][0];
627
628	// pthread_mutex_lock( targ->mutex );
629
630
631
632	len1 = strlen( aseq[m1] );
633	len2 = strlen( aseq[m2] );
634	if( *alloclen <= len1 + len2 )
635	{
636	fprintf( stderr, "\nReallocating (by thread %d) ..", thread_no );
637	*alloclen = ( len1 + len2 ) + 1000;
638	ReallocateCharMtx( aseq, njob, *alloclen + 10 );
639	fprintf( stderr, "done. alloclen = %d\n", alloclen );
640	}
641
642	for( i=0; (j=topol[l][0][i])!=-1; i++ )
643	{
644	localcopy[j] = calloc( *alloclen, sizeof( char ) );
645	strcpy( localcopy[j], aseq[j] );
646	}
647	for( i=0; (j=topol[l][1][i])!=-1; i++ )
648	{
649	localcopy[j] = calloc( *alloclen, sizeof( char ) );
650	strcpy( localcopy[j], aseq[j] );
651	}
652
653	pthread_mutex_unlock( targ->mutex );
654
655
656
657	if( effarr_kozo )
658	{
659	clus1 = fastconjuction_noname_kozo( topol[l][0], localcopy, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
660	clus2 = fastconjuction_noname_kozo( topol[l][1], localcopy, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
661	}
662	else
663	{
664	clus1 = fastconjuction_noname( topol[l][0], localcopy, mseq1, effarr1, effarr, indication1 );
665	clus2 = fastconjuction_noname( topol[l][1], localcopy, mseq2, effarr2, effarr, indication2 );
666	}
667
668
669
670	#if 1
671	fprintf( stderr, "\rSTEP % 5d /%d (thread %4d) ", l+1, njob-1, thread_no );
672	#else
673	fprintf( stderr, "STEP %d /%d (thread %d) \n", l+1, njob-1, thread_no );
674	fprintf( stderr, "group1 = %.66s", indication1 );
675	if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
676	fprintf( stderr, ", child1 = %d\n", dep[l].child0 );
677	fprintf( stderr, "group2 = %.66s", indication2 );
678	if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
679	fprintf( stderr, ", child2 = %d\n", dep[l].child1 );
680
681	fprintf( stderr, "Group1's lengths = " );
682	for( i=0; i<clus1; i++ ) fprintf( stderr, "%d ", strlen( mseq1[i] ) );
683	fprintf( stderr, "\n" );
684	fprintf( stderr, "Group2's lengths = " );
685	for( i=0; i<clus2; i++ ) fprintf( stderr, "%d ", strlen( mseq2[i] ) );
686	fprintf( stderr, "\n" );
687
688	#endif
689
690
691
692	// for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
693
694	if( constraint )
695	{
696	fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
697	// msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
698	// fprintf( stderr, "localhomshrink =\n" );
699	// outlocalhompt( localhomshrink, clus1, clus2 );
700	// weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
701	// fprintf( stderr, "after weight =\n" );
702	// outlocalhompt( localhomshrink, clus1, clus2 );
703	}
704	if( rnakozo && rnaprediction == 'm' )
705	{
706	makegrouprna( grouprna1, singlerna, topol[l][0] );
707	makegrouprna( grouprna2, singlerna, topol[l][1] );
708	}
709
710
711	/*
712	fprintf( stderr, "before align all\n" );
713	display( localcopy, njob );
714	fprintf( stderr, "\n" );
715	fprintf( stderr, "before align 1 %s \n", indication1 );
716	display( mseq1, clus1 );
717	fprintf( stderr, "\n" );
718	fprintf( stderr, "before align 2 %s \n", indication2 );
719	display( mseq2, clus2 );
720	fprintf( stderr, "\n" );
721	*/
722
723	if( !nevermemsave && ( constraint != 2 && alg != 'M' && ( len1 > 30000 \|\| len2 > 30000 ) ) )
724	{
725	fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
726	alg = 'M';
727	if( commonIP ) FreeIntMtx( commonIP );
728	commonIP = NULL;
729	commonAlloc1 = 0;
730	commonAlloc2 = 0;
731	}
732
733
734	// if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
735	if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
736	else ffttry = 0;
737	// ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
738	// fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
739	// fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
740	if( constraint == 2 )
741	{
742	if( alg == 'M' )
743	{
744	fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
745	exit( 1 );
746	}
747	fprintf( stderr, "c" );
748	if( alg == 'A' )
749	{
750	imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
751	if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
752	pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
753	}
754	else if( alg == 'H' )
755	{
756	imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
757	pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
758	}
759	else if( alg == 'Q' )
760	{
761	imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
762	if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
763	pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
764	}
765	else if( alg == 'R' )
766	{
767	imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
768	pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
769	}
770	}
771	else if( force_fft \|\| ( use_fft && ffttry ) )
772	{
773	fprintf( stderr, "f" );
774	if( alg == 'M' )
775	{
776	fprintf( stderr, "m" );
777	pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
778	}
779	else
780	pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
781	}
782	else
783	{
784	fprintf( stderr, "d" );
785	fftlog[m1] = 0;
786	switch( alg )
787	{
788	case( 'a' ):
789	pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
790	break;
791	case( 'M' ):
792	fprintf( stderr, "m" );
793	pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
794	break;
795	case( 'A' ):
796	pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
797	break;
798	case( 'Q' ):
799	pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
800	break;
801	case( 'R' ):
802	pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
803	break;
804	case( 'H' ):
805	pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
806	break;
807	default:
808	ErrorExit( "ERROR IN SOURCE FILE" );
809	}
810	}
811
812	nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
813
814	#if SCOREOUT
815	fprintf( stderr, "score = %10.2f\n", pscore );
816	#endif
817
818	/*
819	fprintf( stderr, "after align 1 %s \n", indication1 );
820	display( mseq1, clus1 );
821	fprintf( stderr, "\n" );
822	fprintf( stderr, "after align 2 %s \n", indication2 );
823	display( mseq2, clus2 );
824	fprintf( stderr, "\n" );
825	*/
826
827	// writePre( njob, name, nlen, localcopy, 0 );
828
829	if( disp ) display( localcopy, njob );
830
831	pthread_mutex_lock( targ->mutex );
832	dep[l].done = 1;
833	(*nrunpt)--;
834	pthread_cond_broadcast( targ->treecond );
835
836	// pthread_mutex_unlock( targ->mutex );
837	// pthread_mutex_lock( targ->mutex );
838
839	for( i=0; (j=topol[l][0][i])!=-1; i++ )
840	strcpy( aseq[j], localcopy[j] );
841	for( i=0; (j=topol[l][1][i])!=-1; i++ )
842	strcpy( aseq[j], localcopy[j] );
843	pthread_mutex_unlock( targ->mutex );
844
845	for( i=0; (j=topol[l][0][i])!=-1; i++ )
846	free( localcopy[j] );
847	for( i=0; (j=topol[l][1][i])!=-1; i++ )
848	free( localcopy[j] );
849	free( topol[l][0] );
850	free( topol[l][1] );
851	free( topol[l] );
852
853	}
854	}
855	#endif
856
857	void treebase( int nlen, char aseq, int nadd, char mergeoralign, char mseq1, char mseq2, int **topol, double effarr, int alloclen, LocalHom localhomtable, RNApair *singlerna, double effarr_kozo )
858	{
859	int i, l, m;
860	int len1nocommongap, len2nocommongap;
861	int len1, len2;
862	int clus1, clus2;
863	float pscore, tscore;
864	static char indication1, indication2;
865	static double *effarr1 = NULL;
866	static double *effarr2 = NULL;
867	static double *effarr1_kozo = NULL;
868	static double *effarr2_kozo = NULL;
869	static LocalHom ***localhomshrink = NULL;
870	static int *fftlog;
871	int m1, m2;
872	static int *gaplen;
873	static int *gapmap;
874	static int *alreadyaligned;
875	float dumfl = 0.0;
876	int ffttry;
877	RNApair *grouprna1 = NULL, *grouprna2 = NULL;
878
879	if( rnakozo && rnaprediction == 'm' )
880	{
881	grouprna1 = (RNApair *)calloc( njob, sizeof( RNApair ) );
882	grouprna2 = (RNApair *)calloc( njob, sizeof( RNApair ) );
883	}
884	else
885	{
886	grouprna1 = grouprna2 = NULL;
887	}
888
889	if( effarr1 == NULL )
890	{
891	fftlog = AllocateIntVec( njob );
892	effarr1 = AllocateDoubleVec( njob );
893	effarr2 = AllocateDoubleVec( njob );
894	indication1 = AllocateCharVec( 150 );
895	indication2 = AllocateCharVec( 150 );
896	gaplen = AllocateIntVec( *alloclen+10 );
897	gapmap = AllocateIntVec( *alloclen+10 );
898	alreadyaligned = AllocateIntVec( njob );
899	#if 0
900	#else
901	if( constraint )
902	{
903	localhomshrink = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
904	for( i=0; i<njob; i++ )
905	localhomshrink[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
906	}
907	#endif
908	effarr1_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
909	effarr2_kozo = AllocateDoubleVec( njob ); //tsuneni allocate sareru.
910	for( i=0; i<njob; i++ ) effarr1_kozo[i] = 0.0;
911	for( i=0; i<njob; i++ ) effarr2_kozo[i] = 0.0;
912	}
913
914	for( i=0; i<njob-nadd; i++ ) alreadyaligned[i] = 1;
915	for( i=njob-nadd; i<njob; i++ ) alreadyaligned[i] = 0;
916
917	for( l=0; l<njob; l++ ) fftlog[l] = 1;
918
919	#if 0
920	fprintf( stderr, "##### fftwinsize = %d, fftthreshold = %d\n", fftWinSize, fftThreshold );
921	#endif
922
923	#if 0
924	for( i=0; i<njob; i++ )
925	fprintf( stderr, "TBFAST effarr[%d] = %f\n", i, effarr[i] );
926	#endif
927
928
929	if( constraint )
930	calcimportance( njob, effarr, aseq, localhomtable );
931	// dontcalcimportance( njob, effarr, aseq, localhomtable ); // CHUUIII!!!!!
932
933
934	// writePre( njob, name, nlen, aseq, 0 );
935
936	tscore = 0.0;
937	for( l=0; l<njob-1; l++ )
938	{
939	if( mergeoralign[l] == 'n' )
940	{
941	// fprintf( stderr, "SKIP!\n" );
942	free( topol[l][0] );
943	free( topol[l][1] );
944	free( topol[l] );
945	continue;
946	}
947
948	m1 = topol[l][0][0];
949	m2 = topol[l][1][0];
950	len1 = strlen( aseq[m1] );
951	len2 = strlen( aseq[m2] );
952	if( *alloclen < len1 + len2 )
953	{
954	fprintf( stderr, "\nReallocating.." );
955	*alloclen = ( len1 + len2 ) + 1000;
956	ReallocateCharMtx( aseq, njob, *alloclen + 10 );
957	gaplen = realloc( gaplen, ( alloclen + 10 ) sizeof( int ) );
958	if( gaplen == NULL )
959	{
960	fprintf( stderr, "Cannot realloc gaplen\n" );
961	exit( 1 );
962	}
963	gapmap = realloc( gapmap, ( alloclen + 10 ) sizeof( int ) );
964	if( gapmap == NULL )
965	{
966	fprintf( stderr, "Cannot realloc gapmap\n" );
967	exit( 1 );
968	}
969	fprintf( stderr, "done. alloclen = %d\n", alloclen );
970	}
971
972	if( effarr_kozo )
973	{
974	clus1 = fastconjuction_noname_kozo( topol[l][0], aseq, mseq1, effarr1, effarr, effarr1_kozo, effarr_kozo, indication1 );
975	clus2 = fastconjuction_noname_kozo( topol[l][1], aseq, mseq2, effarr2, effarr, effarr2_kozo, effarr_kozo, indication2 );
976	}
977	else
978	{
979	clus1 = fastconjuction_noname( topol[l][0], aseq, mseq1, effarr1, effarr, indication1 );
980	clus2 = fastconjuction_noname( topol[l][1], aseq, mseq2, effarr2, effarr, indication2 );
981	}
982
983	if( mergeoralign[l] == '1' \|\| mergeoralign[l] == '2' ) // only in serial version
984	{
985	newgapstr = "=";
986	}
987	else
988	newgapstr = "-";
989
990
991	len1nocommongap = len1;
992	len2nocommongap = len2;
993	if( mergeoralign[l] == '1' ) // nai
994	{
995	findcommongaps( clus2, mseq2, gapmap );
996	commongappick( clus2, mseq2 );
997	len2nocommongap = strlen( mseq2[0] );
998	}
999	else if( mergeoralign[l] == '2' )
1000	{
1001	findcommongaps( clus1, mseq1, gapmap );
1002	commongappick( clus1, mseq1 );
1003	len1nocommongap = strlen( mseq1[0] );
1004	}
1005
1006
1007	fprintf( trap_g, "\nSTEP-%d\n", l );
1008	fprintf( trap_g, "group1 = %s\n", indication1 );
1009	fprintf( trap_g, "group2 = %s\n", indication2 );
1010
1011	#if 1
1012	fprintf( stderr, "\rSTEP % 5d /%d ", l+1, njob-1 );
1013	fflush( stderr );
1014	#else
1015	fprintf( stdout, "STEP %d /%d\n", l+1, njob-1 );
1016	fprintf( stderr, "STEP %d /%d\n", l+1, njob-1 );
1017	fprintf( stderr, "group1 = %.66s", indication1 );
1018	if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
1019	fprintf( stderr, "\n" );
1020	fprintf( stderr, "group2 = %.66s", indication2 );
1021	if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
1022	fprintf( stderr, "\n" );
1023	#endif
1024
1025
1026
1027	// for( i=0; i<clus1; i++ ) fprintf( stderr, "## STEP%d-eff for mseq1-%d %f\n", l+1, i, effarr1[i] );
1028
1029	if( constraint )
1030	{
1031	fastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
1032	// msfastshrinklocalhom( topol[l][0], topol[l][1], localhomtable, localhomshrink );
1033	// fprintf( stderr, "localhomshrink =\n" );
1034	// outlocalhompt( localhomshrink, clus1, clus2 );
1035	// weightimportance4( clus1, clus2, effarr1, effarr2, localhomshrink );
1036	// fprintf( stderr, "after weight =\n" );
1037	// outlocalhompt( localhomshrink, clus1, clus2 );
1038	}
1039	if( rnakozo && rnaprediction == 'm' )
1040	{
1041	makegrouprna( grouprna1, singlerna, topol[l][0] );
1042	makegrouprna( grouprna2, singlerna, topol[l][1] );
1043	}
1044
1045
1046	/*
1047	fprintf( stderr, "before align all\n" );
1048	display( aseq, njob );
1049	fprintf( stderr, "\n" );
1050	fprintf( stderr, "before align 1 %s \n", indication1 );
1051	display( mseq1, clus1 );
1052	fprintf( stderr, "\n" );
1053	fprintf( stderr, "before align 2 %s \n", indication2 );
1054	display( mseq2, clus2 );
1055	fprintf( stderr, "\n" );
1056	*/
1057
1058	if( !nevermemsave && ( constraint != 2 && alg != 'M' && ( len1 > 30000 \|\| len2 > 30000 ) ) )
1059	{
1060	fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode.\n", len1, len2 );
1061	alg = 'M';
1062	if( commonIP ) FreeIntMtx( commonIP );
1063	commonIP = NULL;
1064	commonAlloc1 = 0;
1065	commonAlloc2 = 0;
1066	}
1067
1068
1069	// if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
1070	if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 1000 && clus2 < 1000 );
1071	else ffttry = 0;
1072	// ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 && clus1 < 5000 && clus2 < 5000 ); // v6.708
1073	// fprintf( stderr, "f=%d, len1/fftlog[m1]=%f, clus1=%d, len2/fftlog[m2]=%f, clus2=%d\n", ffttry, (float)len1/fftlog[m1], clus1, (float)len2/fftlog[m2], clus2 );
1074	// fprintf( stderr, "f=%d, clus1=%d, fftlog[m1]=%d, clus2=%d, fftlog[m2]=%d\n", ffttry, clus1, fftlog[m1], clus2, fftlog[m2] );
1075	if( constraint == 2 )
1076	{
1077	if( alg == 'M' )
1078	{
1079	fprintf( stderr, "\n\nMemory saving mode is not supported.\n\n" );
1080	exit( 1 );
1081	}
1082	fprintf( stderr, "c" );
1083	if( alg == 'A' )
1084	{
1085	imp_match_init_strict( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, effarr1_kozo, effarr2_kozo, localhomshrink, 1 );
1086	if( rnakozo ) imp_rna( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1087	pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1088	}
1089	else if( alg == 'H' )
1090	{
1091	imp_match_init_strictH( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1092	pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1093	}
1094	else if( alg == 'Q' )
1095	{
1096	imp_match_init_strictQ( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1097	if( rnakozo ) imp_rnaQ( clus1, clus2, mseq1, mseq2, effarr1, effarr2, grouprna1, grouprna2, NULL, NULL, NULL );
1098	pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1099	}
1100	else if( alg == 'R' )
1101	{
1102	imp_match_init_strictR( NULL, clus1, clus2, strlen( mseq1[0] ), strlen( mseq2[0] ), mseq1, mseq2, effarr1, effarr2, localhomshrink, 1 );
1103	pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, localhomshrink, &dumfl, NULL, NULL, NULL, NULL );
1104	}
1105	}
1106	else if( force_fft \|\| ( use_fft && ffttry ) )
1107	{
1108	fprintf( stderr, "f" );
1109	if( alg == 'M' )
1110	{
1111	fprintf( stderr, "m" );
1112	pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
1113	}
1114	else
1115	pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
1116	}
1117	else
1118	{
1119	fprintf( stderr, "d" );
1120	fftlog[m1] = 0;
1121	switch( alg )
1122	{
1123	case( 'a' ):
1124	pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
1125	break;
1126	case( 'M' ):
1127	fprintf( stderr, "m" );
1128	pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1129	break;
1130	case( 'A' ):
1131	pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1132	break;
1133	case( 'Q' ):
1134	pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1135	break;
1136	case( 'R' ):
1137	pscore = R__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1138	break;
1139	case( 'H' ):
1140	pscore = H__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1141	break;
1142	default:
1143	ErrorExit( "ERROR IN SOURCE FILE" );
1144	}
1145	}
1146
1147	nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
1148
1149	#if SCOREOUT
1150	fprintf( stderr, "score = %10.2f\n", pscore );
1151	#endif
1152	tscore += pscore;
1153	/*
1154	fprintf( stderr, "after align 1 %s \n", indication1 );
1155	display( mseq1, clus1 );
1156	fprintf( stderr, "\n" );
1157	fprintf( stderr, "after align 2 %s \n", indication2 );
1158	display( mseq2, clus2 );
1159	fprintf( stderr, "\n" );
1160	*/
1161
1162	// writePre( njob, name, nlen, aseq, 0 );
1163
1164	if( disp ) display( aseq, njob );
1165
1166	if( mergeoralign[l] == '1' ) // jissainiha nai. atarashii hairetsu ha saigo dakara.
1167	{
1168	adjustgapmap( strlen( mseq2[0] )-len2nocommongap+len2, gapmap, mseq2[0] );
1169	restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen, '-' );
1170	findnewgaps( clus2, 0, mseq2, gaplen );
1171	insertnewgaps( njob, alreadyaligned, aseq, topol[l][1], topol[l][0], gaplen, gapmap, *alloclen, alg, '-' );
1172	for( i=0; i<njob; i++ ) eq2dash( aseq[i] );
1173	for( i=0; (m=topol[l][0][i])>-1; i++ ) alreadyaligned[m] = 1;
1174	}
1175	if( mergeoralign[l] == '2' )
1176	{
1177	// fprintf( stderr, ">mseq1[0] = \n%s\n", mseq1[0] );
1178	// fprintf( stderr, ">mseq2[0] = \n%s\n", mseq2[0] );
1179	adjustgapmap( strlen( mseq1[0] )-len1nocommongap+len1, gapmap, mseq1[0] );
1180	restorecommongaps( njob, aseq, topol[l][0], topol[l][1], gapmap, *alloclen, '-' );
1181	findnewgaps( clus1, 0, mseq1, gaplen );
1182	insertnewgaps( njob, alreadyaligned, aseq, topol[l][0], topol[l][1], gaplen, gapmap, *alloclen, alg, '-' );
1183	for( i=0; i<njob; i++ ) eq2dash( aseq[i] );
1184	for( i=0; (m=topol[l][1][i])>-1; i++ ) alreadyaligned[m] = 1;
1185	}
1186
1187	free( topol[l][0] );
1188	free( topol[l][1] );
1189	free( topol[l] );
1190	}
1191	#if SCOREOUT
1192	fprintf( stderr, "totalscore = %10.2f\n\n", tscore );
1193	#endif
1194	if( rnakozo && rnaprediction == 'm' )
1195	{
1196	if( grouprna1 ) free( grouprna1 ); // nakami ha?
1197	if( grouprna2 ) free( grouprna2 ); // nakami ha?
1198	grouprna1 = grouprna2 = NULL;
1199	}
1200	if( constraint )
1201	{
1202	if( localhomshrink ) // nen no tame
1203	{
1204	for( i=0; i<njob; i++ )
1205	{
1206	free( localhomshrink[i] );
1207	localhomshrink[i] = NULL;
1208	}
1209	free( localhomshrink );
1210	localhomshrink = NULL;
1211	}
1212	}
1213	}
1214
1215	static void WriteOptions( FILE *fp )
1216	{
1217
1218	if( dorp == 'd' ) fprintf( fp, "DNA\n" );
1219	else
1220	{
1221	if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
1222	else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
1223	else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
1224	}
1225	fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1226	if( use_fft ) fprintf( fp, "FFT on\n" );
1227
1228	fprintf( fp, "tree-base method\n" );
1229	if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
1230	else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
1231	if( tbitr \|\| tbweight )
1232	{
1233	fprintf( fp, "iterate at each step\n" );
1234	if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
1235	if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
1236	if( tbweight ) fprintf( fp, " weighted\n" );
1237	fprintf( fp, "\n" );
1238	}
1239
1240	fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1241
1242	if( alg == 'a' )
1243	fprintf( fp, "Algorithm A\n" );
1244	else if( alg == 'A' )
1245	fprintf( fp, "Algorithm A+\n" );
1246	else if( alg == 'C' )
1247	fprintf( fp, "Apgorithm A+/C\n" );
1248	else
1249	fprintf( fp, "Unknown algorithm\n" );
1250
1251	if( treemethod == 'X' )
1252	fprintf( fp, "Tree = UPGMA (mix).\n" );
1253	else if( treemethod == 'E' )
1254	fprintf( fp, "Tree = UPGMA (average).\n" );
1255	else if( treemethod == 'q' )
1256	fprintf( fp, "Tree = Minimum linkage.\n" );
1257	else
1258	fprintf( fp, "Unknown tree.\n" );
1259
1260	if( use_fft )
1261	{
1262	fprintf( fp, "FFT on\n" );
1263	if( dorp == 'd' )
1264	fprintf( fp, "Basis : 4 nucleotides\n" );
1265	else
1266	{
1267	if( fftscore )
1268	fprintf( fp, "Basis : Polarity and Volume\n" );
1269	else
1270	fprintf( fp, "Basis : 20 amino acids\n" );
1271	}
1272	fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
1273	fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
1274	}
1275	else
1276	fprintf( fp, "FFT off\n" );
1277	fflush( fp );
1278	}
1279
1280
1281	int main( int argc, char *argv[] )
1282	{
1283	static int *nlen;
1284	static float *selfscore;
1285	int nogaplen;
1286	static char name, seq;
1287	static char mseq1, mseq2;
1288	static char **bseq;
1289	static float iscore, iscore_kozo;
1290	static double eff, eff_kozo, *eff_kozo_mapped = NULL;
1291	int i, j, ien, ik, jk;
1292	static int *topol, *topol_kozo;
1293	static int *addmem;
1294	static Treedep *dep;
1295	static float len, len_kozo;
1296	FILE *prep;
1297	FILE *infp;
1298	FILE *orderfp;
1299	FILE *hat2p;
1300	double unweightedspscore;
1301	int alignmentlength;
1302	char *mergeoralign;
1303	int foundthebranch;
1304	int nsubalignments, maxmem;
1305	int **subtable;
1306	int *insubtable;
1307	int *preservegaps;
1308	char ***subalnpt;
1309
1310	char c;
1311	int alloclen;
1312	LocalHom **localhomtable = NULL;
1313	RNApair ***singlerna = NULL;
1314	float ssi, ssj, bunbo;
1315	static char *kozoarivec;
1316	int nkozo;
1317
1318	arguments( argc, argv );
1319	#ifndef enablemultithread
1320	nthread = 0;
1321	#endif
1322
1323	if( fastathreshold < 0.0001 ) constraint = 0;
1324
1325	if( inputfile )
1326	{
1327	infp = fopen( inputfile, "r" );
1328	if( !infp )
1329	{
1330	fprintf( stderr, "Cannot open %s\n", inputfile );
1331	exit( 1 );
1332	}
1333	}
1334	else
1335	infp = stdin;
1336
1337	getnumlen( infp );
1338	rewind( infp );
1339
1340
1341	nkozo = 0;
1342
1343	if( njob < 2 )
1344	{
1345	fprintf( stderr, "At least 2 sequences should be input!\n"
1346	"Only %d sequence found.\n", njob );
1347	exit( 1 );
1348	}
1349
1350	if( subalignment )
1351	{
1352	readsubalignmentstable( njob, NULL, NULL, &nsubalignments, &maxmem );
1353	fprintf( stderr, "nsubalignments = %d\n", nsubalignments );
1354	fprintf( stderr, "maxmem = %d\n", maxmem );
1355	subtable = AllocateIntMtx( nsubalignments, maxmem+1 );
1356	insubtable = AllocateIntVec( njob );
1357	for( i=0; i<njob; i++ ) insubtable[i] = 0;
1358	preservegaps = AllocateIntVec( njob );
1359	for( i=0; i<njob; i++ ) preservegaps[i] = 0;
1360	subalnpt = AllocateCharCub( nsubalignments, maxmem, 0 );
1361	readsubalignmentstable( njob, subtable, preservegaps, NULL, NULL );
1362	}
1363
1364	seq = AllocateCharMtx( njob, nlenmax+1 );
1365	mseq1 = AllocateCharMtx( njob, 0 );
1366	mseq2 = AllocateCharMtx( njob, 0 );
1367
1368	name = AllocateCharMtx( njob, B+1 );
1369	nlen = AllocateIntVec( njob );
1370	selfscore = AllocateFloatVec( njob );
1371
1372	topol = AllocateIntCub( njob, 2, 0 );
1373	len = AllocateFloatMtx( njob, 2 );
1374	iscore = AllocateFloatHalfMtx( njob );
1375	eff = AllocateDoubleVec( njob );
1376	kozoarivec = AllocateCharVec( njob );
1377
1378	mergeoralign = AllocateCharVec( njob );
1379
1380	dep = (Treedep *)calloc( njob, sizeof( Treedep ) );
1381	if( nadd ) addmem = AllocateIntVec( nadd+1 );
1382
1383	if( constraint )
1384	{
1385	localhomtable = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1386	for( i=0; i<njob; i++ )
1387	{
1388	localhomtable[i] = (LocalHom *)calloc( njob, sizeof( LocalHom ) );
1389	for( j=0; j<njob; j++ )
1390	{
1391	localhomtable[i][j].start1 = -1;
1392	localhomtable[i][j].end1 = -1;
1393	localhomtable[i][j].start2 = -1;
1394	localhomtable[i][j].end2 = -1;
1395	localhomtable[i][j].overlapaa = -1.0;
1396	localhomtable[i][j].opt = -1.0;
1397	localhomtable[i][j].importance = -1.0;
1398	localhomtable[i][j].next = NULL;
1399	localhomtable[i][j].korh = 'h';
1400	}
1401	}
1402
1403	fprintf( stderr, "Loading 'hat3' ... " );
1404	prep = fopen( "hat3", "r" );
1405	if( prep == NULL ) ErrorExit( "Make hat3." );
1406	readlocalhomtable( prep, njob, localhomtable, kozoarivec );
1407	fclose( prep );
1408	fprintf( stderr, "\ndone.\n" );
1409
1410
1411	nkozo = 0;
1412	for( i=0; i<njob; i++ )
1413	{
1414	// fprintf( stderr, "kozoarivec[%d] = %d\n", i, kozoarivec[i] );
1415	if( kozoarivec[i] ) nkozo++;
1416	}
1417	if( nkozo )
1418	{
1419	topol_kozo = AllocateIntCub( nkozo, 2, 0 );
1420	len_kozo = AllocateFloatMtx( nkozo, 2 );
1421	iscore_kozo = AllocateFloatHalfMtx( nkozo );
1422	eff_kozo = AllocateDoubleVec( nkozo );
1423	eff_kozo_mapped = AllocateDoubleVec( njob );
1424	}
1425
1426
1427	// outlocalhom( localhomtable, njob );
1428
1429	#if 0
1430	fprintf( stderr, "Extending localhom ... " );
1431	extendlocalhom2( njob, localhomtable );
1432	fprintf( stderr, "done.\n" );
1433	#endif
1434	}
1435
1436	#if 0
1437	readData( infp, name, nlen, seq );
1438	#else
1439	readData_pointer( infp, name, nlen, seq );
1440	fclose( infp );
1441	#endif
1442
1443	constants( njob, seq );
1444
1445	#if 0
1446	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
1447	#endif
1448
1449	initSignalSM();
1450
1451	initFiles();
1452
1453	WriteOptions( trap_g );
1454
1455	c = seqcheck( seq );
1456	if( c )
1457	{
1458	fprintf( stderr, "Illegal character %c\n", c );
1459	exit( 1 );
1460	}
1461
1462	// writePre( njob, name, nlen, seq, 0 );
1463
1464	if( treein )
1465	{
1466	#if 0
1467	if( nkozo )
1468	{
1469	fprintf( stderr, "Both structure and user tree have been given. Not yet supported!\n" );
1470	exit( 1 );
1471	}
1472	#endif
1473	fprintf( stderr, "Loading a tree ... " );
1474	loadtree( njob, topol, len, name, nlen, dep );
1475	fprintf( stderr, "\ndone.\n\n" );
1476	}
1477	else
1478	{
1479	if( tbutree == 0 )
1480	{
1481	for( i=1; i<njob; i++ )
1482	{
1483	if( nlen[i] != nlen[0] )
1484	{
1485	fprintf( stderr, "Input pre-aligned seqences or make hat2.\n" );
1486	exit( 1 );
1487	}
1488	}
1489
1490	fprintf( stderr, "Making a distance matrix .. \n" );
1491	fflush( stderr );
1492	ien = njob-1;
1493	for( i=0; i<njob; i++ )
1494	{
1495	selfscore[i] = naivepairscore11( seq[i], seq[i], penalty );
1496	}
1497	#ifdef enablemultithread
1498	if( nthread > 0 )
1499	{
1500	distancematrixthread_arg_t *targ;
1501	Jobtable jobpos;
1502	pthread_t *handle;
1503	pthread_mutex_t mutex;
1504
1505	jobpos.i = 0;
1506	jobpos.j = 0;
1507
1508	targ = calloc( nthread, sizeof( distancematrixthread_arg_t ) );
1509	handle = calloc( nthread, sizeof( pthread_t ) );
1510	pthread_mutex_init( &mutex, NULL );
1511
1512	for( i=0; i<nthread; i++ )
1513	{
1514	targ[i].thread_no = i;
1515	targ[i].njob = njob;
1516	targ[i].selfscore = selfscore;
1517	targ[i].iscore = iscore;
1518	targ[i].seq = seq;
1519	targ[i].jobpospt = &jobpos;
1520	targ[i].mutex = &mutex;
1521
1522	pthread_create( handle+i, NULL, distancematrixthread, (void *)(targ+i) );
1523	}
1524
1525	for( i=0; i<nthread; i++ )
1526	{
1527	pthread_join( handle[i], NULL );
1528	}
1529	pthread_mutex_destroy( &mutex );
1530	free( handle );
1531	free( targ );
1532	}
1533	else
1534	#endif
1535	{
1536	for( i=0; i<ien; i++ )
1537	{
1538	if( i % 10 == 0 )
1539	{
1540	fprintf( stderr, "\r% 5d / %d", i, ien );
1541	fflush( stderr );
1542	}
1543	ssi = selfscore[i];
1544	for( j=i+1; j<njob; j++ )
1545	{
1546	ssj = selfscore[j];
1547	bunbo = MIN( ssi, ssj );
1548	if( bunbo == 0.0 )
1549	iscore[i][j-i] = 1.0;
1550	else
1551	// iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / MIN( selfscore[i], selfscore[j] );
1552	iscore[i][j-i] = 1.0 - naivepairscore11( seq[i], seq[j], penalty ) / bunbo;
1553
1554	#if 0
1555	fprintf( stderr, "### ssj = %f\n", ssj );
1556	fprintf( stderr, "### selfscore[i] = %f\n", selfscore[i] );
1557	fprintf( stderr, "### selfscore[j] = %f\n", selfscore[j] );
1558	fprintf( stderr, "### rawscore = %f\n", naivepairscore11( seq[i], seq[j], penalty ) );
1559	#endif
1560	}
1561	}
1562	}
1563	fprintf( stderr, "\ndone.\n\n" );
1564	fflush( stderr );
1565	}
1566	else
1567	{
1568	if( multidist )
1569	{
1570	fprintf( stderr, "Loading 'hat2n' (aligned sequences - new sequences) ... " );
1571	prep = fopen( "hat2n", "r" );
1572	if( prep == NULL ) ErrorExit( "Make hat2." );
1573	readhat2_floathalf_pointer( prep, njob, name, iscore );
1574	fclose( prep );
1575	fprintf( stderr, "done.\n" );
1576
1577	fprintf( stderr, "Loading 'hat2i' (aligned sequences) ... " );
1578	prep = fopen( "hat2i", "r" );
1579	if( prep == NULL ) ErrorExit( "Make hat2i." );
1580	readhat2_floathalf_pointer( prep, njob-nadd, name, iscore );
1581	fclose( prep );
1582	fprintf( stderr, "done.\n" );
1583	}
1584	else
1585	{
1586	fprintf( stderr, "Loading 'hat2' ... " );
1587	prep = fopen( "hat2", "r" );
1588	if( prep == NULL ) ErrorExit( "Make hat2." );
1589	readhat2_floathalf_pointer( prep, njob, name, iscore );
1590	fclose( prep );
1591	fprintf( stderr, "done.\n" );
1592	}
1593	}
1594	#if 1
1595	if( distout )
1596	{
1597	hat2p = fopen( "hat2", "w" );
1598	WriteFloatHat2_pointer_halfmtx( hat2p, njob, name, iscore );
1599	fclose( hat2p );
1600	}
1601	#endif
1602	if( nkozo )
1603	{
1604	ien = njob-1;
1605	ik = 0;
1606	for( i=0; i<ien; i++ )
1607	{
1608	jk = ik+1;
1609	for( j=i+1; j<njob; j++ )
1610	{
1611	if( kozoarivec[i] && kozoarivec[j] )
1612	{
1613	iscore_kozo[ik][jk-ik] = iscore[i][j-i];
1614	}
1615	if( kozoarivec[j] ) jk++;
1616	}
1617	if( kozoarivec[i] ) ik++;
1618	}
1619	}
1620
1621	fprintf( stderr, "Constructing a UPGMA tree ... " );
1622	fflush( stderr );
1623	if( topin )
1624	{
1625	fprintf( stderr, "--topin has been disabled\n" );
1626	exit( 1 );
1627	// fprintf( stderr, "Loading a topology ... " );
1628	// loadtop( njob, iscore, topol, len );
1629	// fprintf( stderr, "\ndone.\n\n" );
1630	}
1631	else if( subalignment ) // merge error no tame
1632	{
1633	fixed_supg_float_realloc_nobk_halfmtx_treeout_constrained( njob, iscore, topol, len, name, nlen, dep, nsubalignments, subtable );
1634	}
1635	else if( treeout ) // merge error no tame
1636	{
1637	fixed_musclesupg_float_realloc_nobk_halfmtx_treeout( njob, iscore, topol, len, name, nlen, dep );
1638	}
1639	else
1640	{
1641	fixed_musclesupg_float_realloc_nobk_halfmtx( njob, iscore, topol, len, dep, 1 );
1642	}
1643	// else
1644	// ErrorExit( "Incorrect tree\n" );
1645
1646	if( nkozo )
1647	{
1648	// for( i=0; i<nkozo-1; i++ )
1649	// for( j=i+1; j<nkozo; j++ )
1650	// fprintf( stderr, "iscore_kozo[%d][%d] =~ %f\n", i, j, iscore_kozo[i][j-i] );
1651	fixed_musclesupg_float_realloc_nobk_halfmtx( nkozo, iscore_kozo, topol_kozo, len_kozo, NULL, 1 );
1652	}
1653	fprintf( stderr, "\ndone.\n\n" );
1654	fflush( stderr );
1655	}
1656
1657
1658	orderfp = fopen( "order", "w" );
1659	if( !orderfp )
1660	{
1661	fprintf( stderr, "Cannot open 'order'\n" );
1662	exit( 1 );
1663	}
1664	for( i=0; (j=topol[njob-2][0][i])!=-1; i++ )
1665	{
1666	fprintf( orderfp, "%d\n", j );
1667	}
1668	for( i=0; (j=topol[njob-2][1][i])!=-1; i++ )
1669	{
1670	fprintf( orderfp, "%d\n", j );
1671	}
1672	fclose( orderfp );
1673
1674	if( treeout && noalign )
1675	{
1676	writeData_pointer( prep_g, njob, name, nlen, seq );
1677	fprintf( stderr, "\n" );
1678	SHOWVERSION;
1679	return( 0 );
1680	}
1681
1682	// countnode( njob, topol, node0 );
1683	if( tbrweight )
1684	{
1685	weight = 3;
1686	#if 0
1687	utree = 0; counteff( njob, topol, len, eff ); utree = 1;
1688	#else
1689	counteff_simple_float_nostatic( njob, topol, len, eff );
1690	for( i=njob-nadd; i<njob; i++ ) eff[i] /= (double)100;
1691	#if 0
1692	fprintf( stderr, "###### WEIGHT = \n" );
1693	for( i=0; i<njob; i++ )
1694	{
1695	fprintf( stderr, "w[%d] = %f\n", i, eff[i] );
1696	}
1697	exit( 1 );
1698	#endif
1699	if( nkozo )
1700	{
1701	// counteff_simple_float( nkozo, topol_kozo, len_kozo, eff_kozo ); // single weight nanode iranai
1702	for( i=0,j=0; i<njob; i++ )
1703	{
1704	if( kozoarivec[i] )
1705	{
1706	// eff_kozo_mapped[i] = eff_kozo[j]; //
1707	eff_kozo_mapped[i] = eff[i]; // single weight
1708	j++;
1709	}
1710	else
1711	eff_kozo_mapped[i] = 0.0;
1712	// fprintf( stderr, "eff_kozo_mapped[%d] = %f\n", i, eff_kozo_mapped[i] );
1713	// fprintf( stderr, " eff[%d] = %f\n", i, eff[i] );
1714	}
1715	}
1716
1717
1718	#endif
1719	}
1720	else
1721	{
1722	for( i=0; i<njob; i++ ) eff[i] = 1.0;
1723	if( nkozo )
1724	{
1725	for( i=0; i<njob; i++ )
1726	{
1727	if( kozoarivec[i] )
1728	eff_kozo_mapped[i] = 1.0;
1729	else
1730	eff_kozo_mapped[i] = 0.0;
1731	}
1732	}
1733	}
1734
1735	FreeFloatHalfMtx( iscore, njob );
1736	FreeFloatMtx( len );
1737
1738	alloclen = nlenmax*2+1; //chuui!
1739	bseq = AllocateCharMtx( njob, alloclen );
1740
1741
1742	if( nadd )
1743	{
1744	alignmentlength = strlen( seq[0] );
1745	for( i=0; i<njob-nadd; i++ )
1746	{
1747	if( alignmentlength != strlen( seq[i] ) )
1748	{
1749	fprintf( stderr, "#################################################################################\n" );
1750	fprintf( stderr, "# ERROR! #\n" );
1751	fprintf( stderr, "# The original%4d sequences must be aligned #\n", njob-nadd );
1752	fprintf( stderr, "#################################################################################\n" );
1753	exit( 1 );
1754	}
1755	}
1756	if( addprofile )
1757	{
1758	alignmentlength = strlen( seq[njob-nadd] );
1759	for( i=njob-nadd; i<njob; i++ )
1760	{
1761	if( alignmentlength != strlen( seq[i] ) )
1762	{
1763	fprintf( stderr, "###############################################################################\n" );
1764	fprintf( stderr, "# ERROR! #\n" );
1765	fprintf( stderr, "# The%4d additional sequences must be aligned #\n", nadd );
1766	fprintf( stderr, "# Otherwise, try the '--add' option, instead of '--addprofile' option. #\n" );
1767	fprintf( stderr, "###############################################################################\n" );
1768	exit( 1 );
1769	}
1770	}
1771	for( i=0; i<nadd; i++ ) addmem[i] = njob-nadd+i;
1772	addmem[nadd] = -1;
1773	foundthebranch = 0;
1774	for( i=0; i<njob-1; i++ )
1775	{
1776	if( samemember( topol[i][0], addmem ) ) // jissainiha nai
1777	{
1778	mergeoralign[i] = '1';
1779	foundthebranch = 1;
1780	}
1781	else if( samemember( topol[i][1], addmem ) )
1782	{
1783	mergeoralign[i] = '2';
1784	foundthebranch = 1;
1785	}
1786	else
1787	{
1788	mergeoralign[i] = 'n';
1789	}
1790	}
1791	if( !foundthebranch )
1792	{
1793	fprintf( stderr, "###############################################################################\n" );
1794	fprintf( stderr, "# ERROR! #\n" );
1795	fprintf( stderr, "# There is no appropriate position to add the%4d sequences in the guide tree.#\n", nadd );
1796	fprintf( stderr, "# Check whether the%4d sequences form a monophyletic cluster. #\n", nadd );
1797	fprintf( stderr, "# If not, try the '--add' option, instead of the '--addprofile' option. #\n" );
1798	fprintf( stderr, "############################################################################### \n" );
1799	exit( 1 );
1800	}
1801	commongappick( nadd, seq+njob-nadd );
1802	for( i=njob-nadd; i<njob; i++ ) strcpy( bseq[i], seq[i] );
1803	}
1804	else
1805	{
1806	for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'n';
1807	for( j=njob-nadd; j<njob; j++ )
1808	{
1809	addmem[0] = j;
1810	addmem[1] = -1;
1811	for( i=0; i<njob-1; i++ )
1812	{
1813	if( samemember( topol[i][0], addmem ) ) // arieru
1814	{
1815	// fprintf( stderr, "HIT!\n" );
1816	if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1817	else mergeoralign[i] = '1';
1818	}
1819	else if( samemember( topol[i][1], addmem ) )
1820	{
1821	// fprintf( stderr, "HIT!\n" );
1822	if( mergeoralign[i] != 'n' ) mergeoralign[i] = 'w';
1823	else mergeoralign[i] = '2';
1824	}
1825	}
1826	}
1827
1828	for( i=0; i<nadd; i++ ) addmem[i] = njob-nadd+i;
1829	addmem[nadd] = -1;
1830	for( i=0; i<njob-1; i++ )
1831	{
1832	if( includemember( topol[i][0], addmem ) && includemember( topol[i][1], addmem ) )
1833	{
1834	mergeoralign[i] = 'w';
1835	}
1836	else if( includemember( topol[i][0], addmem ) )
1837	{
1838	mergeoralign[i] = '1';
1839	}
1840	else if( includemember( topol[i][1], addmem ) )
1841	{
1842	mergeoralign[i] = '2';
1843	}
1844	}
1845	#if 0
1846	for( i=0; i<njob-1; i++ )
1847	{
1848	fprintf( stderr, "mem0 = " );
1849	for( j=0; topol[i][0][j]>-1; j++ ) fprintf( stderr, "%d ", topol[i][0][j] );
1850	fprintf( stderr, "\n" );
1851	fprintf( stderr, "mem1 = " );
1852	for( j=0; topol[i][1][j]>-1; j++ ) fprintf( stderr, "%d ", topol[i][1][j] );
1853	fprintf( stderr, "\n" );
1854	fprintf( stderr, "i=%d, mergeoralign[] = %c\n", i, mergeoralign[i] );
1855	}
1856	#endif
1857	for( i=njob-nadd; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1858	}
1859
1860	commongappick( njob-nadd, seq );
1861	for( i=0; i<njob-nadd; i++ ) strcpy( bseq[i], seq[i] );
1862	}
1863	//--------------- kokokara ----
1864	else if( subalignment )
1865	{
1866	for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'a';
1867	for( i=0; i<nsubalignments; i++ )
1868	{
1869	fprintf( stderr, "Checking subalignment %d:\n", i+1 );
1870	alignmentlength = strlen( seq[subtable[i][0]] );
1871	// for( j=0; subtable[i][j]!=-1; j++ )
1872	// fprintf( stderr, " %d. %-30.30s\n", subtable[i][j]+1, name[subtable[i][j]]+1 );
1873	for( j=0; subtable[i][j]!=-1; j++ )
1874	{
1875	if( subtable[i][j] >= njob )
1876	{
1877	fprintf( stderr, "No such sequence, %d.\n", subtable[i][j]+1 );
1878	exit( 1 );
1879	}
1880	if( alignmentlength != strlen( seq[subtable[i][j]] ) )
1881	{
1882	fprintf( stderr, "\n" );
1883	fprintf( stderr, "###############################################################################\n" );
1884	fprintf( stderr, "# ERROR!\n" );
1885	fprintf( stderr, "# Subalignment %d must be aligned.\n", i+1 );
1886	fprintf( stderr, "# Please check the alignment lengths of following sequences.\n" );
1887	fprintf( stderr, "#\n" );
1888	fprintf( stderr, "# %d. %-10.10s -> %d letters (including gaps)\n", subtable[i][0]+1, name[subtable[i][0]]+1, alignmentlength );
1889	fprintf( stderr, "# %d. %-10.10s -> %d letters (including gaps)\n", subtable[i][j]+1, name[subtable[i][j]]+1, (int)strlen( seq[subtable[i][j]] ) );
1890	fprintf( stderr, "#\n" );
1891	fprintf( stderr, "# See http://mafft.cbrc.jp/alignment/software/merge.html for details.\n" );
1892	if( subalignmentoffset )
1893	{
1894	fprintf( stderr, "#\n" );
1895	fprintf( stderr, "# You specified seed alignment(s) consisting of %d sequences.\n", subalignmentoffset );
1896	fprintf( stderr, "# In this case, the rule of numbering is:\n" );
1897	fprintf( stderr, "# The aligned seed sequences are numbered as 1 .. %d\n", subalignmentoffset );
1898	fprintf( stderr, "# The input sequences to be aligned are numbered as %d .. %d\n", subalignmentoffset+1, subalignmentoffset+njob );
1899	}
1900	fprintf( stderr, "###############################################################################\n" );
1901	fprintf( stderr, "\n" );
1902	exit( 1 );
1903	}
1904	insubtable[subtable[i][j]] = 1;
1905	}
1906	for( j=0; j<njob-1; j++ )
1907	{
1908	if( includemember( topol[j][0], subtable[i] ) && includemember( topol[j][1], subtable[i] ) )
1909	{
1910	mergeoralign[j] = 'n';
1911	}
1912	}
1913	foundthebranch = 0;
1914	for( j=0; j<njob-1; j++ )
1915	{
1916	if( samemember( topol[j][0], subtable[i] ) \|\| samemember( topol[j][1], subtable[i] ) )
1917	{
1918	foundthebranch = 1;
1919	fprintf( stderr, " -> OK\n" );
1920	break;
1921	}
1922	}
1923	if( !foundthebranch )
1924	{
1925	system( "cp infile.tree GuideTree" ); // tekitou
1926	fprintf( stderr, "\n" );
1927	fprintf( stderr, "###############################################################################\n" );
1928	fprintf( stderr, "# ERROR!\n" );
1929	fprintf( stderr, "# Subalignment %d does not form a monophyletic cluster\n", i+1 );
1930	fprintf( stderr, "# in the guide tree ('GuideTree' in this directory) internally computed.\n" );
1931	fprintf( stderr, "# If you really want to use this subalignment, pelase give a tree with --treein \n" );
1932	fprintf( stderr, "# http://mafft.cbrc.jp/alignment/software/treein.html\n" );
1933	fprintf( stderr, "# http://mafft.cbrc.jp/alignment/software/merge.html\n" );
1934	if( subalignmentoffset )
1935	{
1936	fprintf( stderr, "#\n" );
1937	fprintf( stderr, "# You specified seed alignment(s) consisting of %d sequences.\n", subalignmentoffset );
1938	fprintf( stderr, "# In this case, the rule of numbering is:\n" );
1939	fprintf( stderr, "# The aligned seed sequences are numbered as 1 .. %d\n", subalignmentoffset );
1940	fprintf( stderr, "# The input sequences to be aligned are numbered as %d .. %d\n", subalignmentoffset+1, subalignmentoffset+njob );
1941	}
1942	fprintf( stderr, "############################################################################### \n" );
1943	fprintf( stderr, "\n" );
1944	exit( 1 );
1945	}
1946	// commongappick( seq[subtable[i]], subalignment[i] ); // irukamo
1947	}
1948	#if 0
1949	for( i=0; i<njob-1; i++ )
1950	{
1951	fprintf( stderr, "STEP %d\n", i+1 );
1952	fprintf( stderr, "group1 = " );
1953	for( j=0; topol[i][0][j] != -1; j++ )
1954	fprintf( stderr, "%d ", topol[i][0][j]+1 );
1955	fprintf( stderr, "\n" );
1956	fprintf( stderr, "group2 = " );
1957	for( j=0; topol[i][1][j] != -1; j++ )
1958	fprintf( stderr, "%d ", topol[i][1][j]+1 );
1959	fprintf( stderr, "\n" );
1960	fprintf( stderr, "%d -> %c\n\n", i, mergeoralign[i] );
1961	}
1962	#endif
1963
1964	for( i=0; i<njob; i++ )
1965	{
1966	if( insubtable[i] ) strcpy( bseq[i], seq[i] );
1967	else gappick0( bseq[i], seq[i] );
1968	}
1969
1970	for( i=0; i<nsubalignments; i++ )
1971	{
1972	for( j=0; subtable[i][j]!=-1; j++ ) subalnpt[i][j] = bseq[subtable[i][j]];
1973	if( !preservegaps[i] ) commongappick( j, subalnpt[i] );
1974	}
1975
1976	FreeIntMtx( subtable );
1977	free( insubtable );
1978	for( i=0; i<nsubalignments; i++ ) free( subalnpt[i] );
1979	free( subalnpt );
1980	free( preservegaps );
1981	}
1982	//--------------- kokomade ----
1983	else
1984	{
1985	for( i=0; i<njob; i++ ) gappick0( bseq[i], seq[i] );
1986	for( i=0; i<njob-1; i++ ) mergeoralign[i] = 'a';
1987	}
1988
1989	if( rnakozo && rnaprediction == 'm' )
1990	{
1991	singlerna = (RNApair *)calloc( njob, sizeof( RNApair ) );
1992	prep = fopen( "hat4", "r" );
1993	if( prep == NULL ) ErrorExit( "Make hat4 using mccaskill." );
1994	fprintf( stderr, "Loading 'hat4' ... " );
1995	for( i=0; i<njob; i++ )
1996	{
1997	nogaplen = strlen( bseq[i] );
1998	singlerna[i] = (RNApair *)calloc( nogaplen+1, sizeof( RNApair ) );
1999	for( j=0; j<nogaplen; j++ )
2000	{
2001	singlerna[i][j] = (RNApair *)calloc( 1, sizeof( RNApair ) );
2002	singlerna[i][j][0].bestpos = -1;
2003	singlerna[i][j][0].bestscore = -1.0;
2004	}
2005	singlerna[i][nogaplen] = NULL;
2006	// fprintf( stderr, "### reading bpp %d ...\n", i );
2007	readmccaskill( prep, singlerna[i], nogaplen );
2008	}
2009	fclose( prep );
2010	fprintf( stderr, "\ndone.\n" );
2011	}
2012	else
2013	singlerna = NULL;
2014
2015
2016	fprintf( stderr, "Progressive alignment ... \n" );
2017
2018	#ifdef enablemultithread
2019	if( nthread > 0 && nadd == 0 )
2020	{
2021	treebasethread_arg_t *targ;
2022	int jobpos;
2023	pthread_t *handle;
2024	pthread_mutex_t mutex;
2025	pthread_cond_t treecond;
2026	int *fftlog;
2027	int nrun;
2028	int nthread_yoyu;
2029
2030	nthread_yoyu = nthread * 1;
2031	nrun = 0;
2032	jobpos = 0;
2033	targ = calloc( nthread_yoyu, sizeof( treebasethread_arg_t ) );
2034	fftlog = AllocateIntVec( njob );
2035	handle = calloc( nthread_yoyu, sizeof( pthread_t ) );
2036	pthread_mutex_init( &mutex, NULL );
2037	pthread_cond_init( &treecond, NULL );
2038
2039	for( i=0; i<njob; i++ ) dep[i].done = 0;
2040	for( i=0; i<njob; i++ ) fftlog[i] = 1;
2041
2042	if( constraint )
2043	calcimportance( njob, eff, bseq, localhomtable );
2044	// dontcalcimportance( njob, eff, bseq, localhomtable ); // CHUUUUIIII!!!
2045
2046	for( i=0; i<nthread_yoyu; i++ )
2047	{
2048	targ[i].thread_no = i;
2049	targ[i].nrunpt = &nrun;
2050	targ[i].njob = njob;
2051	targ[i].nlen = nlen;
2052	targ[i].jobpospt = &jobpos;
2053	targ[i].topol = topol;
2054	targ[i].dep = dep;
2055	targ[i].aseq = bseq;
2056	targ[i].effarr = eff;
2057	targ[i].alloclenpt = &alloclen;
2058	targ[i].localhomtable = localhomtable;
2059	targ[i].singlerna = singlerna;
2060	targ[i].effarr_kozo = eff_kozo_mapped;
2061	targ[i].fftlog = fftlog;
2062	targ[i].mergeoralign = mergeoralign;
2063	targ[i].mutex = &mutex;
2064	targ[i].treecond = &treecond;
2065
2066	pthread_create( handle+i, NULL, treebasethread, (void *)(targ+i) );
2067	}
2068
2069	for( i=0; i<nthread_yoyu; i++ )
2070	{
2071	pthread_join( handle[i], NULL );
2072	}
2073	pthread_mutex_destroy( &mutex );
2074	pthread_cond_destroy( &treecond );
2075	free( handle );
2076	free( targ );
2077	free( fftlog );
2078	}
2079	else
2080	#endif
2081	treebase( nlen, bseq, nadd, mergeoralign, mseq1, mseq2, topol, eff, &alloclen, localhomtable, singlerna, eff_kozo_mapped );
2082	fprintf( stderr, "\ndone.\n" );
2083	if( scoreout )
2084	{
2085	unweightedspscore = plainscore( njob, bseq );
2086	fprintf( stderr, "\nSCORE %s = %.0f, ", "(treebase)", unweightedspscore );
2087	fprintf( stderr, "SCORE / residue = %f", unweightedspscore / ( njob * strlen( bseq[0] ) ) );
2088	fprintf( stderr, "\n\n" );
2089	}
2090
2091	#if 0
2092	if( constraint )
2093	{
2094	LocalHom tmppt1, tmppt2;
2095	for( i=0; i<njob; i++ )
2096	{
2097	for( j=0; j<njob; j++ )
2098	{
2099	tmppt1 = localhomtable[i]+j;
2100	while( tmppt2 = tmppt1->next )
2101	{
2102	free( (void *)tmppt1 );
2103	tmppt1 = tmppt2;
2104	}
2105	free( (void *)tmppt1 );
2106	}
2107	free( (void *)(localhomtable[i]+j) );
2108	}
2109	free( (void *)localhomtable );
2110	}
2111	#endif
2112
2113	fprintf( trap_g, "done.\n" );
2114	fclose( trap_g );
2115	free( mergeoralign );
2116	if( rnakozo && rnaprediction == 'm' )
2117	{
2118	if( singlerna ) // nen no tame
2119	{
2120	for( i=0; i<njob; i++ )
2121	{
2122	for( j=0; singlerna[i][j]!=NULL; j++ )
2123	{
2124	if( singlerna[i][j] ) free( singlerna[i][j] );
2125	}
2126	if( singlerna[i] ) free( singlerna[i] );
2127	}
2128	free( singlerna );
2129	singlerna = NULL;
2130	}
2131	}
2132
2133	writeData_pointer( prep_g, njob, name, nlen, bseq );
2134	#if 0
2135	writeData( stdout, njob, name, nlen, bseq );
2136	writePre( njob, name, nlen, bseq, !contin );
2137	writeData_pointer( prep_g, njob, name, nlen, aseq );
2138	#endif
2139	#if IODEBUG
2140	fprintf( stderr, "OSHIMAI\n" );
2141	#endif
2142
2143	if( constraint ) FreeLocalHomTable( localhomtable, njob );
2144
2145	SHOWVERSION;
2146	return( 0 );
2147	}

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

Context Navigation

source: branches/items/GDE/MAFFT/mafft-7.055-with-extensions/core/tbfast.c

Download in other formats: