Context Navigation

splittbfast.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: 82.1 KB

Line
1	#include "mltaln.h"
2
3
4	#define TREE 1
5	#define PICKSIZE 50 // must be >= 3
6	#define WEIGHT 0
7	#define TOKYORIPARA 0.70 // 0.70
8	#define TOKYORIPARA_A 0.70 // changed
9	#define LENFAC 1
10	#define HUKINTOTREE 1
11	#define DIANA 0
12	#define MAX6DIST 10.0
13
14	// kouzoutai ni sasareru pointer ha static
15
16	#define DEBUG 0
17	#define IODEBUG 0
18	#define SCOREOUT 0
19
20	#define END_OF_VEC -1
21
22	static char *fastapath;
23	static int doalign;
24	static int fromaln;
25	static int uselongest;
26	static int treeout;
27	static int classsize;
28	static int picksize;
29	static int maxl;
30	static int tsize;
31	static int reorder;
32	static int pid;
33	static int maxdepth = 0;
34	static double tokyoripara;
35
36	static double lenfaca, lenfacb, lenfacc, lenfacd;
37	#define PLENFACA 0.01
38	#define PLENFACB 10000
39	#define PLENFACC 10000
40	#define PLENFACD 0.1
41	#define DLENFACA 0.01
42	#define DLENFACB 2500
43	#define DLENFACC 2500
44	#define DLENFACD 0.1
45
46	static char datafile[1000];
47	static char queryfile[1000];
48	static char resultfile[1000];
49
50	typedef struct _scores
51	{
52	double score;
53	int selfscore;
54	int orilen;
55	int *pointt;
56	int numinseq;
57	char *name;
58	// char *seq; // reallo
59	// char **seqpt;
60	int shimon;
61	} Scores;
62
63	int intcompare( const int a, const int b )
64	{
65	return( a - b );
66	}
67
68	int lcompare( const Scores a, const Scores b )
69	{
70	if( a->orilen < b->orilen ) return 1;
71	else if( a->orilen > b->orilen ) return -1;
72	else return 0;
73	}
74
75	int dcompare( const Scores a, const Scores b )
76	{
77	if( a->score > b->score ) return 1;
78	else if( a->score < b->score ) return -1;
79	else
80	{
81	if( a->selfscore < b->selfscore ) return 1;
82	else if( a->selfscore > b->selfscore ) return -1;
83	else
84	{
85	if( a->orilen < b->orilen ) return 1;
86	else if( a->orilen > b->orilen ) return -1;
87	else return 0;
88	}
89	}
90	}
91
92	#if 0
93	static void gappickandx0( char out, char in )
94	{
95	char c;
96	if( scoremtx == -1 )
97	{
98	while( *in )
99	{
100	if( (c=*in++) == '-' )
101	;
102	else if( c == 'u' )
103	*out++ = 't';
104	else if( amino_n[c] < 4 && amino_n[c] > -1 )
105	*out++ = c;
106	else
107	*out++ = 'n';
108	}
109	}
110	else
111	{
112	while( *in )
113	{
114	if( (c=*in++) == '-' )
115	;
116	else if( amino_n[c] < 20 && amino_n[c] > -1 )
117	*out++ = c;
118	else
119	*out++ = 'X';
120	}
121	}
122	*out = 0;
123	}
124
125	static int getkouho( int pickkouho, double prob, int nin, Scores scores, char **seq ) // 0 < prob < 1
126	{
127	int nkouho = 0;
128	int i, j;
129	int *iptr = pickkouho;
130	for( i=1; i<nin; i++ )
131	{
132	if( ( nkouho==0 \|\| rnd() < prob ) && ( scores[i].shimon != scores->shimon \|\| strcmp( seq[scores->numinseq], seq[scores[i].numinseq] ) ) )
133	{
134	#if 0
135	for( j=0; j<nkouho; j++ )
136	{
137	if( scores[i].shimon == scores[pickkouho[j]].shimon \|\| !strcmp( seq[scores[pickkouho[j]].numinseq], seq[scores[i].numinseq] ) )
138	break;
139	}
140	if( j == nkouho )
141	#endif
142	{
143	*iptr++ = i;
144	nkouho++;
145	// fprintf( stderr, "ok! nkouho=%d\n", nkouho );
146	}
147	}
148	else
149	{
150	;
151	// fprintf( stderr, "no! %d-%d\n", 0, scores[i].numinseq );
152	}
153	}
154	fprintf( stderr, "\ndone\n\n" );
155	return nkouho;
156	}
157
158	#endif
159
160	static void getfastascoremtx( int **tmpaminodis )
161	{
162	FILE *qfp;
163	FILE *dfp;
164	FILE *rfp;
165	int i, j;
166	char aa;
167	int slen;
168	int res;
169	char com[10000];
170	static char *tmpseq;
171	static char *tmpname;
172	double *resvec;
173
174	if( scoremtx == -1 )
175	{
176	tmpaminodis['a']['a'] = 5;
177	tmpaminodis['g']['g'] = 5;
178	tmpaminodis['c']['c'] = 5;
179	tmpaminodis['t']['t'] = 5;
180	tmpaminodis['n']['n'] = -1;
181
182	return;
183	}
184
185
186	tmpseq = calloc( 2000, sizeof( char ) );
187	tmpname = calloc( B, sizeof( char ) );
188	resvec = calloc( 1, sizeof( double ) );
189
190	// fprintf( stderr, "xformatting .. " );
191	dfp = fopen( datafile, "w" );
192	if( !dfp ) ErrorExit( "Cannot open datafile." );
193	sprintf( tmpname, ">+===========+%d ", 0 );
194	strcpy( tmpseq, "AAAAAAXXXXXX" );
195	strcat( tmpseq, "CCCCCCXXXXXX" );
196	strcat( tmpseq, "DDDDDDXXXXXX" );
197	strcat( tmpseq, "EEEEEEXXXXXX" );
198	strcat( tmpseq, "FFFFFFXXXXXX" );
199	strcat( tmpseq, "GGGGGGXXXXXX" );
200	strcat( tmpseq, "HHHHHHXXXXXX" );
201	strcat( tmpseq, "IIIIIIXXXXXX" );
202	strcat( tmpseq, "KKKKKKXXXXXX" );
203	strcat( tmpseq, "LLLLLLXXXXXX" );
204	strcat( tmpseq, "MMMMMMXXXXXX" );
205	strcat( tmpseq, "NNNNNNXXXXXX" );
206	strcat( tmpseq, "PPPPPPXXXXXX" );
207	strcat( tmpseq, "QQQQQQXXXXXX" );
208	strcat( tmpseq, "RRRRRRXXXXXX" );
209	strcat( tmpseq, "SSSSSSXXXXXX" );
210	strcat( tmpseq, "TTTTTTXXXXXX" );
211	strcat( tmpseq, "VVVVVVXXXXXX" );
212	strcat( tmpseq, "WWWWWWXXXXXX" );
213	strcat( tmpseq, "YYYYYYXXXXXX" );
214	slen = strlen( tmpseq );
215	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
216	fclose( dfp );
217	fprintf( stderr, "done.\n" );
218
219	for( i=0; i<20; i++ )
220	{
221	aa = amino[i];
222	// fprintf( stderr, "checking %c\n", aa );
223	*tmpseq = 0;
224	sprintf( tmpname, ">+===========+%d ", 0 );
225	for( j=0; j<6; j++ )
226	sprintf( tmpseq+strlen( tmpseq ), "%c", aa );
227	qfp = fopen( queryfile, "w" );
228	if( !qfp ) ErrorExit( "Cannot open queryfile." );
229	writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
230	fclose( qfp );
231
232	if( scoremtx == -1 )
233	sprintf( com, "%s -z3 -m10 -n -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, M, M, 0, queryfile, datafile, 6, resultfile );
234	else
235	sprintf( com, "%s -z3 -m10 -p -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, M, M, 0, queryfile, datafile, 2, resultfile );
236	res = system( com );
237	if( res )
238	{
239	fprintf( stderr, "error in %s", fastapath );
240	exit( 1 );
241	}
242
243	rfp = fopen( resultfile, "r" );
244	if( rfp == NULL )
245	ErrorExit( "file 'fasta.$$' does not exist\n" );
246	res = ReadFasta34m10_scoreonly( rfp, resvec, 1 );
247	fprintf( stderr, "%c: %f\n", 'A'+i, *resvec/6 );
248	fclose( rfp );
249	if( ( (int)*resvec % 6 ) > 0.0 )
250	{
251	fprintf( stderr, "Error in blast, resvec=%f\n", resvec );
252	fprintf( stderr, "Error in blast, resvec/6=%f\n", resvec/6 );
253	exit( 1 );
254	}
255	tmpaminodis[(int)aa][(int)aa] = (int)( *resvec / 6 );
256	// fprintf( stderr, "resvec=%f, tmpaminodis[aa][aa] = %d\n", resvec, tmpaminodis[aa][aa] );
257	}
258	tmpaminodis['X']['X'] = -1;
259	free( tmpname );
260	free( tmpseq );
261	free( resvec );
262	}
263
264	#if 0
265	static void getblastscoremtx( int **tmpaminodis )
266	{
267	FILE *qfp;
268	FILE *dfp;
269	FILE *rfp;
270	int i, j;
271	char aa;
272	int slen;
273	int res;
274	char com[10000];
275	static char *tmpseq;
276	static char *tmpname;
277	double *resvec;
278
279	if( scoremtx == -1 )
280	{
281	tmpaminodis['a']['a'] = 1;
282	tmpaminodis['g']['g'] = 1;
283	tmpaminodis['c']['c'] = 1;
284	tmpaminodis['t']['t'] = 1;
285
286	return;
287	}
288
289
290	tmpseq = calloc( 2000, sizeof( char ) );
291	tmpname = calloc( B, sizeof( char ) );
292	resvec = calloc( 1, sizeof( double ) );
293
294	// fprintf( stderr, "xformatting .. " );
295	dfp = fopen( datafile, "w" );
296	if( !dfp ) ErrorExit( "Cannot open datafile." );
297	sprintf( tmpname, "\0", i ); // BUG!!
298	strcpy( tmpseq, "AAAAAAXXXXXX" );
299	strcat( tmpseq, "CCCCCCXXXXXX" );
300	strcat( tmpseq, "DDDDDDXXXXXX" );
301	strcat( tmpseq, "EEEEEEXXXXXX" );
302	strcat( tmpseq, "FFFFFFXXXXXX" );
303	strcat( tmpseq, "GGGGGGXXXXXX" );
304	strcat( tmpseq, "HHHHHHXXXXXX" );
305	strcat( tmpseq, "IIIIIIXXXXXX" );
306	strcat( tmpseq, "KKKKKKXXXXXX" );
307	strcat( tmpseq, "LLLLLLXXXXXX" );
308	strcat( tmpseq, "MMMMMMXXXXXX" );
309	strcat( tmpseq, "NNNNNNXXXXXX" );
310	strcat( tmpseq, "PPPPPPXXXXXX" );
311	strcat( tmpseq, "QQQQQQXXXXXX" );
312	strcat( tmpseq, "RRRRRRXXXXXX" );
313	strcat( tmpseq, "SSSSSSXXXXXX" );
314	strcat( tmpseq, "TTTTTTXXXXXX" );
315	strcat( tmpseq, "VVVVVVXXXXXX" );
316	strcat( tmpseq, "WWWWWWXXXXXX" );
317	strcat( tmpseq, "YYYYYYXXXXXX" );
318	slen = strlen( tmpseq );
319	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
320	fclose( dfp );
321	if( scoremtx == -1 )
322	sprintf( com, "formatdb -p f -i %s -o F", datafile );
323	else
324	sprintf( com, "formatdb -i %s -o F", datafile );
325	system( com );
326	fprintf( stderr, "done.\n" );
327
328	for( i=0; i<20; i++ )
329	{
330	aa = amino[i];
331	fprintf( stderr, "checking %c\n", aa );
332	*tmpseq = 0;
333	for( j=0; j<6; j++ )
334	sprintf( tmpseq+strlen( tmpseq ), "%c", aa );
335	qfp = fopen( queryfile, "w" );
336	if( !qfp ) ErrorExit( "Cannot open queryfile." );
337	writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
338	fclose( qfp );
339
340	sprintf( com, "blastall -b %d -G 10 -E 1 -e 1e10 -p blastp -m 7 -i %s -d %s > %s\0", 1, queryfile, datafile, resultfile );
341	res = system( com );
342	if( res )
343	{
344	fprintf( stderr, "error in %s", "blastall" );
345	exit( 1 );
346	}
347
348	rfp = fopen( resultfile, "r" );
349	if( rfp == NULL )
350	ErrorExit( "file 'fasta.$$' does not exist\n" );
351	res = ReadBlastm7_scoreonly( rfp, resvec, 1 );
352	fprintf( stdout, "%c: %f\n", 'A'+i, *resvec/6 );
353	fclose( rfp );
354	if( ( (int)*resvec % 6 ) > 0.0 )
355	{
356	fprintf( stderr, "Error in blast, resvec=%f\n", resvec );
357	fprintf( stderr, "Error in blast, resvec/6=%f\n", resvec/6 );
358	exit( 1 );
359	}
360	tmpaminodis[aa][aa] = (int)( *resvec / 6 );
361	}
362	tmpaminodis['X']['X'] = 0;
363	free( tmpname );
364	free( tmpseq );
365	free( resvec );
366
367	}
368	#endif
369
370	static double callfasta( char seq, Scores scores, int nin, int *picks, int query, int rewritedata )
371	{
372	double *val;
373	FILE *qfp;
374	FILE *dfp;
375	FILE *rfp;
376	int i;
377	char com[10000];
378	static char datafile[1000];
379	static char queryfile[1000];
380	static char resultfile[1000];
381	static int pid;
382	static char *tmpseq;
383	static char *tmpname;
384	int slen;
385	int res;
386	static Scores *scoresbk = NULL;
387	static int ninbk = 0;
388
389	if( pid == 0 )
390	{
391	pid = (int)getpid();
392	sprintf( datafile, "/tmp/data-%d", pid );
393	sprintf( queryfile, "/tmp/query-%d", pid );
394	sprintf( resultfile, "/tmp/fasta-%d", pid );
395
396	tmpseq = calloc( nlenmax+1, sizeof( char ) );
397	tmpname = calloc( B+1, sizeof( char ) );
398	}
399
400	val = calloc( nin, sizeof( double ) );
401	// fprintf( stderr, "nin=%d, q=%d\n", nin, query );
402
403	if( rewritedata )
404	{
405	scoresbk = scores;
406	ninbk = nin;
407	// fprintf( stderr, "\nformatting .. " );
408	dfp = fopen( datafile, "w" );
409	if( !dfp ) ErrorExit( "Cannot open datafile." );
410	if( picks == NULL ) for( i=0; i<nin; i++ )
411	{
412	// fprintf( stderr, "i=%d / %d / %d\n", i, nin, njob );
413	// fprintf( stderr, "nlenmax = %d\n", nlenmax );
414	// fprintf( stderr, "scores[i].orilen = %d\n", scores[i].orilen );
415	// fprintf( stderr, "strlen( seq[scores[i].numinseq] = %d\n", strlen( seq[scores[i].numinseq] ) );
416	gappick0( tmpseq, seq[scores[i].numinseq] );
417	sprintf( tmpname, ">+===========+%d ", i );
418	slen = scores[i].orilen;
419	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
420	}
421	else for( i=0; i<nin; i++ )
422	{
423	gappick0( tmpseq, seq[scores[picks[i]].numinseq] );
424	sprintf( tmpname, ">+===========+%d ", i );
425	slen = scores[picks[i]].orilen;
426	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
427	}
428	fclose( dfp );
429	}
430
431
432	gappick0( tmpseq, seq[scores[query].numinseq] );
433	sprintf( tmpname, ">+==========+%d ", 0 );
434	slen = scores[query].orilen;
435	qfp = fopen( queryfile, "w" );
436	if( !qfp ) ErrorExit( "Cannot open queryfile." );
437	writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
438	fclose( qfp );
439
440	// fprintf( stderr, "calling fasta, nin=%d\n", nin );
441
442	if( scoremtx == -1 )
443	sprintf( com, "%s -z3 -m10 -n -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, nin, nin, 0, queryfile, datafile, 6, resultfile );
444	else
445	sprintf( com, "%s -z3 -m10 -p -Q -H -b%d -E%d -d%d %s %s %d > %s", fastapath, nin, nin, 0, queryfile, datafile, 2, resultfile );
446	res = system( com );
447	if( res )
448	{
449	fprintf( stderr, "error in %s", fastapath );
450	exit( 1 );
451	}
452	// fprintf( stderr, "fasta done\n" );
453
454	//exit( 1 );
455
456	rfp = fopen( resultfile, "r" );
457	if( rfp == NULL )
458	ErrorExit( "file 'fasta.$$' does not exist\n" );
459
460	// fprintf( stderr, "reading fasta\n" );
461	if( scoremtx == -1 )
462	res = ReadFasta34m10_scoreonly_nuc( rfp, val, nin );
463	else
464	res = ReadFasta34m10_scoreonly( rfp, val, nin );
465	// fprintf( stderr, "done. val[0] = %f\n", val[0] );
466
467
468	fclose( rfp );
469
470	#if 0
471	for( i=0; i<nin; i++ )
472	fprintf( stderr, "r[%d-%d] = %f\n", 0, i, val[i] );
473	exit( 1 );
474	#endif
475
476	return( val );
477	}
478	#if 0
479	static double callblast( char seq, Scores scores, int nin, int query, int rewritedata )
480	{
481	double *val;
482	FILE *qfp;
483	FILE *dfp;
484	FILE *rfp;
485	int i, j;
486	char com[10000];
487	static char datafile[1000];
488	static char queryfile[1000];
489	static char resultfile[1000];
490	static int pid;
491	static char *tmpseq;
492	static char *tmpname;
493	char *seqptr;
494	int slen;
495	int res;
496	static Scores *scoresbk = NULL;
497	static int ninbk = 0;
498
499	if( pid == 0 )
500	{
501	pid = (int)getpid();
502	sprintf( datafile, "/tmp/data-%d\0", pid );
503	sprintf( queryfile, "/tmp/query-%d\0", pid );
504	sprintf( resultfile, "/tmp/fasta-%d\0", pid );
505
506	tmpseq = calloc( nlenmax+1, sizeof( char ) );
507	tmpname = calloc( B+1, sizeof( char ) );
508	}
509
510	val = calloc( nin, sizeof( double ) );
511	// fprintf( stderr, "nin=%d, q=%d\n", nin, query );
512
513	if( rewritedata )
514	{
515	scoresbk = scores;
516	ninbk = nin;
517	fprintf( stderr, "\nformatting .. " );
518	dfp = fopen( datafile, "w" );
519	if( !dfp ) ErrorExit( "Cannot open datafile." );
520	for( i=0; i<nin; i++ )
521	{
522	// fprintf( stderr, "i=%d / %d / %d\n", i, nin, njob );
523	// fprintf( stderr, "nlenmax = %d\n", nlenmax );
524	// fprintf( stderr, "scores[i].orilen = %d\n", scores[i].orilen );
525	// fprintf( stderr, "strlen( seq[scores[i].numinseq] = %d\n", strlen( seq[scores[i].numinseq] ) );
526	gappick0( tmpseq, seq[scores[i].numinseq] );
527	sprintf( tmpname, "+===========+%d \0", i );
528	slen = scores[i].orilen;
529	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
530	}
531	fclose( dfp );
532
533	if( scoremtx == -1 )
534	sprintf( com, "formatdb -p f -i %s -o F", datafile );
535	else
536	sprintf( com, "formatdb -i %s -o F", datafile );
537	system( com );
538	// fprintf( stderr, "done.\n" );
539	}
540
541
542	gappick0( tmpseq, seq[scores[query].numinseq] );
543	sprintf( tmpname, "+==========+%d \0", 0 );
544	slen = scores[query].orilen;
545	qfp = fopen( queryfile, "w" );
546	if( !qfp ) ErrorExit( "Cannot open queryfile." );
547	writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
548	fclose( qfp );
549	// fprintf( stderr, "q=%s\n", tmpseq );
550
551	fprintf( stderr, "\ncalling blast .. \n" );
552	if( scoremtx == -1 )
553	sprintf( com, "blastall -b %d -e 1e10 -p blastn -m 7 -i %s -d %s > %s\0", nin, queryfile, datafile, resultfile );
554	else
555	sprintf( com, "blastall -b %d -G 10 -E 1 -e 1e10 -p blastp -m 7 -i %s -d %s > %s\0", nin, queryfile, datafile, resultfile );
556	res = system( com );
557	if( res ) ErrorExit( "error in blast" );
558
559	rfp = fopen( resultfile, "r" );
560	if( rfp == NULL )
561	ErrorExit( "file 'fasta.$$' does not exist\n" );
562	res = ReadBlastm7_scoreonly( rfp, val, nin );
563	fclose( rfp );
564
565	#if 0
566	for( i=0; i<nin; i++ )
567	fprintf( stderr, "r[%d-%d] = %f\n", 0, i, val[i] );
568	#endif
569
570	return( val );
571	}
572	#endif
573
574	#if 0
575	static void selhead( int *ar, int n )
576	{
577	int min = *ar;
578	int *minptr = ar;
579	int *ptr = ar;
580	int tmp;
581	n--;
582	ar++;
583	while( n-- )
584	{
585	if( ( tmp = *ptr++ ) < min )
586	{
587	min = tmp;
588	minptr = ptr;
589	}
590	}
591	if( minptr != ar )
592	{
593	tmp = *ar;
594	*ar = min;
595	*minptr = tmp;
596	}
597	return;
598	}
599	#endif
600
601	void arguments( int argc, char *argv[] )
602	{
603	int c;
604
605	doalign = 0;
606	fromaln = 0;
607	treeout = 0;
608	uselongest = 1;
609	reorder = 1;
610	nevermemsave = 0;
611	inputfile = NULL;
612	fftkeika = 0;
613	constraint = 0;
614	nblosum = 62;
615	fmodel = 0;
616	calledByXced = 0;
617	devide = 0;
618	use_fft = 0;
619	force_fft = 0;
620	fftscore = 1;
621	fftRepeatStop = 0;
622	fftNoAnchStop = 0;
623	weight = 3;
624	utree = 1;
625	tbutree = 1;
626	refine = 0;
627	check = 1;
628	cut = 0.0;
629	disp = 0;
630	outgap = 1;
631	alg = 'A';
632	mix = 0;
633	tbitr = 0;
634	scmtd = 5;
635	tbweight = 0;
636	tbrweight = 3;
637	checkC = 0;
638	treemethod = 'X';
639	contin = 0;
640	scoremtx = 1;
641	kobetsubunkatsu = 0;
642	dorp = NOTSPECIFIED;
643	ppenalty = -1530;
644	ppenalty_ex = NOTSPECIFIED;
645	poffset = -123;
646	kimuraR = NOTSPECIFIED;
647	pamN = NOTSPECIFIED;
648	geta2 = GETA2;
649	fftWinSize = NOTSPECIFIED;
650	fftThreshold = NOTSPECIFIED;
651	TMorJTT = JTT;
652	classsize = NOTSPECIFIED;
653	picksize = NOTSPECIFIED;
654	tokyoripara = NOTSPECIFIED;
655
656	while( --argc > 0 && (*++argv)[0] == '-' )
657	{
658	while ( ( c = *++argv[0] ) )
659	{
660	switch( c )
661	{
662	case 'p':
663	picksize = atoi( *++argv );
664	fprintf( stderr, "picksize = %d\n", picksize );
665	--argc;
666	goto nextoption;
667	case 's':
668	classsize = atoi( *++argv );
669	fprintf( stderr, "groupsize = %d\n", classsize );
670	--argc;
671	goto nextoption;
672	case 'i':
673	inputfile = *++argv;
674	fprintf( stderr, "inputfile = %s\n", inputfile );
675	--argc;
676	goto nextoption;
677	case 'f':
678	ppenalty = (int)( atof( ++argv ) 1000 - 0.5 );
679	// fprintf( stderr, "ppenalty = %d\n", ppenalty );
680	--argc;
681	goto nextoption;
682	case 'g':
683	ppenalty_ex = (int)( atof( ++argv ) 1000 - 0.5 );
684	fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
685	--argc;
686	goto nextoption;
687	case 'h':
688	poffset = (int)( atof( ++argv ) 1000 - 0.5 );
689	// fprintf( stderr, "poffset = %d\n", poffset );
690	--argc;
691	goto nextoption;
692	case 'k':
693	kimuraR = atoi( *++argv );
694	fprintf( stderr, "kimuraR = %d\n", kimuraR );
695	--argc;
696	goto nextoption;
697	case 'b':
698	nblosum = atoi( *++argv );
699	scoremtx = 1;
700	// fprintf( stderr, "blosum %d\n", nblosum );
701	--argc;
702	goto nextoption;
703	case 'j':
704	pamN = atoi( *++argv );
705	scoremtx = 0;
706	TMorJTT = JTT;
707	fprintf( stderr, "jtt %d\n", pamN );
708	--argc;
709	goto nextoption;
710	case 'm':
711	pamN = atoi( *++argv );
712	scoremtx = 0;
713	TMorJTT = TM;
714	fprintf( stderr, "tm %d\n", pamN );
715	--argc;
716	goto nextoption;
717	case 'T':
718	tokyoripara = (double)atof( *++argv );
719	--argc;
720	goto nextoption;
721	case 'l':
722	uselongest = 0;
723	break;
724	#if 1
725	case 'a':
726	fmodel = 1;
727	break;
728	#endif
729	case 'S':
730	doalign = 'f';
731	break;
732	case 'Z':
733	fromaln = 1;
734	break;
735	case 'L':
736	doalign = 1;
737	break;
738	case 'x':
739	reorder = 0;
740	break;
741	case 't':
742	treeout = 1;
743	break;
744	case 'r':
745	fmodel = -1;
746	break;
747	case 'D':
748	dorp = 'd';
749	break;
750	case 'P':
751	dorp = 'p';
752	break;
753	case 'e':
754	fftscore = 0;
755	break;
756	case 'O':
757	fftNoAnchStop = 1;
758	break;
759	#if 0
760	case 'R':
761	fftRepeatStop = 1;
762	break;
763	case 'Q':
764	calledByXced = 1;
765	break;
766	case 'a':
767	alg = 'a';
768	break;
769	#endif
770	case 'R':
771	alg = 'R';
772	break;
773	case 'Q':
774	alg = 'Q';
775	break;
776	case 'A':
777	alg = 'A';
778	break;
779	case 'N':
780	nevermemsave = 1;
781	break;
782	case 'M':
783	alg = 'M';
784	break;
785	case 'C':
786	alg = 'C';
787	break;
788	case 'F':
789	use_fft = 1;
790	break;
791	case 'G':
792	use_fft = 1;
793	force_fft = 1;
794	break;
795	case 'v':
796	tbrweight = 3;
797	break;
798	case 'd':
799	disp = 1;
800	break;
801	case 'o':
802	outgap = 0;
803	break;
804	case 'J':
805	tbutree = 0;
806	break;
807	case 'X':
808	treemethod = 'X'; // mix
809	break;
810	case 'E':
811	treemethod = 'E'; // upg (average)
812	break;
813	case 'q':
814	treemethod = 'q'; // minimum
815	break;
816	case 'z':
817	fftThreshold = atoi( *++argv );
818	--argc;
819	goto nextoption;
820	case 'w':
821	fftWinSize = atoi( *++argv );
822	--argc;
823	goto nextoption;
824	default:
825	fprintf( stderr, "illegal option %c\n", c );
826	argc = 0;
827	break;
828	}
829	}
830	nextoption:
831	;
832	}
833	if( argc == 1 )
834	{
835	cut = atof( (*argv) );
836	argc--;
837	}
838	if( argc != 0 )
839	{
840	fprintf( stderr, "options: Check source file !\n" );
841	exit( 1 );
842	}
843	if( tbitr == 1 && outgap == 0 )
844	{
845	fprintf( stderr, "conflicting options : o, m or u\n" );
846	exit( 1 );
847	}
848	if( alg == 'C' && outgap == 0 )
849	{
850	fprintf( stderr, "conflicting options : C, o\n" );
851	exit( 1 );
852	}
853	}
854
855	static int maxl;
856	static int tsize;
857	static int nunknown = 0;
858
859	int seq_grp_nuc( int grp, char seq )
860	{
861	int tmp;
862	int *grpbk = grp;
863	while( *seq )
864	{
865	tmp = amino_grp[(int)*seq++];
866	if( tmp < 4 )
867	*grp++ = tmp;
868	else
869	nunknown++;
870	}
871	*grp = END_OF_VEC;
872	return( grp-grpbk );
873	}
874
875	int seq_grp( int grp, char seq )
876	{
877	int tmp;
878	int *grpbk = grp;
879	while( *seq )
880	{
881	tmp = amino_grp[(int)*seq++];
882	if( tmp < 6 )
883	*grp++ = tmp;
884	else
885	nunknown++;
886	}
887	*grp = END_OF_VEC;
888	return( grp-grpbk );
889	}
890
891	void makecompositiontable_p( short table, int pointt )
892	{
893	int point;
894
895	while( ( point = *pointt++ ) != END_OF_VEC )
896	table[point]++;
897	}
898
899	int commonsextet_p( short table, int pointt )
900	{
901	int value = 0;
902	short tmp;
903	int point;
904	static short *memo = NULL;
905	static int *ct = NULL;
906	static int *cp;
907
908	if( !memo )
909	{
910	memo = (short *)calloc( tsize, sizeof( short ) );
911	if( !memo ) ErrorExit( "Cannot allocate memo\n" );
912	ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) );
913	if( !ct ) ErrorExit( "Cannot allocate memo\n" );
914	}
915
916	cp = ct;
917	while( ( point = *pointt++ ) != END_OF_VEC )
918	{
919	tmp = memo[point]++;
920	if( tmp < table[point] )
921	value++;
922	if( tmp == 0 ) *cp++ = point;
923	}
924	*cp = END_OF_VEC;
925
926	cp = ct;
927	while( *cp != END_OF_VEC )
928	memo[*cp++] = 0;
929
930	return( value );
931	}
932
933	void makepointtable_nuc( int pointt, int n )
934	{
935	int point;
936	register int *p;
937
938	p = n;
939	point = n++ 1024;
940	point += n++ 256;
941	point += n++ 64;
942	point += n++ 16;
943	point += n++ 4;
944	point += *n++;
945	*pointt++ = point;
946
947	while( *n != END_OF_VEC )
948	{
949	point -= p++ 1024;
950	point *= 4;
951	point += *n++;
952	*pointt++ = point;
953	}
954	*pointt = END_OF_VEC;
955	}
956
957	void makepointtable( int pointt, int n )
958	{
959	int point;
960	register int *p;
961
962	p = n;
963	point = n++ 7776;
964	point += n++ 1296;
965	point += n++ 216;
966	point += n++ 36;
967	point += n++ 6;
968	point += *n++;
969	*pointt++ = point;
970
971	while( *n != END_OF_VEC )
972	{
973	point -= p++ 7776;
974	point *= 6;
975	point += *n++;
976	*pointt++ = point;
977	}
978	*pointt = END_OF_VEC;
979	}
980
981	#if 1
982	static void pairalign( int nseq, int nlen, char seq, int mem1, int mem2, double weight, int *alloclen )
983	{
984	int l, len1, len2;
985	int clus1, clus2;
986	float pscore, tscore;
987	static int *fftlog;
988	static char indication1, indication2;
989	static double *effarr1 = NULL;
990	static double *effarr2 = NULL;
991	static char mseq1, mseq2;
992	float dumfl = 0.0;
993	int ffttry;
994	int m1, m2;
995	#if 0
996	int i, j;
997	#endif
998
999
1000	if( effarr1 == NULL )
1001	{
1002	fftlog = AllocateIntVec( nseq );
1003	effarr1 = AllocateDoubleVec( nseq );
1004	effarr2 = AllocateDoubleVec( nseq );
1005	indication1 = AllocateCharVec( 150 );
1006	indication2 = AllocateCharVec( 150 );
1007	mseq1 = AllocateCharMtx( nseq, 0 );
1008	mseq2 = AllocateCharMtx( nseq, 0 );
1009	for( l=0; l<nseq; l++ ) fftlog[l] = 1;
1010	}
1011
1012	tscore = 0.0;
1013	m1 = mem1[0];
1014	m2 = mem2[0];
1015	len1 = strlen( seq[m1] );
1016	len2 = strlen( seq[m2] );
1017	if( *alloclen < len1 + len2 )
1018	{
1019	fprintf( stderr, "\nReallocating.." );
1020	*alloclen = ( len1 + len2 ) + 1000;
1021	ReallocateCharMtx( seq, nseq, *alloclen + 10 );
1022	fprintf( stderr, "done. alloclen = %d\n", alloclen );
1023	}
1024
1025	#if WEIGHT
1026	clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1 );
1027	clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2 );
1028	#else
1029	clus1 = fastconjuction_noweight( mem1, seq, mseq1, effarr1, indication1 );
1030	clus2 = fastconjuction_noweight( mem2, seq, mseq2, effarr2, indication2 );
1031	#endif
1032
1033	#if 0
1034	for( i=0; i<clus1; i++ )
1035	fprintf( stderr, "in p seq[%d] = %s\n", mem1[i], seq[mem1[i]] );
1036	for( i=0; i<clus2; i++ )
1037	fprintf( stderr, "in p seq[%d] = %s\n", mem2[i], seq[mem2[i]] );
1038	#endif
1039
1040	#if 0
1041	fprintf( stderr, "group1 = %.66s", indication1 );
1042	if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
1043	fprintf( stderr, "\n" );
1044	fprintf( stderr, "group2 = %.66s", indication2 );
1045	if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
1046	fprintf( stderr, "\n" );
1047	#endif
1048
1049	// fprintf( stdout, "mseq1 = %s\n", mseq1[0] );
1050	// fprintf( stdout, "mseq2 = %s\n", mseq2[0] );
1051
1052	if( !nevermemsave && ( alg != 'M' && ( len1 > 10000 \|\| len2 > 10000 ) ) )
1053	{
1054	fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode\n", len1, len2 );
1055	alg = 'M';
1056	if( commonIP ) FreeIntMtx( commonIP );
1057	commonIP = 0;
1058	commonAlloc1 = 0;
1059	commonAlloc2 = 0;
1060	}
1061
1062	if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
1063	else ffttry = 0;
1064
1065	if( force_fft \|\| ( use_fft && ffttry ) )
1066	{
1067	fprintf( stderr, "\bf" );
1068	if( alg == 'M' )
1069	{
1070	fprintf( stderr, "\bm" );
1071	// fprintf( stderr, "%d-%d", clus1, clus2 );
1072	pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
1073	}
1074	else
1075	{
1076	// fprintf( stderr, "%d-%d", clus1, clus2 );
1077	pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
1078	}
1079	}
1080	else
1081	{
1082	fprintf( stderr, "\bd" );
1083	fftlog[m1] = 0;
1084	switch( alg )
1085	{
1086	case( 'a' ):
1087	pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
1088	break;
1089	case( 'M' ):
1090	fprintf( stderr, "\bm" );
1091	// fprintf( stderr, "%d-%d", clus1, clus2 );
1092	pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1093	break;
1094	case( 'Q' ):
1095	if( clus1 == 1 && clus2 == 1 )
1096	{
1097	// fprintf( stderr, "%d-%d", clus1, clus2 );
1098	pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
1099	}
1100	else
1101	{
1102	// fprintf( stderr, "%d-%d", clus1, clus2 );
1103	pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
1104	}
1105	break;
1106	case( 'A' ):
1107	if( clus1 == 1 && clus2 == 1 )
1108	{
1109	// fprintf( stderr, "%d-%d", clus1, clus2 );
1110	pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
1111	}
1112	else
1113	{
1114	// fprintf( stderr, "%d-%d", clus1, clus2 );
1115	pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
1116	}
1117	break;
1118	default:
1119	ErrorExit( "ERROR IN SOURCE FILE" );
1120	}
1121	}
1122	#if SCOREOUT
1123	fprintf( stderr, "score = %10.2f\n", pscore );
1124	#endif
1125	nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
1126	return;
1127	}
1128	#endif
1129
1130	#if 0
1131	static void treebase( int nseq, int nlen, char aseq, double eff, int nalign, int **topol, int alloclen ) // topol
1132	{
1133	int i, l;
1134	int nlim;
1135	int clus1, clus2;
1136
1137	nlim = nalign-1;
1138	for( l=0; l<nlim; l++ )
1139	{
1140	fprintf( stderr, "in treebase, l = %d\n", l );
1141	fprintf( stderr, "aseq[0] = %s\n", aseq[0] );
1142	fprintf( stderr, "aseq[topol[l][0][0]] = %s\n", aseq[topol[l][0][0]] );
1143	pairalign( nseq, nlen, aseq, topol[l][0], topol[l][1], eff, alloclen );
1144	free( topol[l][0] );
1145	free( topol[l][1] );
1146	}
1147	}
1148	#endif
1149
1150	static void WriteOptions( FILE *fp )
1151	{
1152
1153	if( dorp == 'd' ) fprintf( fp, "DNA\n" );
1154	else
1155	{
1156	if ( scoremtx == 0 ) fprintf( fp, "JTT %dPAM\n", pamN );
1157	else if( scoremtx == 1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
1158	else if( scoremtx == 2 ) fprintf( fp, "M-Y\n" );
1159	}
1160	fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1161	if( use_fft ) fprintf( fp, "FFT on\n" );
1162
1163	fprintf( fp, "tree-base method\n" );
1164	if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
1165	else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
1166	if( tbitr \|\| tbweight )
1167	{
1168	fprintf( fp, "iterate at each step\n" );
1169	if( tbitr && tbrweight == 0 ) fprintf( fp, " unweighted\n" );
1170	if( tbitr && tbrweight == 3 ) fprintf( fp, " reversely weighted\n" );
1171	if( tbweight ) fprintf( fp, " weighted\n" );
1172	fprintf( fp, "\n" );
1173	}
1174
1175	fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
1176
1177	if( alg == 'a' )
1178	fprintf( fp, "Algorithm A\n" );
1179	else if( alg == 'A' )
1180	fprintf( fp, "Algorithm A+\n" );
1181	else if( alg == 'S' )
1182	fprintf( fp, "Apgorithm S\n" );
1183	else if( alg == 'C' )
1184	fprintf( fp, "Apgorithm A+/C\n" );
1185	else
1186	fprintf( fp, "Unknown algorithm\n" );
1187
1188	if( treemethod == 'X' )
1189	fprintf( fp, "Tree = UPGMA (mix).\n" );
1190	else if( treemethod == 'E' )
1191	fprintf( fp, "Tree = .UPGMA (average)\n" );
1192	else if( treemethod == 'q' )
1193	fprintf( fp, "Tree = Minimum linkage.\n" );
1194	else
1195	fprintf( fp, "Unknown tree.\n" );
1196
1197	if( use_fft )
1198	{
1199	fprintf( fp, "FFT on\n" );
1200	if( dorp == 'd' )
1201	fprintf( fp, "Basis : 4 nucleotides\n" );
1202	else
1203	{
1204	if( fftscore )
1205	fprintf( fp, "Basis : Polarity and Volume\n" );
1206	else
1207	fprintf( fp, "Basis : 20 amino acids\n" );
1208	}
1209	fprintf( fp, "Threshold of anchors = %d%%\n", fftThreshold );
1210	fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
1211	}
1212	else
1213	fprintf( fp, "FFT off\n" );
1214	fflush( fp );
1215	}
1216
1217	#if 1
1218	static int splitseq_mq( Scores scores, int nin, int nlen, char seq, char orialn, char *name, char inputfile, int uniform, char *tree, int alloclen, int order, int whichgroup, double weight, int depthpt, int qinoya )
1219	{
1220	int val;
1221	int ii, jj;
1222	int *mptr;
1223	int treelen = 0; // by Mathog, a guess
1224	static int groupid = 0;
1225	static int branchid = 0;
1226	int i, j;
1227	int selfscore0;
1228	double **dfromc;
1229	double **dfromcp;
1230	float **pickmtx;
1231	float **yukomtx;
1232	static short *table1;
1233	Scores *outs, ptr;
1234	int *numin;
1235	int *tsukau;
1236	int belongto;
1237	char **children = NULL; // by Mathog, a guess
1238	char *tmptree;
1239	static int *orderpos = NULL;
1240	int rn;
1241	int npick;
1242	int nyuko;
1243	int *picks;
1244	int *yukos;
1245	int *s_p_map;
1246	int *p_o_map;
1247	int *s_y_map;
1248	int *y_o_map;
1249	int *closeh;
1250	int nkouho;
1251	int *pickkouho;
1252	int *iptr;
1253	int *jptr;
1254	int aligned;
1255	int ***topol;
1256	int *treeorder;
1257	int picktmp;
1258	float **len;
1259	double minscore;
1260	// double *minscoreinpick;
1261	float *hanni;
1262	double lenfac;
1263	double longer;
1264	double shorter;
1265	static char **mseq1 = NULL;
1266	static char **mseq2 = NULL;
1267	double *blastresults = NULL; // by Mathog, a guess
1268	static int palloclen = 0;
1269	float maxdist;
1270
1271	if( orderpos == NULL )
1272	orderpos = order;
1273	if( palloclen == 0 )
1274	palloclen = alloclen 2;
1275	if( mseq1 == NULL && doalign == 1 )
1276	{
1277	mseq1 = AllocateCharMtx( 1, palloclen );
1278	mseq2 = AllocateCharMtx( 1, palloclen );
1279	}
1280
1281	if( nin == 0 )
1282	{
1283	#if TREE
1284	if( treeout )
1285	{
1286	tree = (char )calloc( 1, sizeof( char ) );
1287	**tree = 0;
1288	}
1289	#endif
1290	return 1;
1291	}
1292
1293	if( nin < 2 \|\| uniform == -1 ) // kokodato muda deha nai ga
1294	{
1295	fprintf( stderr, "\nLeaf %d / %d ", ++branchid, njob );
1296	#if 0
1297	outputfile = AllocateCharVec( strlen( inputfile ) + 100 );
1298	sprintf( outputfile, "%s-%d", inputfile, branchid );
1299	if( uniform > 0 )
1300	// sprintf( outputfile, "%su%d", outputfile, uniform );
1301	sprintf( outputfile + strlen(outputfile), "u%d", uniform );
1302	fprintf( stderr, "GROUP %d: %d member(s) (%d) %s\n", branchid, nin, scores[0].numinseq, outputfile );
1303	outfp = fopen( outputfile, "w" );
1304	free( outputfile );
1305	if( outfp == NULL )
1306	{
1307	fprintf( stderr, "Cannot open %s\n", outputfile );
1308	exit( 1 );
1309	}
1310	for( j=0; j<nin; j++ )
1311	fprintf( outfp, ">G%d %s\n%s\n", branchid, scores[j].name+1, seq[scores[j].numinseq] );
1312	fclose( outfp );
1313	#endif
1314
1315
1316	#if TREE
1317	if( treeout )
1318	{
1319	treelen = 0;
1320	tmptree = calloc( 100, sizeof( char ) );
1321	for( j=0; j<nin; j++ )
1322	{
1323	treelen += sprintf( tmptree, "%d", scores[j].numinseq+1 );
1324	}
1325	free( tmptree );
1326
1327	tree = (char )calloc( treelen + nin + 15, sizeof( char ) );
1328	**tree = '\n';
1329	if( nin > 1 )
1330	{
1331	(tree+1) = '(';
1332	(tree+2) = '\0';
1333	}
1334	else
1335	{
1336	(tree+1) = '\0';
1337	}
1338	for( j=0; j<nin-1; j++ )
1339	{
1340	sprintf( tree+strlen( tree ), "%d,", scores[j].numinseq+1 );
1341	}
1342	sprintf( tree+strlen( tree ), "%d", scores[j].numinseq+1 );
1343	if( nin > 1 ) strcat( *tree, ")\n" );
1344	else strcat( *tree, "\n" );
1345	// fprintf( stdout, "tree = %s\n", tree );
1346	}
1347
1348	#endif
1349	for( j=0; j<nin; j++ )
1350	{
1351	*orderpos++ = scores[j].numinseq;
1352	// fprintf( stderr, "*order = %d\n", scores[j].numinseq );
1353	}
1354
1355	return 1;
1356	}
1357
1358
1359
1360	if( uselongest )
1361	{
1362	i = nin;
1363	ptr = scores;
1364	selfscore0 = scores->selfscore;
1365	belongto = 0;
1366	while( i-- )
1367	{
1368	// fprintf( stderr, "ptr-scores=%d, numinseq = %d, score = %f\n", ptr-scores, ptr->numinseq+1, ptr->score );
1369	if( ptr->selfscore > selfscore0 )
1370	{
1371	selfscore0 = ptr->selfscore;
1372	belongto = ptr-scores;
1373	}
1374	ptr++;
1375	}
1376	#if 1
1377	if( belongto != 0 )
1378	{
1379	// fprintf( stderr, "swap %d %s\n<->\n%d %s\n", 0, scores->name, belongto, (scores+belongto)->name );
1380	ptr = calloc( 1, sizeof( Scores ) );
1381	*ptr = scores[belongto];
1382	scores[belongto] = *scores;
1383	scores = ptr;
1384	free( ptr );
1385	}
1386	#endif
1387	}
1388	else
1389	{
1390	qsort( scores, nin, sizeof( Scores ), (int (*)())lcompare );
1391	belongto = (int)( 0.5 * nin );
1392	// fprintf( stderr, "lengths = %d, %d, %d\n", scores->orilen, scores[belongto].orilen, scores[nin-1].orilen );
1393	if( belongto != 0 )
1394	{
1395	// fprintf( stderr, "swap %d %s\n<->\n%d %s\n", 0, scores->name, belongto, (scores+belongto)->name );
1396	ptr = calloc( 1, sizeof( Scores ) );
1397	*ptr = scores[belongto];
1398	scores[belongto] = *scores;
1399	scores = ptr;
1400	free( ptr );
1401	}
1402	}
1403
1404	if( qinoya != scores->numinseq )
1405	// if( 1 \|\| qinoya != scores->numinseq )
1406	{
1407	// fprintf( stdout, "### scores->numinseq = %d, qinoya=%d, depth=%d\n", scores->numinseq, qinoya, *depthpt );
1408
1409
1410	if( doalign )
1411	{
1412	if( doalign == 'f' )
1413	{
1414	blastresults = callfasta( seq, scores, nin, NULL, 0, 1 );
1415	if( scores->selfscore != (int)blastresults[0] )
1416	{
1417	fprintf( stderr, "\n\nWARNING1: selfscore\n" );
1418	fprintf( stderr, "scores->numinseq = %d\n", scores->numinseq+1 );
1419	fprintf( stderr, "scores->orilen = %d\n", scores->orilen );
1420	fprintf( stderr, "scores->selfscore = %d, but blastresults[0] = %f\n", scores->selfscore, blastresults[0] );
1421	// if( abs( scores->selfscore - (int)blastresults[0] ) > 2 )
1422	// exit( 1 );
1423	// scores->selfscore = (int)blastresults[0]; //iinoka?
1424
1425	// exit( 1 );
1426	}
1427	}
1428	else
1429	gappick0( mseq1[0], seq[scores->numinseq] );
1430	}
1431	else
1432	{
1433	table1 = (short *)calloc( tsize, sizeof( short ) );
1434	if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
1435	makecompositiontable_p( table1, scores[0].pointt );
1436	}
1437
1438	selfscore0 = scores[0].selfscore;
1439	for( i=0; i<nin; i++ )
1440	{
1441	if( scores->orilen > scores[i].orilen )
1442	{
1443	longer = (double)scores->orilen;
1444	shorter = (double)scores[i].orilen;
1445	}
1446	else
1447	{
1448	longer = (double)scores[i].orilen; // nai
1449	shorter = (double)scores->orilen; //nai
1450	}
1451
1452	#if LENFAC
1453	lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
1454	// lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
1455	// fprintf( stderr, "lenfac = %f l=%d,%d\n", lenfac,scores->orilen, scores[i].orilen );
1456	#else
1457	lenfac = 1.0;
1458	#endif
1459
1460	if( doalign )
1461	{
1462	if( doalign == 'f' )
1463	{
1464	scores[i].score = ( 1.0 - blastresults[i] / MIN( scores->selfscore, scores[i].selfscore ) ) * 1;
1465	if( scores[i].score < 0.0 ) scores[i].score = 0.0;
1466	}
1467	else
1468	{
1469	if( fromaln )
1470	{
1471	// scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
1472	scores[i].score = ( 1.0 - (double)naivepairscore11( orialn[scores[i].numinseq], orialn[scores->numinseq], penalty ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
1473	}
1474	else
1475	{
1476	if( *depthpt == 0 ) fprintf( stderr, "\r%d / %d ", i, nin );
1477	gappick0( mseq2[0], seq[scores[i].numinseq] );
1478	// fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
1479	scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
1480	// fprintf( stderr, "scores[i] = %f\n", scores[i].score );
1481	// fprintf( stderr, "m1=%s\n", seq[scores[0].numinseq] );
1482	// fprintf( stderr, "m2=%s\n", seq[scores[i].numinseq] );
1483	// fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
1484	}
1485	}
1486	}
1487	else
1488	{
1489	scores[i].score = ( 1.0 - (double)commonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
1490	if( scores[i].score > MAX6DIST ) scores[i].score = MAX6DIST;
1491	}
1492	// if( i ) fprintf( stderr, "%d-%d d %4.2f len %d %d\n", 1, i+1, scores[i].score, scores->orilen, scores[i].orilen );
1493	}
1494	if( doalign == 'f' ) free( blastresults );
1495	if( doalign == 0 ) free( table1 );
1496	//exit( 1 );
1497	}
1498
1499	// fprintf( stderr, "sorting .. " );
1500	qsort( scores, nin, sizeof( Scores ), (int (*)())dcompare );
1501	// fprintf( stderr, "done.\n" );
1502
1503
1504	maxdist = scores[nin-1].score;
1505	if( fromaln ) // kanzen itch ga misalign sareteiru kamoshirenai.
1506	{
1507	if( scores[0].shimon == scores[nin-1].shimon && !strcmp( seq[scores[0].numinseq], seq[scores[nin-1].numinseq] ) )
1508	{
1509	maxdist = 0.0;
1510	}
1511	// fprintf( stderr, "maxdist?? = %f, nin=%d, %d inori\n", scores[nin-1].score, nin, scores[nin-1].numinseq+1 );
1512	}
1513
1514	// fprintf( stderr, "maxdist? = %f, nin=%d\n", scores[nin-1].score, nin );
1515
1516	if( nin == 1 ) fprintf( stderr, "nin=1, scores[0].score = %f\n", scores[0].score );
1517
1518	// kokoni if( nin < 2 \|\| ... )
1519
1520	picks = AllocateIntVec( nin+1 );
1521	s_p_map = AllocateIntVec( nin+1 );
1522	s_y_map = AllocateIntVec( nin+1 );
1523	pickkouho = AllocateIntVec( nin+1 );
1524	closeh = AllocateIntVec( nin+1 );
1525
1526	// nkouho = getkouho( pickkouho, (picksize+100)/nin, nin, scores, seq );
1527	// nkouho = getkouho( pickkouho, 1.0, nin, scores, seq ); // zenbu
1528	// fprintf( stderr, "selecting kouhos phase 2\n" );
1529	// if( nkouho == 0 )
1530	// {
1531	// fprintf( stderr, "selecting kouhos, phase 2\n" );
1532	// nkouho = getkouho( pickkouho, 1.0, nin, scores, seq );
1533	// }
1534	// fprintf( stderr, "\ndone\n\n" );
1535	for( i=0; i<nin; i++ ) pickkouho[i] = i+1; nkouho = nin-1; // zenbu
1536
1537
1538
1539	iptr = picks;
1540	*iptr++ = 0;
1541	npick = 1;
1542	if( nkouho > 0 )
1543	{
1544	// fprintf( stderr, "pickkouho[0] = %d\n", pickkouho[0] );
1545	// fprintf( stderr, "pickkouho[nin-1] = %d\n", pickkouho[nin-1] );
1546	picktmp = pickkouho[nkouho-1];
1547	// fprintf( stderr, "\nMOST DISTANT kouho=%d, nin=%d, nkouho=%d\n", picktmp, nin, nkouho );
1548	nkouho--;
1549	if( ( scores[picktmp].shimon == scores[0].shimon ) && ( !strcmp( seq[scores[0].numinseq], seq[scores[picktmp].numinseq] ) ) )
1550	{
1551	// fprintf( stderr, "known, j=%d (%d inori)\n", 0, scores[picks[0]].numinseq );
1552	// fprintf( stderr, "%s\n%s\n", seq[scores[picktmp].numinseq], seq[scores[picks[0]].numinseq] );
1553	;
1554	}
1555	else
1556	{
1557	*iptr++ = picktmp;
1558	npick++;
1559	// fprintf( stderr, "ok, %dth pick = %d (%d inori)\n", npick, picktmp, scores[picktmp].numinseq );
1560	}
1561	}
1562	i = 1;
1563	while( npick<picksize && nkouho>0 )
1564	{
1565	if( i )
1566	{
1567	i = 0;
1568	rn = nkouho * 0.5;
1569	// fprintf( stderr, "rn = %d\n", rn );
1570	}
1571	else
1572	{
1573	rn = rnd() * (nkouho);
1574	}
1575	picktmp = pickkouho[rn];
1576	// fprintf( stderr, "rn=%d/%d (%d inori), kouho=%d, nin=%d, nkouho=%d\n", rn, nkouho, scores[pickkouho[rn]].numinseq, pickkouho[rn], nin, nkouho );
1577
1578	// fprintf( stderr, "#kouho before swap\n" );
1579	// for( i=0; i<nkouho; i++ ) fprintf( stderr, "%d ", pickkouho[i] ); fprintf( stderr, "\n" );
1580
1581	nkouho--;
1582	pickkouho[rn] = pickkouho[nkouho];
1583	#if 1
1584	// fprintf( stderr, "#kouho after swap\n" );
1585	// for( i=0; i<nkouho; i++ ) fprintf( stderr, "%d ", pickkouho[i] ); fprintf( stderr, "\n" );
1586	for( j=0; j<npick; j++ )
1587	{
1588	if( scores[picktmp].shimon == scores[picks[j]].shimon && !strcmp( seq[scores[picks[j]].numinseq], seq[scores[picktmp].numinseq] ) )
1589	break;
1590	}
1591	if( j == npick )
1592	#endif
1593	{
1594	// fprintf( stderr, "ok, %dth pick = %d (%d inori)\n", npick, picktmp, scores[picktmp].numinseq );
1595	npick++;
1596	*iptr++ = picktmp;
1597	}
1598	else
1599	{
1600	// fprintf( stderr, "known, j=%d (%d inori)\n", j, scores[picks[j]].numinseq );
1601	}
1602	}
1603	#if 0
1604	for( i=0; i<nin; i++ )
1605	{
1606	fprintf( stderr, "i=%d/%d, scores[%d].score = %f, inori=%d\n", i, nin, i, scores[i].score, scores[i].numinseq );
1607	}
1608	fprintf( stderr, "range:nin=%d scores[%d].score <= %f\n", nin, npick, scores[nin-1].score);
1609	for( i=0; i<npick; i++ )
1610	{
1611	fprintf( stderr, "i=%d/%d, scores[%d].score = %f, inori=%d\n", i, npick, picks[i], scores[picks[i]].score, scores[picks[i]].numinseq );
1612	}
1613	exit( 1 );
1614	#endif
1615
1616	// fprintf( stderr, "\nnkouho=%d, defaultq2 = %d (%d inori)\n", nkouho, picks[npick-1], scores[picks[npick-1]].numinseq );
1617
1618	qsort( picks, npick, sizeof( int ), (int (*)())intcompare );
1619
1620	// fprintf( stderr, "allocating..\n" );
1621
1622	// fprintf( stderr, "allocating outs, npick = %d\n", npick );
1623	numin = calloc( npick, sizeof( int ) );
1624	tsukau = calloc( npick, sizeof( int ) );
1625	outs = calloc( npick, sizeof( Scores * ) );
1626	for( i=0; i<npick; i++ ) outs[i] = NULL;
1627	topol = AllocateIntCub( npick, 2, 0 );
1628	treeorder = AllocateIntVec( npick + 1 );
1629	len = AllocateFloatMtx( npick, 2 );
1630	pickmtx = AllocateFloatHalfMtx( npick );
1631	// yukomtx = AllocateFloatHalfMtx( npick );
1632	// minscoreinpick = AllocateDoubleVec( npick );
1633	yukos = AllocateIntVec( npick );
1634	p_o_map = AllocateIntVec( npick+1 );
1635	y_o_map = AllocateIntVec( npick+1 );
1636	hanni = AllocateFloatVec( npick );
1637
1638	for( i=0; i<nin; i++ ) s_p_map[i] = -1;
1639	// fprintf( stderr, "npick = %d\n", npick );
1640	// fprintf( stderr, "picks =" );
1641	for( i=0; i<npick; i++ )
1642	{
1643	s_p_map[picks[i]] = i;
1644	p_o_map[i] = scores[picks[i]].numinseq;
1645	// fprintf( stderr, " %d (%dinori)\n", picks[i], scores[picks[i]].numinseq+1 );
1646	}
1647	// fprintf( stderr, "\n" );
1648
1649	#if 0
1650	fprintf( stderr, "p_o_map =" );
1651	for( i=0; i<npick; i++ )
1652	{
1653	fprintf( stderr, " %d", p_o_map[i]+1 );
1654	}
1655	fprintf( stderr, "\n" );
1656	fprintf( stderr, "picks =" );
1657	for( i=0; i<npick; i++ )
1658	{
1659	fprintf( stderr, " %d", picks[i] );
1660	}
1661	fprintf( stderr, "\n" );
1662	#endif
1663
1664	for( j=0; j<nin; j++ )
1665	{
1666	if( s_p_map[j] != -1 )
1667	{
1668	pickmtx[0][s_p_map[j]] = (float)scores[j].score;
1669	// fprintf( stderr, "pickmtx[0][%d] = %f\n", s_p_map[j], pickmtx[0][s_p_map[j]] );
1670	}
1671	}
1672
1673	for( j=1; j<npick; j++ )
1674	{
1675	if( doalign )
1676	{
1677	if( doalign == 'f' )
1678	{
1679	// blastresults = callfasta( seq, scores, npick-j+1, picks+j-1, picks[j], 1 );
1680	blastresults = callfasta( seq, scores, npick, picks, picks[j], (j==1) );
1681	if( scores[picks[j]].selfscore != (int)blastresults[j] )
1682	{
1683	fprintf( stderr, "\n\nWARNING2: selfscore j=%d/%d\n", j, npick );
1684	fprintf( stderr, "scores[picks[j]].numinseq = %d\n", scores[picks[j]].numinseq+1 );
1685	fprintf( stderr, "scores[picks[j]].orilen = %d\n", scores[picks[j]].orilen );
1686	fprintf( stderr, "scores[picks[j]].selfscore = %d, but blastresults[j] = %f\n", scores[picks[j]].selfscore, blastresults[j] );
1687	// if( abs( scores[picks[j]].selfscore - (int)blastresults[j] ) > 2 )
1688	// exit( 1 );
1689	// scores->selfscore = (int)blastresults[0]; //iinoka?
1690	}
1691	}
1692	else
1693	gappick0( mseq1[0], seq[scores[picks[j]].numinseq] );
1694	}
1695	else
1696	{
1697	table1 = (short *)calloc( tsize, sizeof( short ) );
1698	if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
1699	makecompositiontable_p( table1, scores[picks[j]].pointt );
1700	}
1701
1702	selfscore0 = scores[picks[j]].selfscore;
1703	pickmtx[j][0] = 0.0;
1704	for( i=j+1; i<npick; i++ )
1705	{
1706	if( scores[picks[j]].orilen > scores[picks[i]].orilen )
1707	{
1708	longer = (double)scores[picks[j]].orilen;
1709	shorter = (double)scores[picks[i]].orilen;
1710	}
1711	else
1712	{
1713	longer = (double)scores[picks[i]].orilen;
1714	shorter = (double)scores[picks[j]].orilen;
1715	}
1716
1717	#if LENFAC
1718	lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
1719	// lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
1720	// fprintf( stderr, "lenfac = %f l=%d,%d\n", lenfac,scores->orilen, scores[i].orilen );
1721	#else
1722	lenfac = 1.0;
1723	#endif
1724
1725	if( doalign )
1726	{
1727	if( doalign == 'f' )
1728	{
1729	pickmtx[j][i-j] = ( 1.0 - blastresults[i] / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
1730	if( pickmtx[j][i-j] < 0.0 ) pickmtx[j][i-j] = 0.0;
1731	}
1732	else
1733	{
1734	if( fromaln )
1735	{
1736	fprintf( stderr, "%d-%d/%d\r", j, i, npick );
1737	pickmtx[j][i-j] = ( 1.0 - (double)naivepairscore11( orialn[scores[picks[i]].numinseq], orialn[scores[picks[j]].numinseq], penalty ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
1738	}
1739	else
1740	{
1741	// fprintf( stderr, "\r%d / %d ", i, nin );
1742	gappick0( mseq2[0], seq[scores[picks[i]].numinseq] );
1743	pickmtx[j][i-j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
1744	// fprintf( stderr, "scores[picks[i]] = %f\n", scores[picks[i]].score );
1745	}
1746	}
1747	}
1748	else
1749	{
1750	pickmtx[j][i-j] = ( 1.0 - (double)commonsextet_p( table1, scores[picks[i]].pointt ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * lenfac;
1751	if( pickmtx[j][i-j] > MAX6DIST ) pickmtx[j][i-j] = MAX6DIST;
1752	}
1753
1754	}
1755	if( doalign == 'f' ) free( blastresults );
1756	if( doalign == 0 ) free( table1 );
1757	}
1758
1759	dfromcp = AllocateDoubleMtx( npick, nin );
1760	dfromc = AllocateDoubleMtx( npick, 0 );
1761
1762	for( i=0; i<npick; i++ ) for( j=0; j<nin; j++ )
1763	dfromcp[i][j] = -0.5;
1764	for( j=0; j<nin; j++ )
1765	{
1766	dfromcp[0][j] = ( scores[j].score );
1767	// fprintf( stderr, "j=%d, s_p_map[j]=%d\n", j, s_p_map[j] );
1768	}
1769
1770	for( i=0; i<npick; i++ ) for( j=i; j<npick; j++ )
1771	{
1772	dfromcp[i][picks[j]] = dfromcp[j][picks[i]] = pickmtx[i][j-i];
1773	}
1774
1775	#if 0
1776	fprintf( stderr, "pickmtx = \n" );
1777	for( i=0; i<npick; i++ )
1778	{
1779	for( j=i; j<npick; j++ )
1780	{
1781	fprintf( stderr, "pickmtx[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
1782	}
1783	}
1784	exit( 1 );
1785	#endif
1786
1787
1788
1789	// for( i=0; i<npick-1; i++ ) for( j=i; j<npick; j++ )
1790	// fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
1791
1792	for( i=0; i<npick; i++ ) tsukau[i] = 1;
1793	for( i=0; i<nin; i++ ) closeh[i] = -1;
1794	for( i=0; i<npick; i++ )
1795	{
1796	closeh[picks[i]] = picks[i];
1797	// fprintf( stderr, "i=%d/%d, picks[i]=%d, %d inori, closeh[%d] = %d \n", i, npick, picks[i], p_o_map[i]+1, picks[i], closeh[picks[i]] );
1798	}
1799	#if 0
1800	fprintf( stderr, "closeh = \n" );
1801	for( i=0; i<nin; i++ )
1802	{
1803	fprintf( stderr, "%d ", closeh[i] );
1804	}
1805	fprintf( stderr, "\n" );
1806	#endif
1807	#if DIANA
1808	for( i=0; i<npick-1; i++ ) for( j=i; j<npick; j++ )
1809	fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
1810	fprintf( stderr, "DIANA!!\n" );
1811	if( npick > 2 )
1812	{
1813	float avdist;
1814	float avdist1;
1815	float avdist2;
1816	float maxavdist;
1817	int splinter;
1818	int count;
1819	int dochokoho;
1820	splinter = 0;
1821	int *docholist;
1822	int *docholistbk;
1823	maxavdist = 0.0;
1824	for( i=0; i<npick; i++ )
1825	{
1826	avdist = 0.0;
1827	for( j=i+1; j<npick; j++ )
1828	{
1829	avdist += pickmtx[i][j-i];
1830	}
1831	for( j=0; j<i; j++ )
1832	{
1833	avdist += pickmtx[j][i-j];
1834	}
1835	avdist /= (npick-1);
1836	fprintf( stderr, "avdist[%d] = %f\n", p_o_map[i] + 1, avdist );
1837	if( maxavdist < avdist )
1838	{
1839	maxavdist = avdist;
1840	splinter = i;
1841	}
1842	}
1843	fprintf( stderr, "splinter = %d (%d inori), maxavdist = %f\n", splinter, p_o_map[splinter]+1, maxavdist );
1844
1845	docholist = AllocateIntVec( npick );
1846	docholistbk = AllocateIntVec( npick );
1847	for( i=0; i<npick; i++ ) docholist[i] = 0;
1848	docholist[splinter] = 1;
1849	while( 1 )
1850	{
1851	for( i=0; i<npick; i++ ) docholistbk[i] = docholist[i];
1852	for( dochokoho = 0; dochokoho<npick; dochokoho++ )
1853	{
1854	fprintf( stderr, "dochokoho=%d\n", dochokoho );
1855	if( docholist[dochokoho] ) continue;
1856	count = 0;
1857	avdist1 = 0.0;
1858	i=dochokoho;
1859	{
1860	for( j=i+1; j<npick; j++ )
1861	{
1862	if( docholist[j] \|\| j == dochokoho ) continue;
1863	avdist1 += pickmtx[i][j-i];
1864	count++;
1865	}
1866	for( j=0; j<i; j++ )
1867	{
1868	if( docholist[j] \|\| j == dochokoho ) continue;
1869	avdist1 += pickmtx[j][i-j];
1870	count++;
1871	}
1872	}
1873	if( count < 1 ) avdist1 = 0.0;
1874	else avdist1 /= (float)count;
1875	fprintf( stderr, "docho %d (%dinori), avdist1 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist1 );
1876
1877	count = 0;
1878	avdist2 = 0.0;
1879	i=dochokoho;
1880	{
1881	for( j=i+1; j<npick; j++ )
1882	{
1883	if( !docholist[j] \|\| j == dochokoho ) continue;
1884	avdist2 += pickmtx[i][j-i];
1885	count++;
1886	}
1887	for( j=0; j<i; j++ )
1888	{
1889	if( !docholist[j] \|\| j == dochokoho ) continue;
1890	avdist2 += pickmtx[j][i-j];
1891	count++;
1892	}
1893	}
1894	if( count < 1 ) avdist2 = 0.0;
1895	else avdist2 /= (float)count;
1896	fprintf( stderr, "docho %d (%dinori), avdist2 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist2 );
1897
1898	if( avdist2 < avdist1 )
1899	{
1900	docholist[dochokoho] = 1;
1901	hanni[dochokoho] = avdist2;
1902	}
1903	else
1904	{
1905	docholist[dochokoho] = 0;
1906	hanni[dochokoho] = avdist1;
1907	}
1908	fprintf( stderr, "avdist1=%f, avdist2=%f\n", avdist1, avdist2 );
1909
1910	}
1911	for( i=0; i<npick; i++ ) if( docholist[i] != docholistbk[i] ) break;
1912	if( i == npick ) break;
1913
1914	fprintf( stderr, "docholist = \n" );
1915	for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", docholist[i] );
1916	fprintf( stderr, "\n" );
1917	}
1918	fprintf( stderr, "docholist = \n" );
1919	for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", docholist[i] );
1920	fprintf( stderr, "\n" );
1921
1922	for( i=0; i<npick; i++ ) if( docholist[i] == 0 ) break;
1923	yukos[0] = picks[i];
1924	for( i=0; i<npick; i++ ) if( docholist[i] == 1 ) break;
1925	yukos[1] = picks[splinter];
1926
1927	for( i=0; i<npick; i++ )
1928	{
1929	if( docholist[i] == 0 ) closeh[picks[i]] = yukos[0];
1930	if( docholist[i] == 1 ) closeh[picks[i]] = yukos[1];
1931	}
1932	// for( i=0; i<npick; i++ ) closeh[picks[i]] = -1; // CHUUI !! iminai
1933	nyuko = 2;
1934	free( docholist );
1935	free( docholistbk );
1936	}
1937	else if( npick > 1 )
1938	{
1939	nyuko = 2;
1940	yukos[0] = picks[0]; yukos[1] = picks[1];
1941	closeh[picks[0]] = yukos[0];
1942	closeh[picks[1]] = yukos[1];
1943	}
1944	else
1945	{
1946	nyuko = 1;
1947	yukos[0] = picks[0];
1948	closeh[picks[0]] = yukos[0];
1949	}
1950	#elif HUKINTOTREE
1951	if( npick > 2 )
1952	{
1953	#if 0
1954	float avdist;
1955	float maxavdist;
1956	int count;
1957	int splinter;
1958	maxavdist = 0.0;
1959	splinter=0;
1960	for( i=0; i<npick; i++ )
1961	{
1962	avdist = 0.0;
1963	for( j=i+1; j<npick; j++ )
1964	{
1965	avdist += pickmtx[i][j-i];
1966	}
1967	for( j=0; j<i; j++ )
1968	{
1969	avdist += pickmtx[j][i-j];
1970	}
1971	avdist /= (npick-1);
1972	fprintf( stderr, "avdist[%d] = %f\n", p_o_map[i] + 1, avdist );
1973	if( maxavdist < avdist )
1974	{
1975	maxavdist = avdist;
1976	splinter = i;
1977	}
1978	}
1979	fprintf( stderr, "splinter = %d (%d inori), maxavdist = %f\n", splinter, p_o_map[splinter]+1, maxavdist );
1980	#endif
1981
1982
1983	// fprintf( stderr, "check kaishi =>, npick=%d members = \n", npick );
1984	// for( i=0; i<npick; i++ ) fprintf( stderr, "%d (%d)", p_o_map[i]+1, picks[i] );
1985	// fprintf( stderr, "\n" );
1986	for( i=0; i<npick-1; i++ )
1987	{
1988	if( tsukau[i] == 0 ) continue;
1989	for( j=i+1; j<npick; j++ )
1990	{
1991	// float kijun = maxdist * 1/(npick-2);
1992	// float kijun = maxavdist * tokyoripara;
1993	float kijun;
1994	kijun = maxdist * tokyoripara; // atode kakunin
1995	// fprintf( stderr, "%d-%d\n", i, j );
1996	// fprintf( stderr, "maxdist = %f\n", maxdist );
1997	// if( i==0 && j == 1 ) continue; // machigai!! CHUUI!!
1998	// if( maxdist == pickmtx[i][j-i] ) continue;
1999	if( tsukau[j] == 0 ) continue;
2000	// fprintf( stderr, "checking %d-%d (%d-%d) %f, kijun=%f\n", p_o_map[i]+1, p_o_map[j]+1, i, j, pickmtx[i][j-i], kijun );
2001	if( pickmtx[i][j-i] < kijun )
2002	{
2003	// fprintf( stderr, "dame!! %d => %d, because %f < %f\n", p_o_map[j]+1, p_o_map[i]+1, pickmtx[i][j-i], kijun );
2004	#if 0
2005	if( scores[picks[i]].orilen > scores[picks[j]].orilen )
2006	{
2007	fprintf( stderr, "%d => %d\n", p_o_map[j]+1, p_o_map[i]+1 );
2008	tsukau[j] = 0;
2009	}
2010	else
2011	{
2012	fprintf( stderr, "%d => %d\n", p_o_map[i]+1, p_o_map[j]+1 );
2013	tsukau[i] = 0;
2014	}
2015	if( 0 && j == npick-1 ) tsukau[i] = 0;
2016	else tsukau[j] = 0;
2017	fprintf( stderr, "tsukau[%d] = %d (%d inori)\n", j, tsukau[j], p_o_map[j]+1 );
2018	#else
2019	tsukau[j] = 0;
2020	closeh[picks[j]] = closeh[picks[i]];
2021	// fprintf( stderr, "%d => tsukawanai\n", j );
2022	#endif
2023	}
2024	}
2025	}
2026	}
2027	for( ii=0,i=0; i<npick; i++ )
2028	{
2029	if( tsukau[i] )
2030	{
2031	dfromc[ii] = dfromcp[i];
2032	ii++;
2033	}
2034	else
2035	{
2036	free( dfromcp[i] );
2037	dfromcp[i] = NULL;
2038	}
2039	}
2040	dfromc[ii] = NULL;
2041
2042	for( ii=0,i=0; i<npick; i++ )
2043	{
2044	if( tsukau[i] )
2045	{
2046	for( jj=ii,j=i; j<npick; j++ )
2047	{
2048	if( tsukau[j] )
2049	{
2050	pickmtx[ii][jj-ii] = pickmtx[i][j-i];
2051	jj++;
2052	}
2053	}
2054	ii++;
2055	}
2056	}
2057	for( ; ii<npick; ii++ )
2058	{
2059	free( pickmtx[ii] );
2060	pickmtx[ii] = NULL;
2061	}
2062
2063	for( ii=0,i=0; i<npick; i++ )
2064	{
2065	if( tsukau[i] )
2066	{
2067	yukos[ii++] = picks[i];
2068	}
2069	}
2070
2071
2072	nyuko = ii;
2073	yukomtx = pickmtx;
2074	pickmtx = NULL;
2075
2076	#endif
2077	#if 0
2078	for( i=0; i<npick; i++ ) for( j=i; j<npick; j++ )
2079	{
2080	if( tsukau[i] == 1 && tsukau[j] == 1 )
2081	fprintf( stderr, "dist[%d][%d] = %f (ok)\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
2082	else if( tsukau[i] == 0 && tsukau[j] == 0 )
2083	fprintf( stderr, "dist[%d][%d] = %f (xx)\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
2084	else
2085	fprintf( stderr, "%d-%d, okashii\n", p_o_map[i]+1, p_o_map[j]+1 );
2086	}
2087	#endif
2088	// FreeFloatHalfMtx( pickmtx, npick );
2089
2090
2091	for( i=0; i<nin; i++ ) s_y_map[i] = -1;
2092	// fprintf( stderr, "npick = %d\n", npick );
2093	// fprintf( stderr, "yukos =" );
2094	for( i=0; i<nyuko; i++ )
2095	{
2096	s_y_map[yukos[i]] = i;
2097	y_o_map[i] = scores[yukos[i]].numinseq;
2098	// fprintf( stderr, " %d\n", yukos[i] );
2099	}
2100	// fprintf( stderr, "\n" );
2101	#if 0
2102	for( i=0; i<nyuko; i++ )
2103	{
2104	fprintf( stderr, "y_o_map[%d] = %d\n", i, y_o_map[i]+1 );
2105	}
2106	for( i=0; i<nyuko; i++ ) for( j=i; j<nyuko; j++ )
2107	{
2108	fprintf( stderr, "yukodist[%d][%d] = %f (ok)\n", y_o_map[i]+1, y_o_map[j]+1, yukomtx[i][j-i] );
2109	}
2110	#endif
2111
2112	for( i=0; i<nin; i++ )
2113	{
2114	if( closeh[i] != -1 )
2115	{
2116	// fprintf( stderr, "closeh[%d,%dinori] = %d,%dinori\n", i, scores[i].numinseq+1, closeh[i], scores[closeh[i]].numinseq+1 );
2117	}
2118	}
2119
2120	#if 0
2121	for( i=0; i<nyuko; i++ )
2122	{
2123	minscoreinpick[i] = 99.9;
2124	for( j=i+1; j<nyuko; j++ )
2125	{
2126	if( minscoreinpick[i] > yukomtx[i][j-i] )
2127	minscoreinpick[i] = yukomtx[i][j-i];
2128	}
2129	for( j=0; j<i; j++ )
2130	{
2131	if( minscoreinpick[i] > yukomtx[j][i-j] )
2132	minscoreinpick[i] = yukomtx[j][i-j];
2133	}
2134	fprintf( stderr, "minscoreinpick[%d(%dinori)] = %f\n", i, y_o_map[i]+1, minscoreinpick[i] );
2135	}
2136	#endif
2137
2138
2139	#if TREE
2140	if( treeout )
2141	{
2142	children = calloc( nyuko+1, sizeof( char * ) );
2143	for( i=0; i<nyuko+1; i++ ) children[i] = NULL;
2144	}
2145	#endif
2146	// fprintf( stderr, "done..\n" );
2147
2148	// fprintf( stderr, "classifying, nyuko=%d \n", nyuko );
2149	if( nyuko == 1 )
2150	{
2151	if( npick != 1 )
2152	{
2153	fprintf( stderr, "okashii, nyuko = 1, shikashi npick = %d\n", npick );
2154	exit( 1 );
2155	}
2156	// fprintf( stderr, "### itchi suru hazu, nazenara scores[nin-1].score=%f, selfscores=%d,%d\n", scores[nin-1].score, scores[nin-1].selfscore, scores->selfscore );
2157	// fprintf( stderr, "seq[%d] = scores->seq = \n%s\n", scores->numinseq, seq[scores->numinseq] );
2158
2159	uniform = -1;
2160	for( j=0; j<nin; j++ )
2161	{
2162	belongto = 0;
2163	outs[belongto] = realloc( outs[belongto], sizeof( Scores ) * ( numin[belongto] + 1 ) );
2164	outs[belongto][numin[belongto]] = scores[j];
2165	numin[belongto]++;
2166	}
2167	}
2168	else
2169	{
2170
2171	#if 0
2172	fprintf( stderr, "yukos = \n" );
2173	for( i=0; i<nyuko; i++ ) fprintf( stderr, "%d ", y_o_map[i] + 1 );
2174	fprintf( stderr, "\n" );
2175	#endif
2176	fprintf( stderr, "\n\n%dx%d distance matrix\n", nyuko, nin );
2177
2178	for( i=1; i<nyuko; i++ )
2179	{
2180	fprintf( stderr, "%d / %d \r", i, nyuko );
2181
2182	if( doalign )
2183	{
2184	if( doalign == 'f' )
2185	{
2186	blastresults = callfasta( seq, scores, nin, NULL, yukos[i], (i==1) );
2187	}
2188	else
2189	gappick0( mseq1[0], seq[scores[yukos[i]].numinseq] );
2190	}
2191	else
2192	{
2193	table1 = (short *)calloc( tsize, sizeof( short ) );
2194	if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
2195	makecompositiontable_p( table1, scores[yukos[i]].pointt );
2196	}
2197
2198	selfscore0 = scores[yukos[i]].selfscore;
2199	for( j=0; j<nin; j++ )
2200	{
2201	if( scores[yukos[i]].orilen > scores[j].orilen )
2202	{
2203	longer = scores[yukos[i]].orilen;
2204	shorter = scores[j].orilen;
2205	}
2206	else
2207	{
2208	shorter = scores[yukos[i]].orilen;
2209	longer = scores[j].orilen;
2210	}
2211
2212	#if LENFAC
2213	// lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
2214	lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
2215	// lenfac = 1.0 / ( shorter / longer * LENFACD + LENFACB / ( longer + LENFACC ) + LENFACA );
2216	// fprintf( stderr, "lenfac = %f, l=%d, %d\n", lenfac, scores[yukos[i]].orilen, scores[j].orilen );
2217	#else
2218	lenfac = 1.0;
2219	#endif
2220	#if 0 // iihazu -> dame
2221	ii = s_y_map[j]; jj=s_y_map[yukos[i]];
2222	if( ii != -1 && jj != -1 )
2223	{
2224	if( dfromc[ii][yukos[jj]] != -0.5 )
2225	{
2226	dfromc[i][j] = dfromc[ii][yukos[jj]];
2227	}
2228	else
2229	{
2230	if( ii > jj )
2231	{
2232	kk = jj;
2233	jj = ii;
2234	ii = kk;
2235	}
2236	dfromc[ii][yukos[jj]] =
2237	dfromc[i][j] = yukomtx[ii][jj-ii];
2238	}
2239	}
2240	else
2241	#else
2242	if( dfromc[i][j] == -0.5 )
2243	#endif
2244	{
2245	if( doalign )
2246	{
2247	if( doalign == 'f' )
2248	{
2249	dfromc[i][j] =
2250	( 1.0 - blastresults[j] / MIN( selfscore0, scores[j].selfscore ) ) * 1;
2251	if( dfromc[i][j] < 0.0 ) dfromc[i][j] = 0.0;
2252	}
2253	else
2254	{
2255	if( fromaln )
2256	{
2257	dfromc[i][j] = ( 1.0 - (double)naivepairscore11( orialn[scores[j].numinseq], orialn[scores[yukos[i]].numinseq], penalty ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
2258	}
2259	else
2260	{
2261	gappick0( mseq2[0], seq[scores[j].numinseq] );
2262	dfromc[i][j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
2263	}
2264	}
2265	}
2266	else
2267	{
2268	dfromc[i][j] = ( 1.0 - (double)commonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
2269	if( dfromc[i][j] > MAX6DIST ) dfromc[i][j] = MAX6DIST;
2270	}
2271	}
2272	// fprintf( stderr, "i,j=%d,%d (%d,%d)/ %d,%d, dfromc[][]=%f \n", i, j, scores[yukos[i]].numinseq+1, scores[j].numinseq+1, nyuko, nin, dfromc[i][j] );
2273
2274	// if( i == 1 )
2275	// fprintf( stdout, "&&& dfromc[%d][%d] (%d,%d) = %f\n", i, j, p_o_map[i], scores[j].numinseq, dfromc[i][j] );
2276	}
2277	// fprintf( stderr, "i=%d, freeing\n", i );
2278	if( !doalign ) free( table1 );
2279	if( doalign && doalign == 'f' ) free( blastresults );
2280	}
2281	fprintf( stderr, " \r" );
2282
2283
2284
2285
2286	for( i=0; i<nyuko; i++ ) numin[i] = 0;
2287	// fprintf( stderr, "### itchi shinai hazu, nazenara scores[nin-1].score=%f, selfscores=%d,%d, len=%d,%d, nin=%d\n", scores[nin-1].score, scores[nin-1].selfscore, scores->selfscore, scores->orilen, scores[nin-1].orilen, nin );
2288	for( j=0; j<nin; j++ )
2289	{
2290	#if 0
2291	belongto = s_y_map[j];
2292	{
2293	fprintf( stderr, "belongto = %d (%dinori)\n", belongto, y_o_map[belongto]+1 );
2294	}
2295	if( belongto == -1 && closeh[j] != -1 )
2296	#endif
2297	#if 0
2298	if( closeh[j] != -1 )
2299	{
2300	belongto = s_y_map[closeh[j]];
2301	// if( belongto != -1 )
2302	// fprintf( stderr, "known, %d(%dinori)->%d(%dinori)\n", j, scores[j].numinseq+1, belongto, y_o_map[belongto]+1 );
2303	}
2304	else
2305	// if( belongto == -1 )
2306	#else
2307	belongto = s_y_map[j];
2308	if( belongto == -1 )
2309	#endif
2310	{
2311	belongto = 0; // default ha horyu
2312	minscore = dfromc[0][j];
2313	for( i=0; i<nyuko; i++ )
2314	{
2315	// fprintf( stderr, "checking %d/%d,%d/%d (%d-%d inori) minscore=%f, dfromc[0][j]=%f, dfromc[i][j]=%f\n", i, nyuko, j, nin, y_o_map[i], scores[j].numinseq, minscore, dfromc[0][j], dfromc[i][j] );
2316	if( scores[j].shimon == scores[yukos[i]].shimon && !strcmp( seq[scores[j].numinseq], seq[y_o_map[i]] ) )
2317	{
2318	// fprintf( stderr, "yuko-%d (%d in ori) to score-%d (%d inori) ha kanzen itch\n", i, y_o_map[i], j, scores[j].numinseq );
2319	belongto = i;
2320	break;
2321	}
2322	if( dfromc[i][j] < minscore )
2323	// if( dfromc[i][j] < minscore && minscore-dfromc[i][j] > ( minscoreinpick[yukos[i]] + minscoreinpick[j] ) * 1.0 )
2324	// if( rnd() < 0.5 ) // CHUUI !!!!!
2325	{
2326	// fprintf( stderr, "yuko-%d (%d in ori) to score-%d (%d inori) ha tikai, %f>%f\n", i, y_o_map[i]+1, j, scores[j].numinseq+1, minscore, dfromc[i][j] );
2327	minscore = dfromc[i][j];
2328	belongto = i;
2329	}
2330	}
2331	}
2332	#if 0
2333	if( dfromc[belongto][j] > minscoreinpick[belongto] )
2334	{
2335	fprintf( stderr, "dame, %f > %f\n", dfromc[belongto][j], minscoreinpick[belongto] );
2336	belongto = npick;
2337	}
2338	else
2339	fprintf( stderr, "ok, %f < %f\n", dfromc[belongto][j], minscoreinpick[belongto] );
2340	#endif
2341	// fprintf( stderr, "j=%d (%d inori) -> %d (%d inori) d=%f\n", j, scores[j].numinseq+1, belongto, y_o_map[belongto]+1, dfromc[belongto][j] );
2342	// fprintf( stderr, "numin = %d\n", numin[belongto] );
2343	outs[belongto] = realloc( outs[belongto], sizeof( Scores ) * ( numin[belongto] + 1 ) );
2344	outs[belongto][numin[belongto]] = scores[j];
2345	numin[belongto]++;
2346
2347	}
2348	free( dfromcp );
2349	FreeDoubleMtx( dfromc );
2350
2351	// fprintf( stderr, "##### npick = %d\n", npick );
2352	// fprintf( stderr, "##### nyuko = %d\n", nyuko );
2353
2354
2355	if( nyuko > 2 )
2356	{
2357	fprintf( stderr, "upgma " );
2358	// veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
2359	fixed_musclesupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len, NULL, 1 );
2360	fprintf( stderr, "\r \r" );
2361	}
2362	else
2363	{
2364	topol[0][0] = (int )realloc( topol[0][0], 2 sizeof( int ) );
2365	topol[0][1] = (int )realloc( topol[0][1], 2 sizeof( int ) );
2366	topol[0][0][0] = 0;
2367	topol[0][0][1] = -1;
2368	topol[0][1][0] = 1;
2369	topol[0][1][1] = -1;
2370	}
2371	FreeFloatHalfMtx( yukomtx, npick );
2372
2373	#if 0
2374	ii = nyuko-1;
2375	fprintf( stderr, "nyuko = %d, topol[][] = \n", nyuko );
2376	for( j=0; j<nyuko-1; j++ )
2377	{
2378	fprintf( stderr, "STEP%d \n", j );
2379	for( i=0; ; i++ )
2380	{
2381	fprintf( stderr, "%d ", ( topol[j][0][i] )+0 );
2382	if( topol[j][0][i] == -1 ) break;
2383	}
2384	fprintf( stderr, "\n" );
2385	for( i=0; ; i++ )
2386	{
2387	fprintf( stderr, "%d ", ( topol[j][1][i] )+0 );
2388	if( topol[j][1][i] == -1 ) break;
2389	}
2390	fprintf( stderr, "\n" );
2391	fprintf( stderr, "\n" );
2392	}
2393	exit( 1 );
2394	#endif
2395	jptr = treeorder;
2396	iptr = topol[nyuko-2][0]; while( iptr != -1 ) jptr++ = *iptr++;
2397	iptr = topol[nyuko-2][1]; while( iptr != -1 ) jptr++ = *iptr++;
2398	*jptr++ = -1;
2399	for( j=0; j<nyuko; j++ )
2400	{
2401	// fprintf( stderr, "treeorder[%d] = %d\n", j, treeorder[j] );
2402	if( treeorder[j] == -1 ) break;
2403	}
2404	}
2405
2406	aligned = 1;
2407	for( i=0; i<nyuko; i++ )
2408	{
2409	ii = treeorder[i];
2410	#if 0
2411	if( numin[ii] > 1 )
2412	{
2413	fprintf( stderr, "\ncalling a child, pick%d (%d inori): # of mem=%d\n", i, p_o_map[ii]+1, numin[ii] );
2414	for( j=0; j<numin[ii]; j++ )
2415	{
2416	fprintf( stderr, "%d ", outs[ii][j].numinseq+1 );
2417	}
2418	fprintf( stderr, "\n" );
2419	}
2420	#endif
2421	aligned *= splitseq_mq( outs[ii], numin[ii], nlen, seq, orialn, name, inputfile, uniform, children+ii, alloclen, order, whichgroup, weight, depthpt, scores->numinseq );
2422	}
2423
2424
2425	for( i=0; i<nyuko; i++ )
2426	{
2427	if( !numin[i] )
2428	{
2429	fprintf( stderr, "i=%d/%d, ERROR!\n", i, nyuko );
2430	for( j=0; j<nyuko; j++ )
2431	fprintf( stderr, "numin[%d] = %d (rep=%d inori)\n", j, numin[j], y_o_map[j] );
2432	exit( 1 );
2433	}
2434	}
2435
2436	#if TREE
2437	if( treeout )
2438	{
2439	treelen = 0;
2440	for( i=0; i<nyuko; i++ )
2441	treelen += strlen( children[i] );
2442	tree = calloc( treelen + nin 3, sizeof ( char ) );
2443	}
2444	#endif
2445
2446
2447	if( nin >= classsize \|\| !aligned )
2448	val = 0;
2449	else
2450	val = 1;
2451
2452	if( nyuko > 1 )
2453	{
2454	int mem1p, mem2p;
2455	int mem1size, mem2size;
2456	int v1 = 0, v2 = 0, v3 = 0;
2457	int nlim;
2458	int l;
2459	static int *mem1 = NULL;
2460	static int *mem2 = NULL;
2461	char **parttree = NULL; // by Mathog
2462
2463	#if TREE
2464	if( treeout )
2465	{
2466	parttree = (char *)calloc( nyuko, sizeof( char ) );
2467	for( i=0; i<nyuko; i++ )
2468	{
2469	// fprintf( stderr, "allocating parttree, size = %d\n", treelen + nin * 5 );
2470	parttree[i] = calloc( treelen + nin * 5, sizeof ( char ) );
2471	strcpy( parttree[i], children[i] );
2472	free( children[i] );
2473	}
2474	free( children );
2475	}
2476	#endif
2477	if( mem1 == NULL )
2478	{
2479	mem1 = AllocateIntVec( njob+1 );
2480	mem2 = AllocateIntVec( njob+1 );
2481	}
2482
2483	// veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
2484
2485	// counteff_simple_float( nyuko, topol, len, eff );
2486
2487
2488	nlim = nyuko-1;
2489	for( l=0; l<nlim; l++ )
2490	{
2491	mem1p = topol[l][0];
2492	mptr = mem1;
2493	mem1size = 0;
2494	while( *mem1p != -1 )
2495	{
2496	// fprintf( stderr, "mem1p = %d (%d inori), numin[]=%d\n", mem1p, p_o_map[mem1p], numin[mem1p] );
2497	i = numin[mem1p]; ptr = outs[(mem1p++)];
2498	mem1size += i;
2499	while( i-- )
2500	{
2501	*mptr++ = (ptr++)->numinseq;
2502	}
2503	}
2504	*mptr = -1;
2505
2506	mem2p = topol[l][1];
2507	mptr = mem2;
2508	mem2size = 0;
2509	while( *mem2p != -1 )
2510	{
2511	// fprintf( stderr, "mem2p = %d (%d inori), numin[]=%d\n", mem2p, p_o_map[mem2p], numin[mem2p] );
2512	i = numin[mem2p]; ptr = outs[(mem2p++)];
2513	mem2size += i;
2514	while( i-- )
2515	{
2516	*mptr++ = (ptr++)->numinseq;
2517	}
2518	}
2519	*mptr = -1;
2520
2521	qsort( mem1, mem1size, sizeof( int ), (int (*)())intcompare );
2522	qsort( mem2, mem2size, sizeof( int ), (int (*)())intcompare );
2523	// selhead( mem1, numin[0] );
2524	// selhead( mem2, numin[1] );
2525
2526
2527	#if 0
2528	fprintf( stderr, "\n" );
2529	fprintf( stderr, "mem1 (nin=%d) = \n", nin );
2530	for( i=0; ; i++ )
2531	{
2532	fprintf( stderr, "%d ", mem1[i]+1 );
2533	if( mem1[i] == -1 ) break;
2534	}
2535	fprintf( stderr, "\n" );
2536	fprintf( stderr, "mem2 (nin=%d) = \n", nin );
2537	for( i=0; ; i++ )
2538	{
2539	fprintf( stderr, "%d ", mem2[i]+1 );
2540	if( mem2[i] == -1 ) break;
2541	}
2542	fprintf( stderr, "\n" );
2543	#endif
2544
2545	#if 0
2546	fprintf( stderr, "before pairalign, l = %d, nyuko=%d, mem1size=%d, mem2size=%d\n", l, nyuko, mem1size, mem2size );
2547	fprintf( stderr, "before alignment\n" );
2548	for( j=0; j<mem1size; j++ )
2549	fprintf( stderr, "%s\n", seq[mem1[j]] );
2550	fprintf( stderr, "----\n" );
2551	for( j=0; j<mem2size; j++ )
2552	fprintf( stderr, "%s\n", seq[mem2[j]] );
2553	fprintf( stderr, "----\n\n" );
2554	#endif
2555
2556	if( val )
2557	{
2558	fprintf( stderr, "\r Alignment %d-%d \r", mem1size, mem2size );
2559	if( mem1 < mem2 )
2560	pairalign( njob, nlen, seq, mem1, mem2, weight, alloclen );
2561	else
2562	pairalign( njob, nlen, seq, mem2, mem1, weight, alloclen );
2563	}
2564
2565	#if TREE
2566	if( treeout )
2567	{
2568	v1 = topol[l][0][0];
2569	v2 = topol[l][1][0];
2570
2571	// fprintf( stderr, "nyuko=%d, v1=%d, v2=%d\n", nyuko, v1, v2 );
2572	if( v1 > v2 )
2573	{
2574	v3 = v1;
2575	v1 = v2;
2576	v2 = v3;
2577	}
2578	// fprintf( stderr, "nyuko=%d, v1=%d, v2=%d after sort\n", nyuko, v1, v2 );
2579	// fprintf( stderr, "nyuko=%d, v1=%d, v2=%d\n", nyuko, v1, v2 );
2580	// fprintf( stderr, "v1=%d, v2=%d, parttree[v1]=%s, parttree[v2]=%s\n", v1, v2, parttree[v1], parttree[v2] );
2581	sprintf( *tree, "(%s,%s)", parttree[v1], parttree[v2] );
2582	strcpy( parttree[v1], *tree );
2583	// fprintf( stderr, "parttree[%d] = %s\n", v1, parttree[v1] );
2584	// fprintf( stderr, "tree = %s\n", tree );
2585	free( parttree[v2] ); parttree[v2] = NULL;
2586	}
2587	#endif
2588
2589	#if 0
2590	fprintf( stderr, "after alignment\n" );
2591	for( j=0; j<mem1size; j++ )
2592	fprintf( stderr, "%s\n", seq[mem1[j]] );
2593	fprintf( stderr, "----\n" );
2594	for( j=0; j<mem2size; j++ )
2595	fprintf( stderr, "%s\n", seq[mem2[j]] );
2596	fprintf( stderr, "----\n\n" );
2597	#endif
2598	}
2599	#if TREE
2600	if( treeout )
2601	{
2602	free( parttree[v1] ); parttree[v1] = NULL;
2603	// fprintf( stderr, "tree = %s\n", tree );
2604	// FreeCharMtx( parttree );
2605	free( parttree ); parttree = NULL;
2606	}
2607	#endif
2608
2609	#if 0
2610	fprintf( stderr, "after alignment\n" );
2611	for( i=0; i<nyuko; i++ )
2612	{
2613	for( j=0; j<numin[i]; j++ )
2614	fprintf( stderr, "%s\n", seq[outs[i][j].numinseq] );
2615	}
2616	#endif
2617
2618	if( val )
2619	{
2620	groupid++;
2621	mptr = mem1; while( *mptr != -1 )
2622	{
2623	#if 0
2624	fprintf( stdout, "==g1-%d \n", *mptr+1 );
2625	fprintf( stdout, "%s \n", seq[*mptr] );
2626	#endif
2627	whichgroup[*mptr] = groupid;
2628	weight[mptr++] = 0.5;
2629	}
2630
2631	mptr = mem2; while( *mptr != -1 )
2632	{
2633	#if 0
2634	fprintf( stdout, "=g2-%d ", *mptr+1 );
2635	fprintf( stdout, "%s \n", seq[*mptr] );
2636	#endif
2637	whichgroup[*mptr] = groupid;
2638	weight[mptr++] = 0.5;
2639	}
2640
2641	if( numin[1] == 0 )
2642	{
2643	mptr = mem1; while( *mptr != -1 )
2644	{
2645	whichgroup[*mptr] = groupid;
2646	weight[mptr++] = (double)2.0/numin[0];
2647	}
2648	}
2649	}
2650	{
2651	if( depthpt > maxdepth ) maxdepth = depthpt;
2652	(*depthpt)++;
2653	}
2654	}
2655	else
2656	{
2657	#if TREE
2658	if( treeout )
2659	{
2660	sprintf( *tree, "%s", children[0] );
2661	free( children[0] );
2662	free( children );
2663	}
2664	#endif
2665	}
2666	for( i=0; i<npick; i++ ) free( (void *)outs[i] );
2667	// FreeFloatHalfMtx( pickmtx, npick );
2668	// FreeFloatHalfMtx( yukomtx, npick );
2669	FreeFloatMtx( len );
2670	FreeIntCub( topol );
2671	FreeIntVec( treeorder );
2672	free( outs );
2673	free( numin );
2674	free( picks );
2675	free( yukos );
2676	free( s_p_map );
2677	free( s_y_map );
2678	free( p_o_map );
2679	free( y_o_map );
2680	free( hanni );
2681	free( closeh );
2682	free( pickkouho );
2683	free( tsukau );
2684	// free( minscoreinpick );
2685	return val;
2686	}
2687	#endif
2688
2689	static void alignparaphiles( int nseq, int nlen, double weight, char *seq, int nmem, int members, int *alloclen )
2690	{
2691	int i, ilim;
2692	int *mem1 = AllocateIntVec( nmem );
2693	int *mem2 = AllocateIntVec( 2 );
2694
2695	mem2[1] = -1;
2696	ilim = nmem-1;
2697	for( i=0; i<ilim; i++ )
2698	{
2699	mem1[i] = members[i];
2700	mem1[i+1] = -1;
2701	mem2[0] = members[i+1];
2702	pairalign( nseq, nlen, seq, mem1, mem2, weight, alloclen );
2703	}
2704	free( mem1 );
2705	free( mem2 );
2706	}
2707
2708
2709
2710
2711
2712
2713
2714
2715
2716
2717
2718	int main( int argc, char *argv[] )
2719	{
2720	static char name, seq, **orialn;
2721	static int *grpseq;
2722	static char *tmpseq;
2723	static int **pointt;
2724	static int *nlen;
2725	int i, st, en;
2726	FILE *infp;
2727	FILE *treefp;
2728	char *treefile = NULL; //by Mathog
2729	char c;
2730	int alloclen;
2731	static int *order;
2732	static int *whichgroup;
2733	static double *weight;
2734	static char tmpname[B+100];
2735	int groupnum;
2736	int groupid;
2737	int pos;
2738	int *paramem;
2739	int npara;
2740	int completed;
2741	int orilen;
2742	int pscore;
2743	char *pt;
2744	int **tmpaminodis;
2745	static char com[1000];
2746	int depth;
2747	int aan;
2748
2749	static Scores *scores;
2750	static short *table1;
2751	static char **tree;
2752
2753
2754
2755	arguments( argc, argv );
2756
2757	if( inputfile )
2758	{
2759	infp = fopen( inputfile, "r" );
2760	if( !infp )
2761	{
2762	fprintf( stderr, "Cannot open %s\n", inputfile );
2763	exit( 1 );
2764	}
2765	}
2766	else
2767	infp = stdin;
2768
2769	getnumlen( infp );
2770	rewind( infp );
2771
2772	if( njob < 2 )
2773	{
2774	fprintf( stderr, "At least 2 sequences should be input!\n"
2775	"Only %d sequence found.\n", njob );
2776	exit( 1 );
2777	}
2778
2779
2780	if( picksize == NOTSPECIFIED \|\| picksize < 2 )
2781	picksize = PICKSIZE;
2782
2783	if( classsize == NOTSPECIFIED \|\| classsize < 0 )
2784	{
2785	classsize = njob + 1;
2786	}
2787	else
2788	{
2789	// picksize = MIN( picksize, (int)sqrt( classsize ) + 1);
2790	}
2791
2792	if( tokyoripara == NOTSPECIFIED )
2793	{
2794	if( doalign ) tokyoripara = TOKYORIPARA_A;
2795	else tokyoripara = TOKYORIPARA;
2796	}
2797	if( picksize > njob )
2798	tokyoripara = 0.0;
2799
2800
2801	alloclen = nlenmax * 2;
2802	name = AllocateCharMtx( njob, B+1 );
2803
2804	if( classsize == 1 )
2805	seq = AllocateCharMtx( njob, 0 );
2806	else
2807	seq = AllocateCharMtx( njob, alloclen+1 );
2808
2809
2810	nlen = AllocateIntVec( njob );
2811	tmpseq = calloc( nlenmax+1, sizeof( char ) );
2812	pointt = AllocateIntMtx( njob, 0 );
2813	grpseq = AllocateIntVec( nlenmax + 1 );
2814	order = (int *)calloc( njob + 1, sizeof( int ) );
2815	whichgroup = (int *)calloc( njob, sizeof( int ) );
2816	weight = (double *)calloc( njob, sizeof( double ) );
2817
2818	fprintf( stderr, "alloclen = %d in main\n", alloclen );
2819
2820	for( i=0; i<njob; i++ ) whichgroup[i] = 0;
2821	for( i=0; i<njob; i++ ) weight[i] = 1.0;
2822	for( i=0; i<njob; i++ ) order[i] = -1;
2823
2824	if( classsize == 1 )
2825	readData_varlen( infp, name, nlen, seq );
2826	else
2827	readData_pointer( infp, name, nlen, seq );
2828
2829	fclose( infp );
2830
2831	if( fromaln ) doalign = 1;
2832
2833	if( fromaln )
2834	{
2835	orialn = AllocateCharMtx( njob, alloclen+1 );
2836	for( i=0; i<njob; i++ )
2837	{
2838	if( strlen( seq[i] ) != nlenmax )
2839	{
2840	fprintf( stderr, "Input sequences must be aligned\n" );
2841	exit( 1 );
2842	}
2843	strcpy( orialn[i], seq[i] );
2844	}
2845	}
2846
2847	constants( njob, seq );
2848
2849	if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
2850	else tsize = (int)pow( 6, 6 );
2851
2852	if( dorp == 'd' )
2853	{
2854	lenfaca = DLENFACA;
2855	lenfacb = DLENFACB;
2856	lenfacc = DLENFACC;
2857	lenfacd = DLENFACD;
2858	}
2859	else
2860	{
2861	lenfaca = PLENFACA;
2862	lenfacb = PLENFACB;
2863	lenfacc = PLENFACC;
2864	lenfacd = PLENFACD;
2865	}
2866
2867	maxl = 0;
2868	for( i=0; i<njob; i++ )
2869	{
2870	gappick0( tmpseq, seq[i] );
2871	nlen[i] = strlen( tmpseq );
2872	strcpy( seq[i], tmpseq );
2873	pointt[i] = AllocateIntVec( nlen[i]+1 ); // ??
2874
2875	if( nlen[i] < 6 )
2876	{
2877	fprintf( stderr, "Seq %d, too short, %d characters\n", i+1, nlen[i] );
2878	fprintf( stderr, "name = %s\n", name[i] );
2879	fprintf( stderr, "seq = %s\n", seq[i] );
2880	exit( 1 );
2881	// continue;
2882	}
2883	if( nlen[i] > maxl ) maxl = nlen[i];
2884	if( dorp == 'd' ) /* nuc */
2885	{
2886	if( seq_grp_nuc( grpseq, tmpseq ) < 6 )
2887	{
2888	fprintf( stderr, "Seq %d, too short.\n", i+1 );
2889	fprintf( stderr, "name = %s\n", name[i] );
2890	fprintf( stderr, "seq = %s\n", seq[i] );
2891	exit( 1 );
2892	// continue;
2893	}
2894	makepointtable_nuc( pointt[i], grpseq );
2895	}
2896	else /* amino */
2897	{
2898	if( seq_grp( grpseq, tmpseq ) < 6 )
2899	{
2900	fprintf( stderr, "Seq %d, too short.\n", i+1 );
2901	fprintf( stderr, "name = %s\n", name[i] );
2902	fprintf( stderr, "seq = %s\n", seq[i] );
2903	exit( 1 );
2904	// continue;
2905	}
2906	makepointtable( pointt[i], grpseq );
2907	}
2908	// fprintf( stdout, ">%s\n", name[i] );
2909	// fprintf( stdout, "%s\n", seq[i] );
2910	}
2911	if( nunknown ) fprintf( stderr, "\nWARNING : %d unknown characters\n", nunknown );
2912	// exit( 1 );
2913
2914	#if 0
2915	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
2916	#endif
2917
2918	initSignalSM();
2919
2920	initFiles();
2921
2922	WriteOptions( trap_g );
2923
2924	c = seqcheck( seq );
2925	if( c )
2926	{
2927	fprintf( stderr, "Illeagal character %c\n", c );
2928	exit( 1 );
2929	}
2930
2931	pid = (int)getpid();
2932	sprintf( datafile, "/tmp/data-%d", pid );
2933	sprintf( queryfile, "/tmp/query-%d", pid );
2934	sprintf( resultfile, "/tmp/fasta-%d", pid );
2935
2936	scores = (Scores *)calloc( njob, sizeof( Scores ) );
2937
2938	// fprintf( stderr, "\nCalculating i-i scores ... \n" );
2939	for( i=0; i<njob; i++ )
2940	{
2941	orilen = strlen( seq[i] );
2942	scores[i].numinseq = i; // irimasu
2943	scores[i].orilen = orilen;
2944	}
2945
2946	if( doalign == 'f' )
2947	{
2948	fastapath = getenv( "FASTA_4_MAFFT" );
2949	if( !fastapath ) fastapath = "fasta34";
2950	fprintf( stderr, "fastapath=%s\n", fastapath );
2951	tmpaminodis = AllocateIntMtx( 0x80, 0x80 );
2952	getfastascoremtx( tmpaminodis );
2953	}
2954	else
2955	tmpaminodis = NULL;
2956
2957	for( i=0; i<njob; i++ )
2958	{
2959	scores[i].pointt = pointt[i];
2960	scores[i].shimon = (int)pointt[i][0] + (int)pointt[i][1] + (int)pointt[i][scores[i].orilen-6];
2961	scores[i].name = name[i];
2962	if( doalign )
2963	{
2964	fprintf( stderr, "\r %05d/%05d ", i, njob );
2965	free( scores[i].pointt );
2966	if( doalign == 'f' )
2967	{
2968	#if 0
2969	#define KIZAMI 100
2970	int ipos = (int)( i / KIZAMI ) * KIZAMI;
2971	int iposamari = i % KIZAMI;
2972
2973	fprintf( stderr, "%d / %d\r", i, njob );
2974	// fprintf( stderr, "ipos = %d\n", ipos );
2975	// fprintf( stderr, "iposamari = %d\n", iposamari );
2976
2977	// fprintf( stderr, " calling blast, i=%d\n", i );
2978	// blastresults = callfasta( seq, scores+i, 1, 0, 1 );
2979	blastresults = callfasta( seq, scores+ipos, MIN(KIZAMI,njob-ipos), NULL, iposamari, (iposamari==0) );
2980	// fprintf( stderr, "done., i=%d\n\n", i );
2981	scores[i].selfscore = (int)blastresults[iposamari];
2982	#if 0
2983	for( j=0; j<100; j++ )
2984	{
2985	fprintf( stderr, "res[%d] = %f\n", j, blastresults[j] );
2986	}
2987	#endif
2988	// fprintf( stderr, "%d->selfscore = %d\n", i, scores[i].selfscore );
2989	free( blastresults );
2990	#else
2991	pscore = 0;
2992	if( scoremtx == -1 )
2993	{
2994	st = 1;
2995	en = 0;
2996	for( pt=seq[i]; *pt; pt++ )
2997	{
2998	if( pt == 'u' ) pt = 't';
2999	aan = amino_n[(int)*pt];
3000	if( aan<0 \|\| aan >= 4 ) *pt = 'n';
3001
3002	if( *pt == 'n' )
3003	{
3004	en++;
3005	if( st ) continue;
3006	else pscore += tmpaminodis[(int)pt][(int)pt];
3007	}
3008	else
3009	{
3010	st = 0;
3011	en = 0;
3012	pscore += tmpaminodis[(int)pt][(int)pt];
3013	}
3014	}
3015	scores[i].selfscore = pscore - en * tmpaminodis['n']['n'];
3016	}
3017	else
3018	{
3019	st = 1;
3020	en = 0;
3021	for( pt=seq[i]; *pt; pt++ )
3022	{
3023	aan = amino_n[(int)*pt];
3024	if( aan<0 \|\| aan >= 20 ) *pt = 'X';
3025	if( *pt == 'X' )
3026	{
3027	en++;
3028	if( st ) continue;
3029	else pscore += tmpaminodis[(int)pt][(int)pt];
3030	}
3031	else
3032	{
3033	st = 0;
3034	en = 0;
3035	pscore += tmpaminodis[(int)pt][(int)pt];
3036	}
3037	}
3038	scores[i].selfscore = pscore - en * tmpaminodis['X']['X'];
3039	}
3040	#endif
3041	}
3042	else
3043	{
3044	pscore = 0;
3045	for( pt=seq[i]; *pt; pt++ )
3046	{
3047	pscore += amino_dis[(int)pt][(int)pt];
3048	}
3049	scores[i].selfscore = pscore;
3050	}
3051	// fprintf( stderr, "selfscore[%d] = %d\n", i+1, scores[i].selfscore );
3052	}
3053	else
3054	{
3055	table1 = (short *)calloc( tsize, sizeof( short ) );
3056	if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
3057	makecompositiontable_p( table1, pointt[i] );
3058	scores[i].selfscore = commonsextet_p( table1, pointt[i] );
3059	free( table1 );
3060	}
3061	}
3062	if( tmpaminodis ) FreeIntMtx( tmpaminodis );
3063
3064	depth = 0;
3065	#if TREE
3066	if( treeout )
3067	{
3068	tree = (char *)calloc( 1, sizeof( char ) );
3069	*tree = NULL;
3070	// splitseq_bin( scores, njob, nlen, seq, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight );
3071	completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
3072	treefile = (char *)calloc( strlen( inputfile ) + 10, sizeof( char ) );
3073	if( inputfile )
3074	sprintf( treefile, "%s.tree", inputfile );
3075	else
3076	sprintf( treefile, "splittbfast.tree" );
3077	treefp = fopen( treefile, "w" );
3078	fprintf( treefp, "%s\n", *tree );
3079	fclose( treefp );
3080	}
3081	else
3082	completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
3083	#else
3084	completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
3085	#endif
3086
3087	fprintf( stderr, "\nDone.\n\n" );
3088
3089	#if 1
3090	groupnum = 0;
3091	groupid = -1;
3092	paramem = NULL;
3093	npara = 0;
3094	for( i=0; i<njob; i++ )
3095	{
3096	pos = order[i];
3097	if( whichgroup[pos] != groupid )
3098	{
3099	groupnum++;
3100	groupid = whichgroup[pos];
3101	}
3102	if( whichgroup[pos] )
3103	{
3104	if( paramem )
3105	{
3106	paramem[npara] = -1;
3107	if( npara > 1 && classsize > 2 )
3108	{
3109	qsort( paramem, npara, sizeof( int ), (int ()(const void , const void*))intcompare );
3110	// selhead( paramem, npara );
3111	alignparaphiles( njob, nlen, weight, seq, npara, paramem, &alloclen );
3112	}
3113	free( paramem ); paramem = NULL; npara = 0;
3114	}
3115	sprintf( tmpname, "Group-%d %s", groupnum, name[pos]+1 );
3116	}
3117	else
3118	{
3119	paramem = realloc( paramem, sizeof( int) * ( npara + 2 ) );
3120	paramem[npara++] = pos;
3121	sprintf( tmpname, "Group-para %s", name[pos]+1 );
3122	}
3123	tmpname[B-1] = 0;
3124	if( classsize > 1 && classsize <= njob )
3125	strcpy( name[pos]+1, tmpname );
3126	}
3127	if( paramem )
3128	{
3129	paramem[npara] = -1;
3130	if( npara > 1 && classsize > 2 )
3131	{
3132	qsort( paramem, npara, sizeof( int ), (int ()(const void , const void*))intcompare );
3133	// selhead( paramem, npara );
3134	alignparaphiles( njob, nlen, weight, seq, npara, paramem, &alloclen );
3135	}
3136	free( paramem ); paramem = NULL; npara = 0;
3137	}
3138	#else
3139	for( i=0; i<njob; i++ )
3140	{
3141	sprintf( tmpname, "Group-%d %s", whichgroup[i], name[i]+1 );
3142	strcpy( name[i]+1, tmpname );
3143	}
3144	#endif
3145
3146
3147	// maketanni( name, seq, njob, nlenmax, nlen );
3148
3149	fclose( trap_g );
3150
3151	#if DEBUG
3152	fprintf( stderr, "writing alignment to stdout\n" );
3153	#endif
3154	if( reorder )
3155	writeData_reorder_pointer( stdout, njob, name, nlen, seq, order );
3156	else
3157	writeData_pointer( stdout, njob, name, nlen, seq );
3158	#if IODEBUG
3159	fprintf( stderr, "OSHIMAI\n" );
3160	#endif
3161	if( classsize == 1 )
3162	{
3163	fprintf( stderr, "\n\n" );
3164	fprintf( stderr, "----------------------------------------------------------------------------\n" );
3165	fprintf( stderr, "\n" );
3166	fprintf( stderr, "nseq = %d\n", njob );
3167	fprintf( stderr, "groupsize = %d, picksize=%d\n", classsize, picksize );
3168	fprintf( stderr, "The input sequences have been sorted so that similar sequences are close.\n" );
3169	if( reorder )
3170	fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
3171	#if TREE
3172	if( treeout )
3173	{
3174	fprintf( stderr, "\n" );
3175	fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
3176	}
3177	// else
3178	// {
3179	// fprintf( stderr, "To output guide tree,\n" );
3180	// fprintf( stderr, "%% %s -t -i %s\n", progName( argv[0] ), "inputfile" );
3181	// }
3182	#endif
3183	if( !doalign )
3184	{
3185	fprintf( stderr, "\n" );
3186	fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
3187	fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
3188	}
3189	fprintf( stderr, "\n" );
3190	fprintf( stderr, "----------------------------------------------------------------------------\n" );
3191	}
3192	else if( groupnum > 1 )
3193	{
3194	fprintf( stderr, "\n\n" );
3195	fprintf( stderr, "----------------------------------------------------------------------------\n" );
3196	fprintf( stderr, "\n" );
3197	fprintf( stderr, "groupsize = %d, picksize=%d\n", classsize, picksize );
3198	fprintf( stderr, "The input sequences have been classified into %d groups + some paraphyletic groups\n", groupnum );
3199	fprintf( stderr, "Note that the alignment is not completed.\n" );
3200	if( reorder )
3201	fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
3202	#if TREE
3203	if( treeout )
3204	{
3205	fprintf( stderr, "\n" );
3206	fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
3207	}
3208	// else
3209	// {
3210	// fprintf( stderr, "To output guide tree,\n" );
3211	// fprintf( stderr, "%% %s -t -i %s\n", progName( argv[0] ), "inputfile" );
3212	// }
3213	#endif
3214	if( !doalign )
3215	{
3216	fprintf( stderr, "\n" );
3217	fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
3218	fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
3219	}
3220	fprintf( stderr, "\n" );
3221	fprintf( stderr, "----------------------------------------------------------------------------\n" );
3222	}
3223	else
3224	{
3225	fprintf( stderr, "\n\n" );
3226	fprintf( stderr, "----------------------------------------------------------------------------\n" );
3227	fprintf( stderr, "\n" );
3228	fprintf( stderr, "nseq = %d\n", njob );
3229	fprintf( stderr, "groupsize = %d, partsize=%d\n", classsize, picksize );
3230	// fprintf( stderr, "A single alignment containing all the input sequences has been computed.\n" );
3231	// fprintf( stderr, "If the sequences are highly diverged and you feel there are too many gaps,\n" );
3232	// fprintf( stderr, "please try \n" );
3233	// fprintf( stderr, "%% mafft --parttree --groupsize 100 inputfile\n" );
3234	// fprintf( stderr, "which classifies the sequences into several groups with <~ 100 sequences\n" );
3235	// fprintf( stderr, "and performs only intra-group alignments.\n" );
3236	if( reorder )
3237	fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
3238	#if TREE
3239	if( treeout )
3240	{
3241	fprintf( stderr, "\n" );
3242	fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
3243	}
3244	// else
3245	// {
3246	// fprintf( stderr, "To output guide tree,\n" );
3247	// fprintf( stderr, "%% %s -t -i %s\n", progName( argv[0] ), "inputfile" );
3248	// }
3249	#endif
3250	if( !doalign \|\| fromaln )
3251	{
3252	fprintf( stderr, "\n" );
3253	fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
3254	fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
3255	}
3256	fprintf( stderr, "\n" );
3257	fprintf( stderr, "----------------------------------------------------------------------------\n" );
3258	}
3259	#if TREE
3260	if( treeout ) free( treefile );
3261	#endif
3262
3263	#if 0
3264	fprintf( stdout, "weight =\n" );
3265	for( i=0; i<njob; i++ )
3266	fprintf( stdout, "%d: %f\n", i+1, weight[i] );
3267	#endif
3268
3269	if( doalign == 'f' )
3270	{
3271	strcpy( com, "rm -f" );
3272	strcat( com, " " );
3273	strcat( com, datafile );
3274	strcat( com, "* " );
3275	strcat( com, queryfile );
3276	strcat( com, " " );
3277	strcat( com, resultfile );
3278	fprintf( stderr, "%s\n", com );
3279	system( com );
3280	}
3281
3282	SHOWVERSION;
3283
3284	return( 0 );
3285	}

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

Context Navigation

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

Download in other formats: