Context Navigation

disttbfast.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: 52.8 KB

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

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

Context Navigation

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

Download in other formats: