Context Navigation

addsingle.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: 76.3 KB

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

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

Context Navigation

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

Download in other formats: