text: #include "mltaln.h"
  
  
  #define TREE 1
  #define PICKSIZE 50 // must be >= 3
  #define WEIGHT 0
  #define TOKYORIPARA 0.70 // 0.70
  #define TOKYORIPARA_A 0.70  // changed
  #define LENFAC 1
  #define HUKINTOTREE 1
  #define DIANA 0
  #define MAX6DIST 10.0
  
  // kouzoutai ni sasareru pointer ha static
  
  #define DEBUG 0
  #define IODEBUG 0
  #define SCOREOUT 0
  
  #define END_OF_VEC -1
  
  static char *fastapath;
  static int doalign;
  static int fromaln;
  static int uselongest;
  static int treeout;
  static int classsize;
  static int picksize;
  static int maxl;
  static int tsize;
  static int reorder;
  static int pid;
  static int maxdepth = 0;
  static double tokyoripara;
  
  static double lenfaca, lenfacb, lenfacc, lenfacd;
  #define PLENFACA 0.01
  #define PLENFACB 10000
  #define PLENFACC 10000
  #define PLENFACD 0.1
  #define DLENFACA 0.01
  #define DLENFACB 2500
  #define DLENFACC 2500
  #define DLENFACD 0.1
  
  static char datafile[1000];
  static char queryfile[1000];
  static char resultfile[1000];
  
  typedef struct _scores
  {
  	double score;
  	int selfscore;
  	int orilen;
  	int *pointt;
  	int numinseq;
  	char *name;
  //	char *seq; // reallo
  //	char **seqpt;
  	int shimon;
  } Scores;
  
  int intcompare( const int *a, const int *b )
  {
  	return( *a - *b );
  }
  
  int lcompare( const Scores *a, const Scores *b )
  {
  	if( a->orilen < b->orilen ) return 1;
  	else if( a->orilen > b->orilen ) return -1;
  	else return 0;
  }
  
  int dcompare( const Scores *a, const Scores *b )
  {
  	if( a->score > b->score ) return 1;
  	else if( a->score < b->score ) return -1;
  	else
  	{
  		if( a->selfscore < b->selfscore ) return 1;
  		else if( a->selfscore > b->selfscore ) return -1;
  		else 
  		{
  			if( a->orilen < b->orilen ) return 1;
  			else if( a->orilen > b->orilen ) return -1;
  			else return 0;
  		}
  	}
  }
  
  #if 0
  static void	gappickandx0( char *out, char *in )
  {
  	char c;
  	if( scoremtx == -1 )
  	{
  		while( *in )
  		{
  			if( (c=*in++) == '-' )
  				;	
  			else if( c == 'u' )
  				*out++ = 't';
  			else if( amino_n[c] < 4 && amino_n[c] > -1 )
  				*out++ = c;
  			else
  				*out++ = 'n';
  		}
  	}
  	else
  	{
  		while( *in )
  		{
  			if( (c=*in++) == '-' )
  				;	
  			else if( amino_n[c] < 20 && amino_n[c] > -1 )
  				*out++ = c;
  			else
  				*out++ = 'X';
  		}
  	}
  	*out = 0;
  }	
  
  static int getkouho( int *pickkouho, double prob, int nin, Scores *scores, char **seq ) // 0 < prob < 1
  {
  	int nkouho = 0;
  	int i, j;
  	int *iptr = pickkouho;
  	for( i=1; i<nin; i++ )
  	{
  		if( ( nkouho==0 || rnd() < prob ) && ( scores[i].shimon != scores->shimon || strcmp( seq[scores->numinseq], seq[scores[i].numinseq] ) ) )
  		{
  #if 0
  			for( j=0; j<nkouho; j++ )
  			{
  				if( scores[i].shimon == scores[pickkouho[j]].shimon || !strcmp( seq[scores[pickkouho[j]].numinseq], seq[scores[i].numinseq] ) ) 
  					break;
  			}
  			if( j == nkouho )
  #endif
  			{
  				*iptr++ = i;
  				nkouho++;
  //				fprintf( stderr, "ok! nkouho=%d\n", nkouho );
  			}
  		}
  		else
  		{
  			;
  //			fprintf( stderr, "no! %d-%d\n", 0, scores[i].numinseq );
  		}
  	}
  	fprintf( stderr, "\ndone\n\n"  );
  	return nkouho;
  }
  
  #endif
  
  static void getfastascoremtx( int **tmpaminodis )
  {
  	FILE *qfp;
  	FILE *dfp;
  	FILE *rfp;
  	int i, j;
  	char aa;
  	int slen;
  	int res;
  	char com[10000];
  	static char *tmpseq;
  	static char *tmpname;
  	double *resvec;
  
  	if( scoremtx == -1 )
  	{
  		tmpaminodis['a']['a'] = 5;
  		tmpaminodis['g']['g'] = 5;
  		tmpaminodis['c']['c'] = 5;
  		tmpaminodis['t']['t'] = 5;
  		tmpaminodis['n']['n'] = -1;
  
  		return;
  	}
  
  
  	tmpseq = calloc( 2000, sizeof( char ) );
  	tmpname = calloc( B, sizeof( char ) );
  	resvec = calloc( 1, sizeof( double ) );
  
  //	fprintf( stderr, "xformatting .. " );
  	dfp = fopen( datafile, "w" );
  	if( !dfp ) ErrorExit( "Cannot open datafile." );
  	sprintf( tmpname, ">+===========+%d                      ", 0 );
  	strcpy( tmpseq, "AAAAAAXXXXXX" );
  	strcat( tmpseq, "CCCCCCXXXXXX" );
  	strcat( tmpseq, "DDDDDDXXXXXX" );
  	strcat( tmpseq, "EEEEEEXXXXXX" );
  	strcat( tmpseq, "FFFFFFXXXXXX" );
  	strcat( tmpseq, "GGGGGGXXXXXX" );
  	strcat( tmpseq, "HHHHHHXXXXXX" );
  	strcat( tmpseq, "IIIIIIXXXXXX" );
  	strcat( tmpseq, "KKKKKKXXXXXX" );
  	strcat( tmpseq, "LLLLLLXXXXXX" );
  	strcat( tmpseq, "MMMMMMXXXXXX" );
  	strcat( tmpseq, "NNNNNNXXXXXX" );
  	strcat( tmpseq, "PPPPPPXXXXXX" );
  	strcat( tmpseq, "QQQQQQXXXXXX" );
  	strcat( tmpseq, "RRRRRRXXXXXX" );
  	strcat( tmpseq, "SSSSSSXXXXXX" );
  	strcat( tmpseq, "TTTTTTXXXXXX" );
  	strcat( tmpseq, "VVVVVVXXXXXX" );
  	strcat( tmpseq, "WWWWWWXXXXXX" );
  	strcat( tmpseq, "YYYYYYXXXXXX" );
  	slen = strlen( tmpseq );
  	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
  	fclose( dfp );
  	fprintf( stderr, "done.\n" );
  
  	for( i=0; i<20; i++ )
  	{
  		aa = amino[i];
  //		fprintf( stderr, "checking %c\n", aa );
  		*tmpseq = 0;
  		sprintf( tmpname, ">+===========+%d                      ", 0 );
  		for( j=0; j<6; j++ )
  			sprintf( tmpseq+strlen( tmpseq ), "%c", aa );
  		qfp = fopen( queryfile, "w" );
  		if( !qfp ) ErrorExit( "Cannot open queryfile." );
  		writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
  		fclose( qfp );
  
  		if( scoremtx == -1 ) 
  			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 );
  		else
  			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 );
  		res = system( com );
  		if( res )
  		{
  			fprintf( stderr, "error in %s", fastapath );
  			exit( 1 );
  		}
  
  		rfp = fopen( resultfile, "r" );
  		if( rfp == NULL )  
  			ErrorExit( "file 'fasta.$$' does not exist\n" );
  		res = ReadFasta34m10_scoreonly( rfp, resvec, 1 );
  		fprintf( stderr, "%c: %f\n", 'A'+i, *resvec/6 );
  		fclose( rfp );
  		if( ( (int)*resvec % 6 ) > 0.0 )
  		{
  			fprintf( stderr, "Error in blast, *resvec=%f\n", *resvec );
  			fprintf( stderr, "Error in blast, *resvec/6=%f\n", *resvec/6 );
  			exit( 1 );
  		}
  		tmpaminodis[(int)aa][(int)aa] = (int)( *resvec / 6 );
  //		fprintf( stderr, "*resvec=%f, tmpaminodis[aa][aa] = %d\n", *resvec, tmpaminodis[aa][aa] );
  	}
  	tmpaminodis['X']['X'] = -1;
  	free( tmpname );
  	free( tmpseq );
  	free( resvec );
  }
  
  #if 0
  static void getblastscoremtx( int **tmpaminodis )
  {
  	FILE *qfp;
  	FILE *dfp;
  	FILE *rfp;
  	int i, j;
  	char aa;
  	int slen;
  	int res;
  	char com[10000];
  	static char *tmpseq;
  	static char *tmpname;
  	double *resvec;
  
  	if( scoremtx == -1 )
  	{
  		tmpaminodis['a']['a'] = 1;
  		tmpaminodis['g']['g'] = 1;
  		tmpaminodis['c']['c'] = 1;
  		tmpaminodis['t']['t'] = 1;
  
  		return;
  	}
  
  
  	tmpseq = calloc( 2000, sizeof( char ) );
  	tmpname = calloc( B, sizeof( char ) );
  	resvec = calloc( 1, sizeof( double ) );
  
  //	fprintf( stderr, "xformatting .. " );
  	dfp = fopen( datafile, "w" );
  	if( !dfp ) ErrorExit( "Cannot open datafile." );
  	sprintf( tmpname, "\0", i ); // BUG!!
  	strcpy( tmpseq, "AAAAAAXXXXXX" );
  	strcat( tmpseq, "CCCCCCXXXXXX" );
  	strcat( tmpseq, "DDDDDDXXXXXX" );
  	strcat( tmpseq, "EEEEEEXXXXXX" );
  	strcat( tmpseq, "FFFFFFXXXXXX" );
  	strcat( tmpseq, "GGGGGGXXXXXX" );
  	strcat( tmpseq, "HHHHHHXXXXXX" );
  	strcat( tmpseq, "IIIIIIXXXXXX" );
  	strcat( tmpseq, "KKKKKKXXXXXX" );
  	strcat( tmpseq, "LLLLLLXXXXXX" );
  	strcat( tmpseq, "MMMMMMXXXXXX" );
  	strcat( tmpseq, "NNNNNNXXXXXX" );
  	strcat( tmpseq, "PPPPPPXXXXXX" );
  	strcat( tmpseq, "QQQQQQXXXXXX" );
  	strcat( tmpseq, "RRRRRRXXXXXX" );
  	strcat( tmpseq, "SSSSSSXXXXXX" );
  	strcat( tmpseq, "TTTTTTXXXXXX" );
  	strcat( tmpseq, "VVVVVVXXXXXX" );
  	strcat( tmpseq, "WWWWWWXXXXXX" );
  	strcat( tmpseq, "YYYYYYXXXXXX" );
  	slen = strlen( tmpseq );
  	writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
  	fclose( dfp );
  	if( scoremtx == -1 )
  		sprintf( com, "formatdb  -p f -i %s -o F", datafile );
  	else
  		sprintf( com, "formatdb  -i %s -o F", datafile );
  	system( com );
  	fprintf( stderr, "done.\n" );
  
  	for( i=0; i<20; i++ )
  	{
  		aa = amino[i];
  		fprintf( stderr, "checking %c\n", aa );
  		*tmpseq = 0;
  		for( j=0; j<6; j++ )
  			sprintf( tmpseq+strlen( tmpseq ), "%c", aa );
  		qfp = fopen( queryfile, "w" );
  		if( !qfp ) ErrorExit( "Cannot open queryfile." );
  		writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
  		fclose( qfp );
  
  		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 );
  		res = system( com );
  		if( res )
  		{
  			fprintf( stderr, "error in %s", "blastall" );
  			exit( 1 );
  		}
  
  		rfp = fopen( resultfile, "r" );
  		if( rfp == NULL )  
  			ErrorExit( "file 'fasta.$$' does not exist\n" );
  		res = ReadBlastm7_scoreonly( rfp, resvec, 1 );
  		fprintf( stdout, "%c: %f\n", 'A'+i, *resvec/6 );
  		fclose( rfp );
  		if( ( (int)*resvec % 6 ) > 0.0 )
  		{
  			fprintf( stderr, "Error in blast, *resvec=%f\n", *resvec );
  			fprintf( stderr, "Error in blast, *resvec/6=%f\n", *resvec/6 );
  			exit( 1 );
  		}
  		tmpaminodis[aa][aa] = (int)( *resvec / 6 );
  	}
  	tmpaminodis['X']['X'] = 0;
  	free( tmpname );
  	free( tmpseq );
  	free( resvec );
  	
  }
  #endif
  
  static double *callfasta( char **seq, Scores *scores, int nin, int *picks, int query, int rewritedata )
  {
  	double *val;
  	FILE *qfp;
  	FILE *dfp;
  	FILE *rfp;
  	int i;
  	char com[10000];
  	static char datafile[1000];
  	static char queryfile[1000];
  	static char resultfile[1000];
  	static int pid;
  	static char *tmpseq;
  	static char *tmpname;
  	int slen;
  	int res;
  	static Scores *scoresbk = NULL;
  	static int ninbk = 0;
  
  	if( pid == 0 )
  	{
  		pid = (int)getpid();
  		sprintf( datafile, "/tmp/data-%d", pid );
  		sprintf( queryfile, "/tmp/query-%d", pid );
  		sprintf( resultfile, "/tmp/fasta-%d", pid );
  
  		tmpseq = calloc( nlenmax+1, sizeof( char ) );
  		tmpname = calloc( B+1, sizeof( char ) );
  	}
  
  	val = calloc( nin, sizeof( double ) );
  //	fprintf( stderr, "nin=%d, q=%d\n", nin, query );
  
  	if( rewritedata )
  	{
  		scoresbk = scores;
  		ninbk = nin;
  //		fprintf( stderr, "\nformatting .. " );
  		dfp = fopen( datafile, "w" );
  		if( !dfp ) ErrorExit( "Cannot open datafile." );
  		if( picks == NULL ) for( i=0; i<nin; i++ )
  		{
  //			fprintf( stderr, "i=%d / %d / %d\n", i,  nin, njob );
  //			fprintf( stderr, "nlenmax = %d\n", nlenmax );
  //			fprintf( stderr, "scores[i].orilen = %d\n", scores[i].orilen );
  //			fprintf( stderr, "strlen( seq[scores[i].numinseq] = %d\n", strlen( seq[scores[i].numinseq] ) );
  			gappick0( tmpseq, seq[scores[i].numinseq] );
  			sprintf( tmpname, ">+===========+%d                      ", i );
  			slen = scores[i].orilen;
  			writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
  		}
  		else for( i=0; i<nin; i++ )
  		{
  			gappick0( tmpseq, seq[scores[picks[i]].numinseq] );
  			sprintf( tmpname, ">+===========+%d                      ", i );
  			slen = scores[picks[i]].orilen;
  			writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
  		}
  		fclose( dfp );
  	}
  
  
  	gappick0( tmpseq, seq[scores[query].numinseq] );
  	sprintf( tmpname, ">+==========+%d                      ", 0 );
  	slen = scores[query].orilen;
  	qfp = fopen( queryfile, "w" );
  	if( !qfp ) ErrorExit( "Cannot open queryfile." );
  	writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
  	fclose( qfp );
  
  //	fprintf( stderr, "calling fasta, nin=%d\n", nin );
  
  	if( scoremtx == -1 ) 
  		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 );
  	else
  		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 );
  	res = system( com );
  	if( res )
  	{
  		fprintf( stderr, "error in %s", fastapath );
  		exit( 1 );
  	}
  //	fprintf( stderr, "fasta done\n" );
  
  //exit( 1 );
  
  	rfp = fopen( resultfile, "r" );
  	if( rfp == NULL )  
  		ErrorExit( "file 'fasta.$$' does not exist\n" );
  
  //	fprintf( stderr, "reading fasta\n" );
  	if( scoremtx == -1 ) 
  		res = ReadFasta34m10_scoreonly_nuc( rfp, val, nin );
  	else
  		res = ReadFasta34m10_scoreonly( rfp, val, nin );
  //	fprintf( stderr, "done. val[0] = %f\n", val[0] );
  
  
  	fclose( rfp );
  
  #if 0
  	for( i=0; i<nin; i++ )
  		fprintf( stderr, "r[%d-%d] = %f\n", 0, i, val[i] );
  	exit( 1 );
  #endif
  
  	return( val );
  }
  #if 0
  static double *callblast( char **seq, Scores *scores, int nin, int query, int rewritedata )
  {
  	double *val;
  	FILE *qfp;
  	FILE *dfp;
  	FILE *rfp;
  	int i, j;
  	char com[10000];
  	static char datafile[1000];
  	static char queryfile[1000];
  	static char resultfile[1000];
  	static int pid;
  	static char *tmpseq;
  	static char *tmpname;
  	char *seqptr;
  	int slen;
  	int res;
  	static Scores *scoresbk = NULL;
  	static int ninbk = 0;
  
  	if( pid == 0 )
  	{
  		pid = (int)getpid();
  		sprintf( datafile, "/tmp/data-%d\0", pid );
  		sprintf( queryfile, "/tmp/query-%d\0", pid );
  		sprintf( resultfile, "/tmp/fasta-%d\0", pid );
  
  		tmpseq = calloc( nlenmax+1, sizeof( char ) );
  		tmpname = calloc( B+1, sizeof( char ) );
  	}
  
  	val = calloc( nin, sizeof( double ) );
  //	fprintf( stderr, "nin=%d, q=%d\n", nin, query );
  
  	if( rewritedata )
  	{
  		scoresbk = scores;
  		ninbk = nin;
  		fprintf( stderr, "\nformatting .. " );
  		dfp = fopen( datafile, "w" );
  		if( !dfp ) ErrorExit( "Cannot open datafile." );
  		for( i=0; i<nin; i++ )
  		{
  //			fprintf( stderr, "i=%d / %d / %d\n", i,  nin, njob );
  //			fprintf( stderr, "nlenmax = %d\n", nlenmax );
  //			fprintf( stderr, "scores[i].orilen = %d\n", scores[i].orilen );
  //			fprintf( stderr, "strlen( seq[scores[i].numinseq] = %d\n", strlen( seq[scores[i].numinseq] ) );
  			gappick0( tmpseq, seq[scores[i].numinseq] );
  			sprintf( tmpname, "+===========+%d                      \0", i );
  			slen = scores[i].orilen;
  			writeData_pointer( dfp, 1, &tmpname, &slen, &tmpseq );
  		}
  		fclose( dfp );
  			
  		if( scoremtx == -1 )
  			sprintf( com, "formatdb  -p f -i %s -o F", datafile );
  		else
  			sprintf( com, "formatdb  -i %s -o F", datafile );
  		system( com );
  //		fprintf( stderr, "done.\n" );
  	}
  
  
  	gappick0( tmpseq, seq[scores[query].numinseq] );
  	sprintf( tmpname, "+==========+%d                      \0", 0 );
  	slen = scores[query].orilen;
  	qfp = fopen( queryfile, "w" );
  	if( !qfp ) ErrorExit( "Cannot open queryfile." );
  	writeData_pointer( qfp, 1, &tmpname, &slen, &tmpseq );
  	fclose( qfp );
  //	fprintf( stderr, "q=%s\n", tmpseq );
  
  	fprintf( stderr, "\ncalling blast .. \n" );
  	if( scoremtx == -1 ) 
  		sprintf( com, "blastall -b %d -e 1e10 -p blastn -m 7  -i %s -d %s >  %s\0", nin, queryfile, datafile, resultfile );
  	else
  		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 );
  	res = system( com );
  	if( res ) ErrorExit( "error in blast" );
  
  	rfp = fopen( resultfile, "r" );
  	if( rfp == NULL )  
  		ErrorExit( "file 'fasta.$$' does not exist\n" );
  	res = ReadBlastm7_scoreonly( rfp, val, nin );
  	fclose( rfp );
  
  #if 0
  	for( i=0; i<nin; i++ )
  		fprintf( stderr, "r[%d-%d] = %f\n", 0, i, val[i] );
  #endif
  
  	return( val );
  }
  #endif
  
  #if 0
  static void selhead( int *ar, int n )
  {
  	int min = *ar;
  	int *minptr = ar;
  	int *ptr = ar;
  	int tmp;
  	n--;
  	ar++;
  	while( n-- )
  	{
  		if( ( tmp = *ptr++ ) < min )
  		{
  			min = tmp;
  			minptr = ptr;
  		}
  	}
  	if( minptr != ar )
  	{
  		tmp = *ar;
  		*ar = min;
  		*minptr = tmp;
  	}
  	return;
  }
  #endif
  
  void arguments( int argc, char *argv[] )
  {
      int c;
  
  	doalign = 0;
  	fromaln = 0;
  	treeout = 0;
  	uselongest = 1;
  	reorder = 1;
  	nevermemsave = 0;
  	inputfile = NULL;
  	fftkeika = 0;
  	constraint = 0;
  	nblosum = 62;
  	fmodel = 0;
  	calledByXced = 0;
  	devide = 0;
  	use_fft = 0;
  	force_fft = 0;
  	fftscore = 1;
  	fftRepeatStop = 0;
  	fftNoAnchStop = 0;
      weight = 3;
      utree = 1;
  	tbutree = 1;
      refine = 0;
      check = 1;
      cut = 0.0;
      disp = 0;
      outgap = 1;
      alg = 'A';
      mix = 0;
  	tbitr = 0;
  	scmtd = 5;
  	tbweight = 0;
  	tbrweight = 3;
  	checkC = 0;
  	treemethod = 'X';
  	contin = 0;
  	scoremtx = 1;
  	kobetsubunkatsu = 0;
  	dorp = NOTSPECIFIED;
  	ppenalty = -1530;
  	ppenalty_ex = NOTSPECIFIED;
  	poffset = -123;
  	kimuraR = NOTSPECIFIED;
  	pamN = NOTSPECIFIED;
  	geta2 = GETA2;
  	fftWinSize = NOTSPECIFIED;
  	fftThreshold = NOTSPECIFIED;
  	TMorJTT = JTT;
  	classsize = NOTSPECIFIED;
  	picksize = NOTSPECIFIED;
  	tokyoripara = NOTSPECIFIED;
  
      while( --argc > 0 && (*++argv)[0] == '-' )
  	{
          while ( ( c = *++argv[0] ) )
  		{
              switch( c )
              {
  				case 'p':
  					picksize = atoi( *++argv );
  					fprintf( stderr, "picksize = %d\n", picksize );
  					--argc;
  					goto nextoption;
  				case 's':
  					classsize = atoi( *++argv );
  					fprintf( stderr, "groupsize = %d\n", classsize );
  					--argc;
  					goto nextoption;
  				case 'i':
  					inputfile = *++argv;
  					fprintf( stderr, "inputfile = %s\n", inputfile );
  					--argc;
  					goto nextoption;
  				case 'f':
  					ppenalty = (int)( atof( *++argv ) * 1000 - 0.5 );
  //					fprintf( stderr, "ppenalty = %d\n", ppenalty );
  					--argc;
  					goto nextoption;
  				case 'g':
  					ppenalty_ex = (int)( atof( *++argv ) * 1000 - 0.5 );
  					fprintf( stderr, "ppenalty_ex = %d\n", ppenalty_ex );
  					--argc;
  					goto nextoption;
  				case 'h':
  					poffset = (int)( atof( *++argv ) * 1000 - 0.5 );
  //					fprintf( stderr, "poffset = %d\n", poffset );
  					--argc;
  					goto nextoption;
  				case 'k':
  					kimuraR = atoi( *++argv );
  					fprintf( stderr, "kimuraR = %d\n", kimuraR );
  					--argc;
  					goto nextoption;
  				case 'b':
  					nblosum = atoi( *++argv );
  					scoremtx = 1;
  //					fprintf( stderr, "blosum %d\n", nblosum );
  					--argc;
  					goto nextoption;
  				case 'j':
  					pamN = atoi( *++argv );
  					scoremtx = 0;
  					TMorJTT = JTT;
  					fprintf( stderr, "jtt %d\n", pamN );
  					--argc;
  					goto nextoption;
  				case 'm':
  					pamN = atoi( *++argv );
  					scoremtx = 0;
  					TMorJTT = TM;
  					fprintf( stderr, "tm %d\n", pamN );
  					--argc;
  					goto nextoption;
  				case 'T':
  					tokyoripara = (double)atof( *++argv );
  					--argc;
  					goto nextoption;
  				case 'l':
  					uselongest = 0;
  					break;
  #if 1
  				case 'a':
  					fmodel = 1;
  					break;
  #endif
  				case 'S':
  					doalign = 'f';
  					break;
  				case 'Z':
  					fromaln = 1;
  					break;
  				case 'L':
  					doalign = 1;
  					break;
  				case 'x':
  					reorder = 0;
  					break;
  				case 't':
  					treeout = 1;
  					break;
  				case 'r':
  					fmodel = -1;
  					break;
  				case 'D':
  					dorp = 'd';
  					break;
  				case 'P':
  					dorp = 'p';
  					break;
  				case 'e':
  					fftscore = 0;
  					break;
  				case 'O':
  					fftNoAnchStop = 1;
  					break;
  #if 0
  				case 'R':
  					fftRepeatStop = 1;
  					break;
  				case 'Q':
  					calledByXced = 1;
  					break;
  				case 'a':
  					alg = 'a';
  					break;
  #endif
  				case 'R':
  					alg = 'R';
  					break;
  				case 'Q':
  					alg = 'Q';
  					break;
  				case 'A':
  					alg = 'A';
  					break;
  				case 'N':
  					nevermemsave = 1;
  					break;
  				case 'M':
  					alg = 'M';
  					break;
  				case 'C':
  					alg = 'C';
  					break;
  				case 'F':
  					use_fft = 1;
  					break;
  				case 'G':
  					use_fft = 1;
  					force_fft = 1;
  					break;
  				case 'v':
  					tbrweight = 3;
  					break;
  				case 'd':
  					disp = 1;
  					break;
  				case 'o':
  					outgap = 0;
  					break;
  				case 'J':
  					tbutree = 0;
  					break;
  				case 'X':
  					treemethod = 'X'; // mix
  					break;
  				case 'E':
  					treemethod = 'E'; // upg (average)
  					break;
  				case 'q':
  					treemethod = 'q'; // minimum
  					break;
  				case 'z':
  					fftThreshold = atoi( *++argv );
  					--argc; 
  					goto nextoption;
  				case 'w':
  					fftWinSize = atoi( *++argv );
  					--argc;
  					goto nextoption;
                  default:
                      fprintf( stderr, "illegal option %c\n", c );
                      argc = 0;
                      break;
              }
  		}
  		nextoption:
  			;
  	}
      if( argc == 1 )
      {
          cut = atof( (*argv) );
          argc--;
      }
      if( argc != 0 ) 
      {
          fprintf( stderr, "options: Check source file !\n" );
          exit( 1 );
      }
  	if( tbitr == 1 && outgap == 0 )
  	{
  		fprintf( stderr, "conflicting options : o, m or u\n" );
  		exit( 1 );
  	}
  	if( alg == 'C' && outgap == 0 )
  	{
  		fprintf( stderr, "conflicting options : C, o\n" );
  		exit( 1 );
  	}
  }
  
  static int maxl;
  static int tsize;
  
  int seq_grp_nuc( int *grp, char *seq )
  {
  	int tmp;
  	int *grpbk = grp;
  	while( *seq )
  	{
  		tmp = amino_grp[(int)*seq++];
  		if( tmp < 4 )
  			*grp++ = tmp;
  		else
  			fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
  	}
  	*grp = END_OF_VEC;
  	return( grp-grpbk );
  }
  
  int seq_grp( int *grp, char *seq )
  {
  	int tmp;
  	int *grpbk = grp;
  	while( *seq )
  	{
  		tmp = amino_grp[(int)*seq++];
  		if( tmp < 6 )
  			*grp++ = tmp;
  		else
  			fprintf( stderr, "WARNING : Unknown character %c\r", *(seq-1) );
  	}
  	*grp = END_OF_VEC;
  	return( grp-grpbk );
  }
  
  void makecompositiontable_p( short *table, int *pointt )
  {
  	int point;
  
  	while( ( point = *pointt++ ) != END_OF_VEC )
  		table[point]++;
  }
  
  int commonsextet_p( short *table, int *pointt )
  {
  	int value = 0;
  	short tmp;
  	int point;
  	static short *memo = NULL;
  	static int *ct = NULL;
  	static int *cp;
  
  	if( !memo )
  	{
  		memo = (short *)calloc( tsize, sizeof( short ) );
  		if( !memo ) ErrorExit( "Cannot allocate memo\n" );
  		ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) );
  		if( !ct ) ErrorExit( "Cannot allocate memo\n" );
  	}
  
  	cp = ct;
  	while( ( point = *pointt++ ) != END_OF_VEC )
  	{
  		tmp = memo[point]++;
  		if( tmp < table[point] )
  			value++;
  		if( tmp == 0 ) *cp++ = point;
  	}
  	*cp = END_OF_VEC;
  	
  	cp =  ct;
  	while( *cp != END_OF_VEC )
  		memo[*cp++] = 0;
  
  	return( value );
  }
  
  void makepointtable_nuc( int *pointt, int *n )
  {
  	int point;
  	register int *p;
  
  	p = n;
  	point  = *n++ *  1024;
  	point += *n++ *   256;
  	point += *n++ *    64;
  	point += *n++ *    16;
  	point += *n++ *     4;
  	point += *n++;
  	*pointt++ = point;
  
  	while( *n != END_OF_VEC )
  	{
  		point -= *p++ * 1024;
  		point *= 4;
  		point += *n++;
  		*pointt++ = point;
  	}
  	*pointt = END_OF_VEC;
  }
  
  void makepointtable( int *pointt, int *n )
  {
  	int point;
  	register int *p;
  
  	p = n;
  	point  = *n++ *  7776;
  	point += *n++ *  1296;
  	point += *n++ *   216;
  	point += *n++ *    36;
  	point += *n++ *     6;
  	point += *n++;
  	*pointt++ = point;
  
  	while( *n != END_OF_VEC )
  	{
  		point -= *p++ * 7776;
  		point *= 6;
  		point += *n++;
  		*pointt++ = point;
  	}
  	*pointt = END_OF_VEC;
  }
  
  #if 1
  static void pairalign( int nseq, int *nlen, char **seq, int *mem1, int *mem2, double *weight, int *alloclen )
  {
  	int l, len1, len2;
  	int clus1, clus2;
  	float pscore, tscore;
  	static int *fftlog;
  	static char *indication1, *indication2;
  	static double *effarr1 = NULL;
  	static double *effarr2 = NULL;
  	static char **mseq1, **mseq2;
  	float dumfl = 0.0;
  	int ffttry;
  	int m1, m2;
  #if 0
  	int i, j;
  #endif
  
  
  	if( effarr1 == NULL ) 
  	{
  		fftlog = AllocateIntVec( nseq );
  		effarr1 = AllocateDoubleVec( nseq );
  		effarr2 = AllocateDoubleVec( nseq );
  		indication1 = AllocateCharVec( 150 );
  		indication2 = AllocateCharVec( 150 );
  		mseq1 = AllocateCharMtx( nseq, 0 );
  		mseq2 = AllocateCharMtx( nseq, 0 );
  		for( l=0; l<nseq; l++ ) fftlog[l] = 1;
  	}
  
  	tscore = 0.0;
  	m1 = mem1[0];
  	m2 = mem2[0];
  	len1 = strlen( seq[m1] );
  	len2 = strlen( seq[m2] );
  	if( *alloclen < len1 + len2 )
  	{
  		fprintf( stderr, "\nReallocating.." );
  		*alloclen = ( len1 + len2 ) + 1000;
  		ReallocateCharMtx( seq, nseq, *alloclen + 10 ); 
  		fprintf( stderr, "done. *alloclen = %d\n", *alloclen );
  	}
  
  #if WEIGHT
  	clus1 = fastconjuction_noname( mem1, seq, mseq1, effarr1, weight, indication1 );
  	clus2 = fastconjuction_noname( mem2, seq, mseq2, effarr2, weight, indication2 );
  #else
  	clus1 = fastconjuction_noweight( mem1, seq, mseq1, effarr1, indication1 );
  	clus2 = fastconjuction_noweight( mem2, seq, mseq2, effarr2, indication2 );
  #endif
  
  #if 0
  	for( i=0; i<clus1; i++ )
  		fprintf( stderr, "in p seq[%d] = %s\n", mem1[i], seq[mem1[i]] );
  	for( i=0; i<clus2; i++ )
  		fprintf( stderr, "in p seq[%d] = %s\n", mem2[i], seq[mem2[i]] );
  #endif
  
  #if 0
  	fprintf( stderr, "group1 = %.66s", indication1 );
  	if( strlen( indication1 ) > 66 ) fprintf( stderr, "..." );
  	fprintf( stderr, "\n" );
  	fprintf( stderr, "group2 = %.66s", indication2 );
  	if( strlen( indication2 ) > 66 ) fprintf( stderr, "..." );
  	fprintf( stderr, "\n" );
  #endif
  
  //	fprintf( stdout, "mseq1 = %s\n", mseq1[0] );
  //	fprintf( stdout, "mseq2 = %s\n", mseq2[0] );
  
  	if( !nevermemsave && ( alg != 'M' && ( len1 > 10000 || len2 > 10000  ) ) )
  	{
  		fprintf( stderr, "\nlen1=%d, len2=%d, Switching to the memsave mode\n", len1, len2 );
  		alg = 'M';
  		if( commonIP ) FreeIntMtx( commonIP );
  		commonAlloc1 = 0;
  		commonAlloc2 = 0;
  	}
  
  	if( fftlog[m1] && fftlog[m2] ) ffttry = ( nlen[m1] > clus1 && nlen[m2] > clus2 );
  	else	    				   ffttry = 0;
  
  	if( force_fft || ( use_fft && ffttry ) )
  	{
  		fprintf( stderr, "\bf" );
  		if( alg == 'M' )
  		{
  			fprintf( stderr, "\bm" );
  //			fprintf( stderr, "%d-%d", clus1, clus2 );
  			pscore = Falign_udpari_long( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1 );
  		}
  		else
  		{
  //			fprintf( stderr, "%d-%d", clus1, clus2 );
  			pscore = Falign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, fftlog+m1, NULL, 0, NULL );
  		}
  	}
  	else
  	{
  		fprintf( stderr, "\bd" );
  		fftlog[m1] = 0;
  		switch( alg )
  		{
  			case( 'a' ):
  				pscore = Aalign( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen );
  				break;
  			case( 'M' ):
  				fprintf( stderr, "\bm" );
  //				fprintf( stderr, "%d-%d", clus1, clus2 );
  				pscore = MSalignmm( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
  				break;
  			case( 'Q' ):
  				if( clus1 == 1 && clus2 == 1 )
  				{
  //					fprintf( stderr, "%d-%d", clus1, clus2 );
  					pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
  				}
  				else
  				{
  //					fprintf( stderr, "%d-%d", clus1, clus2 );
  					pscore = Q__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL );
  				}
  				break;
  			case( 'A' ):
  				if( clus1 == 1 && clus2 == 1 )
  				{
  //					fprintf( stderr, "%d-%d", clus1, clus2 );
  					pscore = G__align11( mseq1, mseq2, *alloclen, outgap, outgap );
  				}
  				else
  				{
  //					fprintf( stderr, "%d-%d", clus1, clus2 );
  					pscore = A__align( mseq1, mseq2, effarr1, effarr2, clus1, clus2, *alloclen, NULL, &dumfl, NULL, NULL, NULL, NULL, NULL, 0, NULL, outgap, outgap );
  				}
  				break;
  			default:
  				ErrorExit( "ERROR IN SOURCE FILE" );
  		}
  	}
  #if SCOREOUT
  	fprintf( stderr, "score = %10.2f\n", pscore );
  #endif
  	nlen[m1] = 0.5 * ( nlen[m1] + nlen[m2] );
  	return;
  }
  #endif
  
  #if 0
  static void treebase( int nseq, int *nlen, char **aseq, double *eff, int nalign, int ***topol, int *alloclen ) // topol
  {
  	int i, l;
  	int nlim;
  	int clus1, clus2;
  
  	nlim = nalign-1;
  	for( l=0; l<nlim; l++ )
  	{
  		fprintf( stderr, "in treebase, l = %d\n", l );
  		fprintf( stderr, "aseq[0] = %s\n", aseq[0] );
  		fprintf( stderr, "aseq[topol[l][0][0]] = %s\n", aseq[topol[l][0][0]] );
  		pairalign( nseq, nlen, aseq, topol[l][0], topol[l][1], eff, alloclen );
  		free( topol[l][0] );
  		free( topol[l][1] );
  	}
  }
  #endif
  
  static void WriteOptions( FILE *fp )
  {
  
  	if( dorp == 'd' ) fprintf( fp, "DNA\n" );
  	else
  	{
  		if     ( scoremtx ==  0 ) fprintf( fp, "JTT %dPAM\n", pamN );
  		else if( scoremtx ==  1 ) fprintf( fp, "BLOSUM %d\n", nblosum );
  		else if( scoremtx ==  2 ) fprintf( fp, "M-Y\n" );
  	}
      fprintf( stderr, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
      if( use_fft ) fprintf( fp, "FFT on\n" );
  
  	fprintf( fp, "tree-base method\n" );
  	if( tbrweight == 0 ) fprintf( fp, "unweighted\n" );
  	else if( tbrweight == 3 ) fprintf( fp, "clustalw-like weighting\n" );
  	if( tbitr || tbweight ) 
  	{
  		fprintf( fp, "iterate at each step\n" );
  		if( tbitr && tbrweight == 0 ) fprintf( fp, "  unweighted\n" ); 
  		if( tbitr && tbrweight == 3 ) fprintf( fp, "  reversely weighted\n" ); 
  		if( tbweight ) fprintf( fp, "  weighted\n" ); 
  		fprintf( fp, "\n" );
  	}
  
     	 fprintf( fp, "Gap Penalty = %+5.2f, %+5.2f, %+5.2f\n", (double)ppenalty/1000, (double)ppenalty_ex/1000, (double)poffset/1000 );
  
  	if( alg == 'a' )
  		fprintf( fp, "Algorithm A\n" );
  	else if( alg == 'A' ) 
  		fprintf( fp, "Algorithm A+\n" );
  	else if( alg == 'S' ) 
  		fprintf( fp, "Apgorithm S\n" );
  	else if( alg == 'C' ) 
  		fprintf( fp, "Apgorithm A+/C\n" );
  	else
  		fprintf( fp, "Unknown algorithm\n" );
  
  	if( treemethod == 'X' )
  		fprintf( fp, "Tree = UPGMA (mix).\n" );
  	else if( treemethod == 'E' )
  		fprintf( fp, "Tree = .UPGMA (average)\n" );
  	else if( treemethod == 'q' )
  		fprintf( fp, "Tree = Minimum linkage.\n" );
  	else
  		fprintf( fp, "Unknown tree.\n" );
  
      if( use_fft )
      {
          fprintf( fp, "FFT on\n" );
          if( dorp == 'd' )
              fprintf( fp, "Basis : 4 nucleotides\n" );
          else
          {
              if( fftscore )
                  fprintf( fp, "Basis : Polarity and Volume\n" );
              else
                  fprintf( fp, "Basis : 20 amino acids\n" );
          }
          fprintf( fp, "Threshold   of anchors = %d%%\n", fftThreshold );
          fprintf( fp, "window size of anchors = %dsites\n", fftWinSize );
      }
  	else
          fprintf( fp, "FFT off\n" );
  	fflush( fp );
  }
  	 
  #if 1
  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 )
  {
  	int val;
  	int ii, jj;
  	int *mptr;
  	int treelen = 0; // by Mathog, a guess
  	static int groupid = 0;
  	static int branchid = 0;
  	int i, j;
  	int selfscore0;
  	double **dfromc;
  	double **dfromcp;
  	float **pickmtx;
  	float **yukomtx;
  	static short *table1;
  	Scores **outs, *ptr;
  	int *numin;
  	int *tsukau;
  	int belongto;
  	char **children = NULL; // by Mathog, a guess
  	char *tmptree;
  	static int *orderpos = NULL;
  	int rn;
  	int npick;
  	int nyuko;
  	int *picks;
  	int *yukos;
  	int *s_p_map;
  	int *p_o_map;
  	int *s_y_map;
  	int *y_o_map;
  	int *closeh;
  	int nkouho;
  	int *pickkouho;
  	int *iptr;
  	int *jptr;
  	int aligned;
  	int ***topol;
  	int *treeorder;
  	int picktmp;
  	float **len;
  	double minscore;
  //	double *minscoreinpick;
  	float *hanni;
  	double lenfac;
  	double longer;
  	double shorter;
  	static char **mseq1 = NULL;
  	static char **mseq2 = NULL;
  	double *blastresults = NULL; // by Mathog, a guess
  	static int palloclen = 0;
  	float maxdist;
  
  	if( orderpos == NULL )
  		orderpos = order;
  	if( palloclen == 0 )
  		palloclen = *alloclen * 2;
  	if( mseq1 == NULL && doalign == 1 )
  	{
  		mseq1 = AllocateCharMtx( 1, palloclen );
  		mseq2 = AllocateCharMtx( 1, palloclen );
  	}
  
  	if( nin == 0 ) 
  	{
  #if TREE
  		if( treeout )
  		{
  			*tree = (char *)calloc( 1, sizeof( char ) );
  			**tree = 0;
  		}
  #endif
  		return 1;
  	}
  
  	if( nin < 2 || uniform == -1 ) // kokodato muda deha nai ga
  	{
  		fprintf( stderr, "\nLeaf  %d / %d                ", ++branchid, njob );
  #if 0
  		outputfile = AllocateCharVec( strlen( inputfile ) + 100 );
  		sprintf( outputfile, "%s-%d", inputfile, branchid );
  		if( uniform > 0 )
  //			sprintf( outputfile, "%su%d", outputfile, uniform );
  			sprintf( outputfile + strlen(outputfile), "u%d", uniform );
  		fprintf( stderr, "GROUP %d: %d member(s) (%d) %s\n", branchid, nin, scores[0].numinseq, outputfile );
  		outfp = fopen( outputfile, "w" );
  		free( outputfile );
  		if( outfp == NULL )
  		{
  			fprintf( stderr, "Cannot open %s\n", outputfile );
  			exit( 1 );
  		}
  		for( j=0; j<nin; j++ )
  			fprintf( outfp, ">G%d %s\n%s\n", branchid, scores[j].name+1, seq[scores[j].numinseq] );
  		fclose( outfp );
  #endif
  
  
  #if TREE
  		if( treeout )
  		{
  			treelen = 0;
  			tmptree = calloc( 100, sizeof( char ) );
  			for( j=0; j<nin; j++ )
  			{
  				treelen += sprintf( tmptree, "%d", scores[j].numinseq+1 );
  			}
  			free( tmptree );
  	
  			*tree = (char *)calloc( treelen + nin + 5, sizeof( char ) );
  			if( nin > 1 ) **tree = '(';
  			else		  **tree = '\0';
  //			**tree = '\0';
  			for( j=0; j<nin-1; j++ )
  			{
  				sprintf( *tree+strlen( *tree ), "%d,", scores[j].numinseq+1 );
  			}
  			sprintf( *tree+strlen( *tree ), "%d", scores[j].numinseq+1 );
  			if( nin > 1 ) strcat( *tree, ")" );
  //			fprintf( stdout, "*tree = %s\n", *tree );
  		}
  
  #endif
  		for( j=0; j<nin; j++ )
  		{
  			*orderpos++ = scores[j].numinseq;
  //			fprintf( stderr, "*order = %d\n", scores[j].numinseq );
  		}
  
  		return 1;
  	}
  
  
  
  	if( uselongest )
  	{
  		i = nin;
  		ptr = scores;
  		selfscore0 = scores->selfscore;
  		belongto = 0;
  		while( i-- )
  		{
  //			fprintf( stderr, "ptr-scores=%d, numinseq = %d, score = %f\n", ptr-scores, ptr->numinseq+1, ptr->score );
  			if( ptr->selfscore > selfscore0 )
  			{
  				selfscore0 = ptr->selfscore;
  				belongto = ptr-scores;
  			}
  			ptr++;
  		} 
  #if 1 
  		if( belongto != 0 )
  		{
  //			fprintf( stderr, "swap %d %s\n<->\n%d %s\n", 0, scores->name, belongto, (scores+belongto)->name );
  			ptr = calloc( 1, sizeof( Scores ) );
  			*ptr = scores[belongto];
  			scores[belongto] = *scores;
  			*scores = *ptr;
  			free( ptr );
  		}
  #endif
  	}
  	else
  	{
  		qsort( scores, nin, sizeof( Scores ), (int (*)())lcompare );
  		belongto = (int)( 0.5 * nin );
  //		fprintf( stderr, "lengths = %d, %d, %d\n", scores->orilen, scores[belongto].orilen, scores[nin-1].orilen );
  		if( belongto != 0 )
  		{
  //			fprintf( stderr, "swap %d %s\n<->\n%d %s\n", 0, scores->name, belongto, (scores+belongto)->name );
  			ptr = calloc( 1, sizeof( Scores ) );
  			*ptr = scores[belongto];
  			scores[belongto] = *scores;
  			*scores = *ptr;
  			free( ptr );
  		}
  	}
  
  	if( qinoya != scores->numinseq )
  //	if( 1 || qinoya != scores->numinseq )
  	{
  //		fprintf( stdout, "### scores->numinseq = %d, qinoya=%d, depth=%d\n", scores->numinseq, qinoya, *depthpt );
  
  
  		if( doalign )
  		{
  			if( doalign == 'f' )
  			{
  				blastresults = callfasta( seq, scores, nin, NULL, 0, 1 );
  				if( scores->selfscore != (int)blastresults[0] )
  				{
  					fprintf( stderr, "\n\nWARNING1: selfscore\n" );
  					fprintf( stderr, "scores->numinseq = %d\n", scores->numinseq+1 );
  					fprintf( stderr, "scores->orilen = %d\n", scores->orilen );
  					fprintf( stderr, "scores->selfscore = %d, but blastresults[0] = %f\n", scores->selfscore, blastresults[0] );
  //					if( abs( scores->selfscore - (int)blastresults[0] ) > 2 )
  //						exit( 1 );
  //					scores->selfscore = (int)blastresults[0]; //iinoka?
  	
  //					exit( 1 );
  				}
  			}
  			else
  				gappick0( mseq1[0], seq[scores->numinseq] );
  		}
  		else
  		{
  			table1 = (short *)calloc( tsize, sizeof( short ) );
  			if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
  			makecompositiontable_p( table1, scores[0].pointt );
  		}
  	
  		selfscore0 = scores[0].selfscore;
  		for( i=0; i<nin; i++ ) 
  		{
  			if( scores->orilen > scores[i].orilen )
  			{
  				longer = (double)scores->orilen;
  				shorter = (double)scores[i].orilen;
  			}
  			else
  			{
  				longer = (double)scores[i].orilen; // nai
  				shorter = (double)scores->orilen; //nai
  			}
  
  #if LENFAC
  			lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
  //			lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
  //			fprintf( stderr, "lenfac = %f l=%d,%d\n", lenfac,scores->orilen, scores[i].orilen );
  #else
  			lenfac = 1.0;
  #endif
  
  			if( doalign )
  			{
  				if( doalign == 'f' )
  				{
  					scores[i].score = ( 1.0 - blastresults[i] / MIN( scores->selfscore, scores[i].selfscore ) ) * 1;
  					if( scores[i].score < 0.0 ) scores[i].score = 0.0;
  				}
  				else
  				{
  					if( fromaln )
  					{
  //						scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
  						scores[i].score = ( 1.0 - (double)naivepairscore11( orialn[scores[i].numinseq], orialn[scores->numinseq], penalty ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
  					}
  					else
  					{
  						if( *depthpt == 0 ) fprintf( stderr, "\r%d / %d   ", i, nin );
  						gappick0( mseq2[0], seq[scores[i].numinseq] );
  //						fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
  						scores[i].score = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[i].selfscore ) ) * 1;
  //						fprintf( stderr, "scores[i] = %f\n", scores[i].score );
  //						fprintf( stderr, "m1=%s\n", seq[scores[0].numinseq] );
  //						fprintf( stderr, "m2=%s\n", seq[scores[i].numinseq] );
  //						fprintf( stdout, "### before calc scores[%d] = %f (%c)\n", i, scores[i].score, qinoya == scores->numinseq?'o':'x' );
  					}
  				}
  			}
  			else
  			{
  				scores[i].score = ( 1.0 - (double)commonsextet_p( table1, scores[i].pointt ) / MIN( selfscore0, scores[i].selfscore ) ) * lenfac;
  				if( scores[i].score > MAX6DIST ) scores[i].score = MAX6DIST;
  			}
  //			if( i ) fprintf( stderr, "%d-%d d %4.2f len %d %d\n", 1, i+1, scores[i].score, scores->orilen, scores[i].orilen );
  		}
  		if( doalign == 'f' ) free( blastresults );
  		if( doalign == 0 ) free( table1 );
  //exit( 1 );
  	}
  
  //	fprintf( stderr, "sorting .. " );
  	qsort( scores, nin, sizeof( Scores ), (int (*)())dcompare );
  //	fprintf( stderr, "done.\n" );
  
  
  	maxdist = scores[nin-1].score;
  	if( fromaln ) // kanzen itch ga misalign sareteiru kamoshirenai.
  	{
  		if( scores[0].shimon == scores[nin-1].shimon && !strcmp( seq[scores[0].numinseq], seq[scores[nin-1].numinseq] ) ) 
  		{
  			maxdist = 0.0;
  		}
  //		fprintf( stderr, "maxdist?? = %f, nin=%d, %d inori\n", scores[nin-1].score, nin, scores[nin-1].numinseq+1 );
  	}
  
  //	fprintf( stderr, "maxdist? = %f, nin=%d\n", scores[nin-1].score, nin );
  
  	if( nin == 1 ) fprintf( stderr, "nin=1, scores[0].score = %f\n", scores[0].score );
  
  // kokoni if( nin < 2 || ... )
  
  	picks = AllocateIntVec( nin+1 );
  	s_p_map = AllocateIntVec( nin+1 );
  	s_y_map = AllocateIntVec( nin+1 );
  	pickkouho = AllocateIntVec( nin+1 );
  	closeh = AllocateIntVec( nin+1 );
  
  //	nkouho = getkouho( pickkouho, (picksize+100)/nin, nin, scores, seq );
  //	nkouho = getkouho( pickkouho, 1.0, nin, scores, seq ); // zenbu
  //	fprintf( stderr, "selecting kouhos phase 2\n"  );
  //	if( nkouho == 0 )
  //	{
  //		fprintf( stderr, "selecting kouhos, phase 2\n"  );
  //		nkouho = getkouho( pickkouho, 1.0, nin, scores, seq );
  //	}
  //	fprintf( stderr, "\ndone\n\n"  );
  	for( i=0; i<nin; i++ ) pickkouho[i] = i+1; nkouho = nin-1; // zenbu
  
  
  
  	iptr = picks;
  	*iptr++ = 0;
  	npick = 1;
  	if( nkouho > 0 )
  	{
  //		fprintf( stderr, "pickkouho[0] = %d\n", pickkouho[0] );
  //		fprintf( stderr, "pickkouho[nin-1] = %d\n", pickkouho[nin-1] );
  		picktmp = pickkouho[nkouho-1];
  //		fprintf( stderr, "\nMOST DISTANT kouho=%d, nin=%d, nkouho=%d\n", picktmp, nin, nkouho );
  		nkouho--;
  		if( ( scores[picktmp].shimon == scores[0].shimon ) && ( !strcmp( seq[scores[0].numinseq], seq[scores[picktmp].numinseq] ) ) )
  		{
  //			fprintf( stderr, "known, j=%d (%d inori)\n", 0, scores[picks[0]].numinseq );
  //			fprintf( stderr, "%s\n%s\n", seq[scores[picktmp].numinseq], seq[scores[picks[0]].numinseq] );
  			;
  		}
  		else
  		{
  			*iptr++ = picktmp;
  			npick++;
  //			fprintf( stderr, "ok, %dth pick = %d (%d inori)\n", npick, picktmp, scores[picktmp].numinseq );
  		}
  	}
  	i = 1;
  	while( npick<picksize && nkouho>0 )
  	{
  		if( i )
  		{
  			i = 0;
  			rn = nkouho * 0.5;
  //			fprintf( stderr, "rn = %d\n", rn );
  		}
  		else
  		{
  			rn = rnd() * (nkouho);
  		}
  		picktmp = pickkouho[rn];
  //		fprintf( stderr, "rn=%d/%d (%d inori), kouho=%d, nin=%d, nkouho=%d\n", rn, nkouho, scores[pickkouho[rn]].numinseq, pickkouho[rn], nin, nkouho );
  
  //		fprintf( stderr, "#kouho before swap\n" );
  //		for( i=0; i<nkouho; i++ ) fprintf( stderr, "%d ",  pickkouho[i] ); fprintf( stderr, "\n" );
  
  		nkouho--;
  		pickkouho[rn] = pickkouho[nkouho];
  #if 1
  //		fprintf( stderr, "#kouho after swap\n" ); 
  //		for( i=0; i<nkouho; i++ ) fprintf( stderr, "%d ",  pickkouho[i] ); fprintf( stderr, "\n" );
  		for( j=0; j<npick; j++ )
  		{
  			if( scores[picktmp].shimon == scores[picks[j]].shimon && !strcmp( seq[scores[picks[j]].numinseq], seq[scores[picktmp].numinseq] ) ) 
  				break;
  		}
  		if( j == npick )
  #endif
  		{
  //			fprintf( stderr, "ok, %dth pick = %d (%d inori)\n", npick, picktmp, scores[picktmp].numinseq );
  			npick++;
  			*iptr++ = picktmp;
  		}
  		else
  		{
  //			fprintf( stderr, "known, j=%d (%d inori)\n", j, scores[picks[j]].numinseq );
  		}
  	}
  #if 0
  	for( i=0; i<nin; i++ )
  	{
  		fprintf( stderr, "i=%d/%d, scores[%d].score = %f, inori=%d\n", i, nin, i, scores[i].score, scores[i].numinseq );
  	}
  	fprintf( stderr, "range:nin=%d scores[%d].score <= %f\n", nin, npick, scores[nin-1].score);
  	for( i=0; i<npick; i++ )
  	{
  		fprintf( stderr, "i=%d/%d, scores[%d].score = %f, inori=%d\n", i, npick, picks[i], scores[picks[i]].score, scores[picks[i]].numinseq );
  	}
  exit( 1 );
  #endif
  
  //	fprintf( stderr, "\nnkouho=%d, defaultq2 = %d (%d inori)\n", nkouho, picks[npick-1], scores[picks[npick-1]].numinseq );
  
  	qsort( picks, npick, sizeof( int ), (int (*)())intcompare );
  
  //	fprintf( stderr, "allocating..\n" );
  
  //	fprintf( stderr, "allocating outs, npick = %d\n", npick );
  	numin = calloc( npick, sizeof( int ) );
  	tsukau = calloc( npick, sizeof( int ) );
  	outs = calloc( npick, sizeof( Scores * ) );
  	for( i=0; i<npick; i++ ) outs[i] = NULL;
  	topol = AllocateIntCub( npick, 2, 0 );
  	treeorder = AllocateIntVec( npick + 1 );
  	len = AllocateFloatMtx( npick, 2 );
  	pickmtx = AllocateFloatHalfMtx( npick );
  //	yukomtx = AllocateFloatHalfMtx( npick );
  //	minscoreinpick = AllocateDoubleVec( npick );
  	yukos = AllocateIntVec( npick );
  	p_o_map = AllocateIntVec( npick+1 );
  	y_o_map = AllocateIntVec( npick+1 );
  	hanni = AllocateFloatVec( npick );
  
  	for( i=0; i<nin; i++ ) s_p_map[i] = -1;
  //	fprintf( stderr, "npick = %d\n", npick );
  //	fprintf( stderr, "picks =" );
  	for( i=0; i<npick; i++ )
  	{
  		s_p_map[picks[i]] = i;
  		p_o_map[i] = scores[picks[i]].numinseq;
  //		fprintf( stderr, " %d (%dinori)\n", picks[i], scores[picks[i]].numinseq+1 );
  	}
  //	fprintf( stderr, "\n" );
  
  #if 0
  	fprintf( stderr, "p_o_map =" );
  	for( i=0; i<npick; i++ )
  	{
  		fprintf( stderr, " %d", p_o_map[i]+1 );
  	}
  	fprintf( stderr, "\n" );
  	fprintf( stderr, "picks =" );
  	for( i=0; i<npick; i++ )
  	{
  		fprintf( stderr, " %d", picks[i] );
  	}
  	fprintf( stderr, "\n" );
  #endif
  
  	for( j=0; j<nin; j++ )
  	{
  		if( s_p_map[j] != -1 )
  		{
  			pickmtx[0][s_p_map[j]] = (float)scores[j].score;
  //			fprintf( stderr, "pickmtx[0][%d] = %f\n", s_p_map[j], pickmtx[0][s_p_map[j]] );
  		}
  	}
  
  	for( j=1; j<npick; j++ )
  	{
  		if( doalign )
  		{
  			if( doalign == 'f' )
  			{
  //				blastresults = callfasta( seq, scores, npick-j+1, picks+j-1, picks[j], 1 );
  				blastresults = callfasta( seq, scores, npick, picks, picks[j], (j==1) );
  				if( scores[picks[j]].selfscore != (int)blastresults[j] )
  				{
  					fprintf( stderr, "\n\nWARNING2: selfscore j=%d/%d\n", j, npick );
  					fprintf( stderr, "scores[picks[j]].numinseq = %d\n", scores[picks[j]].numinseq+1 );
  					fprintf( stderr, "scores[picks[j]].orilen = %d\n", scores[picks[j]].orilen );
  					fprintf( stderr, "scores[picks[j]].selfscore = %d, but blastresults[j] = %f\n", scores[picks[j]].selfscore, blastresults[j] );
  //					if( abs( scores[picks[j]].selfscore - (int)blastresults[j] ) > 2 )
  //						exit( 1 );
  //					scores->selfscore = (int)blastresults[0]; //iinoka?
  				}
  			}
  			else
  				gappick0( mseq1[0], seq[scores[picks[j]].numinseq] );
  		}
  		else
  		{
  			table1 = (short *)calloc( tsize, sizeof( short ) );
  			if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
  			makecompositiontable_p( table1, scores[picks[j]].pointt );
  		}
  	
  		selfscore0 = scores[picks[j]].selfscore;
  		pickmtx[j][0] = 0.0;
  	    for( i=j+1; i<npick; i++ ) 
  		{
  			if( scores[picks[j]].orilen > scores[picks[i]].orilen )
  			{
  				longer = (double)scores[picks[j]].orilen;
  				shorter = (double)scores[picks[i]].orilen;
  			}
  			else
  			{
  				longer = (double)scores[picks[i]].orilen;
  				shorter = (double)scores[picks[j]].orilen;
  			}
  	
  	#if LENFAC
  			lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
  	//		lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
  	//		fprintf( stderr, "lenfac = %f l=%d,%d\n", lenfac,scores->orilen, scores[i].orilen );
  	#else
  			lenfac = 1.0;
  	#endif
  	
  			if( doalign )
  			{
  				if( doalign == 'f' )
  				{
  					pickmtx[j][i-j] = ( 1.0 - blastresults[i] / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
  					if( pickmtx[j][i-j] < 0.0 ) pickmtx[j][i-j] = 0.0;
  				}
  				else
  				{
  					if( fromaln )
  					{
  						fprintf( stderr, "%d-%d/%d\r", j, i, npick );
  						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;
  					}
  					else
  					{
  //						fprintf( stderr, "\r%d / %d   ", i, nin );
  						gappick0( mseq2[0], seq[scores[picks[i]].numinseq] );
  						pickmtx[j][i-j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * 1;
  	//					fprintf( stderr, "scores[picks[i]] = %f\n", scores[picks[i]].score );
  					}
  				}
  			}
  			else
  			{
  				pickmtx[j][i-j] = ( 1.0 - (double)commonsextet_p( table1, scores[picks[i]].pointt ) / MIN( selfscore0, scores[picks[i]].selfscore ) ) * lenfac;
  				if( pickmtx[j][i-j] > MAX6DIST ) pickmtx[j][i-j] = MAX6DIST;
  			}
  
  		}
  		if( doalign == 'f' ) free( blastresults );
  		if( doalign == 0 ) free( table1 );
  	}
  
  	dfromcp = AllocateDoubleMtx( npick, nin );
  	dfromc = AllocateDoubleMtx( npick, 0 );
  
  	for( i=0; i<npick; i++ ) for( j=0; j<nin; j++ )
  		dfromcp[i][j] = -0.5;
  	for( j=0; j<nin; j++ )
  	{
  		dfromcp[0][j] = ( scores[j].score );
  //		fprintf( stderr, "j=%d, s_p_map[j]=%d\n", j, s_p_map[j] );
  	}
  
  	for( i=0; i<npick; i++ ) for( j=i; j<npick; j++ )
  	{
  		dfromcp[i][picks[j]] = dfromcp[j][picks[i]] = pickmtx[i][j-i];
  	}
  
  #if 0
  	fprintf( stderr, "pickmtx = \n" );
  	for( i=0; i<npick; i++ )
  	{
  		for( j=i; j<npick; j++ )
  		{
  			fprintf( stderr, "pickmtx[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
  		}
  	}
  	exit( 1 );
  #endif
  
  
  
  //	for( i=0; i<npick-1; i++ ) for( j=i; j<npick; j++ )
  //		fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
  
  	for( i=0; i<npick; i++ ) tsukau[i] = 1;
  	for( i=0; i<nin; i++ ) closeh[i] = -1;
  	for( i=0; i<npick; i++ ) 
  	{
  		closeh[picks[i]] = picks[i];
  //		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]] );
  	}
  #if 0
  	fprintf( stderr, "closeh = \n" );
  	for( i=0; i<nin; i++ )
  	{
  		fprintf( stderr, "%d ", closeh[i] );
  	}
  	fprintf( stderr, "\n" );
  #endif
  #if DIANA
  	for( i=0; i<npick-1; i++ ) for( j=i; j<npick; j++ )
  		fprintf( stderr, "dist[%d][%d] = %f\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i] );
  	fprintf( stderr, "DIANA!!\n" );
  	if( npick > 2 )
  	{
  		float avdist;
  		float avdist1;
  		float avdist2;
  		float maxavdist;
  		int splinter;
  		int count;
  		int dochokoho;
  		splinter = 0;
  		int *docholist;
  		int *docholistbk;
  		maxavdist = 0.0;
  		for( i=0; i<npick; i++ )
  		{
  			avdist = 0.0;
  			for( j=i+1; j<npick; j++ )
  			{
  				avdist += pickmtx[i][j-i];
  			}
  			for( j=0; j<i; j++ )
  			{
  				avdist += pickmtx[j][i-j];
  			}
  			avdist /= (npick-1);
  			fprintf( stderr, "avdist[%d] = %f\n", p_o_map[i] + 1, avdist );
  			if( maxavdist < avdist ) 
  			{
  				maxavdist = avdist;
  				splinter = i;
  			}
  		}
  		fprintf( stderr, "splinter = %d (%d inori), maxavdist = %f\n", splinter, p_o_map[splinter]+1, maxavdist );
  
  		docholist = AllocateIntVec( npick );
  		docholistbk = AllocateIntVec( npick );
  		for( i=0; i<npick; i++ ) docholist[i] = 0;
  		docholist[splinter] = 1;
  		while( 1 )
  		{
  			for( i=0; i<npick; i++ ) docholistbk[i] = docholist[i]; 
  			for( dochokoho = 0; dochokoho<npick; dochokoho++ )
  			{
  				fprintf( stderr, "dochokoho=%d\n", dochokoho );
  				if( docholist[dochokoho] ) continue;
  				count = 0;
  				avdist1 = 0.0;
  				i=dochokoho;
  				{
  					for( j=i+1; j<npick; j++ )
  					{
  						if( docholist[j] || j == dochokoho ) continue;
  						avdist1 += pickmtx[i][j-i];
  						count++;
  					}
  					for( j=0; j<i; j++ )
  					{
  						if( docholist[j] || j == dochokoho ) continue;
  						avdist1 += pickmtx[j][i-j];
  						count++;
  					}
  				}
  				if( count < 1 ) avdist1 = 0.0;
  				else avdist1 /= (float)count;
  				fprintf( stderr, "docho %d (%dinori), avdist1 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist1 );
  
  				count = 0;
  				avdist2 = 0.0;
  				i=dochokoho;
  				{
  					for( j=i+1; j<npick; j++ )
  					{
  						if( !docholist[j] || j == dochokoho ) continue;
  						avdist2 += pickmtx[i][j-i];
  						count++;
  					}
  					for( j=0; j<i; j++ )
  					{
  						if( !docholist[j] || j == dochokoho ) continue;
  						avdist2 += pickmtx[j][i-j];
  						count++;
  					}
  				}
  				if( count < 1 ) avdist2 = 0.0;
  				else avdist2 /= (float)count;
  				fprintf( stderr, "docho %d (%dinori), avdist2 = %f\n", dochokoho, p_o_map[dochokoho] + 1, avdist2 );
  
  				if( avdist2 < avdist1 ) 
  				{
  					docholist[dochokoho] = 1;
  					hanni[dochokoho] = avdist2;
  				}
  				else
  				{
  					docholist[dochokoho] = 0;
  					hanni[dochokoho] = avdist1;
  				}
  				fprintf( stderr, "avdist1=%f, avdist2=%f\n", avdist1, avdist2 );
  
  			}
  			for( i=0; i<npick; i++ ) if( docholist[i] != docholistbk[i] ) break;
  			if( i == npick ) break;
  
  			fprintf( stderr, "docholist = \n" );
  			for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", docholist[i] );
  			fprintf( stderr, "\n" );
  		}
  		fprintf( stderr, "docholist = \n" );
  		for( i=0; i<npick; i++ ) fprintf( stderr, "%d ", docholist[i] );
  		fprintf( stderr, "\n" );
  
  		for( i=0; i<npick; i++ ) if( docholist[i] == 0 ) break;
  		yukos[0] = picks[i];
  		for( i=0; i<npick; i++ ) if( docholist[i] == 1 ) break;
  		yukos[1] = picks[splinter];
  
  		for( i=0; i<npick; i++ ) 
  		{
  			if( docholist[i] == 0 ) closeh[picks[i]] = yukos[0];
  			if( docholist[i] == 1 ) closeh[picks[i]] = yukos[1];
  		}
  //		for( i=0; i<npick; i++ ) closeh[picks[i]] = -1; // CHUUI !! iminai
  		nyuko = 2;
  		free( docholist );
  		free( docholistbk );
  	}
  	else if( npick > 1 )
  	{
  		nyuko = 2;
  		yukos[0] = picks[0]; yukos[1] = picks[1];
  		closeh[picks[0]] = yukos[0];
  		closeh[picks[1]] = yukos[1];
  	}
  	else
  	{
  		nyuko = 1;
  		yukos[0] = picks[0];
  		closeh[picks[0]] = yukos[0];
  	}
  #elif HUKINTOTREE
  	if( npick > 2 )
  	{
  #if 0
  		float avdist;
  		float maxavdist;
  		int count;
  		int splinter;
  		maxavdist = 0.0;
  		splinter=0;
  		for( i=0; i<npick; i++ )
  		{
  			avdist = 0.0;
  			for( j=i+1; j<npick; j++ )
  			{
  				avdist += pickmtx[i][j-i];
  			}
  			for( j=0; j<i; j++ )
  			{
  				avdist += pickmtx[j][i-j];
  			}
  			avdist /= (npick-1);
  			fprintf( stderr, "avdist[%d] = %f\n", p_o_map[i] + 1, avdist );
  			if( maxavdist < avdist ) 
  			{
  				maxavdist = avdist;
  				splinter = i;
  			}
  		}
  		fprintf( stderr, "splinter = %d (%d inori), maxavdist = %f\n", splinter, p_o_map[splinter]+1, maxavdist );
  #endif
  
  
  //		fprintf( stderr, "check kaishi =>, npick=%d members = \n", npick );
  //		for( i=0; i<npick; i++ ) fprintf( stderr, "%d (%d)", p_o_map[i]+1, picks[i] );
  //		fprintf( stderr, "\n" );
  		for( i=0; i<npick-1; i++ ) 
  		{
  			if( tsukau[i] == 0 ) continue;
  			for( j=i+1; j<npick; j++ )
  			{
  //				float kijun = maxdist *  1/(npick-2);
  //				float kijun = maxavdist * tokyoripara;
  				float kijun;
  				kijun = maxdist * tokyoripara;  // atode kakunin
  //				fprintf( stderr, "%d-%d\n", i, j );
  //				fprintf( stderr, "maxdist = %f\n", maxdist );
  //				if( i==0 && j == 1 ) continue; // machigai!! CHUUI!!
  //				if( maxdist == pickmtx[i][j-i] ) continue;
  				if( tsukau[j] == 0 ) continue;
  //				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 );
  				if( pickmtx[i][j-i] < kijun )
  				{
  //					fprintf( stderr, "dame!! %d => %d, because %f < %f\n", p_o_map[j]+1, p_o_map[i]+1, pickmtx[i][j-i], kijun );
  #if 0
  					if( scores[picks[i]].orilen > scores[picks[j]].orilen )
  					{
  						fprintf( stderr, "%d => %d\n", p_o_map[j]+1, p_o_map[i]+1 );
  						tsukau[j] = 0;
  					}
  					else
  					{
  						fprintf( stderr, "%d => %d\n", p_o_map[i]+1, p_o_map[j]+1 );
  						tsukau[i] = 0;
  					}
  					if( 0 && j == npick-1 ) tsukau[i] = 0;
  					else  			   tsukau[j] = 0;
  					fprintf( stderr, "tsukau[%d] = %d (%d inori)\n", j, tsukau[j], p_o_map[j]+1 );
  #else
  					tsukau[j] = 0;
  					closeh[picks[j]] = closeh[picks[i]];
  //					fprintf( stderr, "%d => tsukawanai\n", j );
  #endif
  				}
  			}
  		}
  	}
  	for( ii=0,i=0; i<npick; i++ )
  	{
  		if( tsukau[i] )
  		{
  			dfromc[ii] = dfromcp[i];
  			ii++;
  		}
  		else
  		{
  			free( dfromcp[i] );
  			dfromcp[i] = NULL;
  		}
  	}
  	dfromc[ii] = NULL;
  
  	for( ii=0,i=0; i<npick; i++ )
  	{
  		if( tsukau[i] )
  		{
  			for( jj=ii,j=i; j<npick; j++ )
  			{
  				if( tsukau[j] )
  				{
  					pickmtx[ii][jj-ii] = pickmtx[i][j-i];
  					jj++;
  				}
  			}
  			ii++;
  		}
  	}
  	for( ; ii<npick; ii++ )
  	{
  		free( pickmtx[ii] ); 
  		pickmtx[ii] = NULL;
  	}
  
  	for( ii=0,i=0; i<npick; i++ )
  	{
  		if( tsukau[i] )
  		{
  			yukos[ii++] = picks[i];
  		}
  	}
  
  
  	nyuko = ii;
  	yukomtx = pickmtx;
  	pickmtx = NULL;
  
  #endif
  #if 0
  	for( i=0; i<npick; i++ ) for( j=i; j<npick; j++ )
  	{
  		if( tsukau[i] == 1 && tsukau[j] == 1 )
  			fprintf( stderr, "dist[%d][%d] = %f (ok)\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i]  );
  		else if( tsukau[i] == 0 && tsukau[j] == 0 )
  			fprintf( stderr, "dist[%d][%d] = %f (xx)\n", p_o_map[i]+1, p_o_map[j]+1, pickmtx[i][j-i]  );
  		else	
  			fprintf( stderr, "%d-%d, okashii\n", p_o_map[i]+1, p_o_map[j]+1 );
  	}
  #endif
  //	FreeFloatHalfMtx( pickmtx, npick );
  
  
  	for( i=0; i<nin; i++ ) s_y_map[i] = -1;
  //	fprintf( stderr, "npick = %d\n", npick );
  //	fprintf( stderr, "yukos =" );
  	for( i=0; i<nyuko; i++ )
  	{
  		s_y_map[yukos[i]] = i;
  		y_o_map[i] = scores[yukos[i]].numinseq;
  //		fprintf( stderr, " %d\n", yukos[i] );
  	}
  //	fprintf( stderr, "\n" );
  #if 0
  	for( i=0; i<nyuko; i++ )
  	{
  		fprintf( stderr, "y_o_map[%d] = %d\n", i, y_o_map[i]+1 );
  	}
  	for( i=0; i<nyuko; i++ ) for( j=i; j<nyuko; j++ )
  	{
  		fprintf( stderr, "yukodist[%d][%d] = %f (ok)\n", y_o_map[i]+1, y_o_map[j]+1, yukomtx[i][j-i]  );
  	}
  #endif
  
  	for( i=0; i<nin; i++ )
  	{
  		if( closeh[i] != -1 )
  		{
  //			fprintf( stderr, "closeh[%d,%dinori] = %d,%dinori\n", i, scores[i].numinseq+1, closeh[i], scores[closeh[i]].numinseq+1 );
  		}
  	}
  
  #if 0
  		for( i=0; i<nyuko; i++ )
  		{
  			minscoreinpick[i] = 99.9;
  			for( j=i+1; j<nyuko; j++ )
  			{
  				if( minscoreinpick[i] > yukomtx[i][j-i] )
  					minscoreinpick[i] = yukomtx[i][j-i];
  			}
  			for( j=0; j<i; j++ )
  			{
  				if( minscoreinpick[i] > yukomtx[j][i-j] )
  					minscoreinpick[i] = yukomtx[j][i-j];
  			}
  			fprintf( stderr, "minscoreinpick[%d(%dinori)] = %f\n", i, y_o_map[i]+1, minscoreinpick[i] );
  		}
  #endif
  
  
  #if TREE
  	if( treeout )
  	{
  		children = calloc( nyuko+1, sizeof( char * ) );
  		for( i=0; i<nyuko+1; i++ ) children[i] = NULL;
  	}
  #endif
  //	fprintf( stderr, "done..\n" );
  	
  //	fprintf( stderr, "classifying, nyuko=%d \n", nyuko );
  	if( nyuko == 1 )
  	{
  		if( npick != 1 )
  		{
  			fprintf( stderr, "okashii, nyuko = 1, shikashi npick = %d\n", npick );
  			exit( 1 );
  		}
  //		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 );
  //		fprintf( stderr, "seq[%d] = scores->seq = \n%s\n", scores->numinseq, seq[scores->numinseq] );
  
  		uniform = -1;
  		for( j=0; j<nin; j++ ) 
  		{
  			belongto = 0;
  			outs[belongto] = realloc( outs[belongto], sizeof( Scores ) * ( numin[belongto] + 1 ) );
  			outs[belongto][numin[belongto]] = scores[j];
  			numin[belongto]++;
  		}
  	}
  	else
  	{
  
  #if 0
  		fprintf( stderr, "yukos = \n" );
  		for( i=0; i<nyuko; i++ ) fprintf( stderr, "%d ", y_o_map[i] + 1 );
  		fprintf( stderr, "\n" );
  #endif
  		fprintf( stderr, "\n\n%dx%d distance matrix\n", nyuko, nin );
  
  		for( i=1; i<nyuko; i++ )
  		{
  			fprintf( stderr, "%d / %d \r", i, nyuko );
  
  			if( doalign )
  			{
  				if( doalign == 'f' )
  				{
  					blastresults = callfasta( seq, scores, nin, NULL, yukos[i], (i==1) );
  				}
  				else
  					gappick0( mseq1[0], seq[scores[yukos[i]].numinseq] );
  			}
  			else
  			{
  				table1 = (short *)calloc( tsize, sizeof( short ) );
  				if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
  				makecompositiontable_p( table1, scores[yukos[i]].pointt );
  			}
  		
  			selfscore0 = scores[yukos[i]].selfscore;
  			for( j=0; j<nin; j++ ) 
  			{
  				if( scores[yukos[i]].orilen > scores[j].orilen )
  				{
  					longer = scores[yukos[i]].orilen;
  					shorter = scores[j].orilen;
  				}
  				else
  				{
  					shorter = scores[yukos[i]].orilen;
  					longer = scores[j].orilen;
  				}
  
  #if LENFAC
  //				lenfac = 1.0 / ( (double)LENFACA + (double)LENFACB / ( (double)longer + (double)LENFACC ) + (double)shorter / (double)longer * LENFACD );
  				lenfac = 1.0 / ( shorter / longer * lenfacd + lenfacb / ( longer + lenfacc ) + lenfaca );
  //				lenfac = 1.0 / ( shorter / longer * LENFACD + LENFACB / ( longer + LENFACC ) + LENFACA );
  //				fprintf( stderr, "lenfac = %f, l=%d, %d\n", lenfac, scores[yukos[i]].orilen, scores[j].orilen );
  #else
  				lenfac = 1.0;
  #endif
  #if 0 // iihazu -> dame
  				ii = s_y_map[j]; jj=s_y_map[yukos[i]];
  				if( ii != -1 && jj != -1 )
  				{
  					if( dfromc[ii][yukos[jj]] != -0.5 )
  					{
  						dfromc[i][j] = dfromc[ii][yukos[jj]];
  					}
  					else
  					{
  						if( ii > jj )
  						{
  							kk = jj;
  							jj = ii;
  							ii = kk;
  						}
  						dfromc[ii][yukos[jj]] = 
  						dfromc[i][j] = yukomtx[ii][jj-ii];
  					}
  				}
  				else
  #else
  				if( dfromc[i][j] == -0.5 )
  #endif
  				{
  					if( doalign )
  					{
  						if( doalign == 'f' )
  						{
  							dfromc[i][j] = 
  							( 1.0 - blastresults[j] / MIN( selfscore0, scores[j].selfscore ) ) * 1;
  							if( dfromc[i][j] < 0.0 ) dfromc[i][j] = 0.0;
  						}
  						else
  						{
  							if( fromaln )
  							{
  								dfromc[i][j] = ( 1.0 - (double)naivepairscore11( orialn[scores[j].numinseq], orialn[scores[yukos[i]].numinseq], penalty ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
  							}
  							else
  							{
  								gappick0( mseq2[0], seq[scores[j].numinseq] );
  								dfromc[i][j] = ( 1.0 - (double)G__align11_noalign( amino_disLN, -1200, -60, mseq1, mseq2, palloclen ) / MIN( selfscore0, scores[j].selfscore ) ) * 1;
  							}
  						}
  					}
  					else
  					{
  						dfromc[i][j] = ( 1.0 - (double)commonsextet_p( table1, scores[j].pointt ) / MIN( selfscore0, scores[j].selfscore ) ) * lenfac;
  						if( dfromc[i][j] > MAX6DIST ) dfromc[i][j] = MAX6DIST;
  					}
  				}
  //				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] );
  
  //				if( i == 1 )
  //					fprintf( stdout, "&&& dfromc[%d][%d] (%d,%d) = %f\n", i, j, p_o_map[i], scores[j].numinseq, dfromc[i][j] );
  			}
  //			fprintf( stderr, "i=%d, freeing\n", i );
  			if( !doalign ) free( table1 );
  			if( doalign && doalign == 'f' ) free( blastresults );
  		}
  		fprintf( stderr, "                \r" );
  
  
  
  
  		for( i=0; i<nyuko; i++ ) numin[i] = 0;
  //		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 );
  		for( j=0; j<nin; j++ ) 
  		{
  #if 0
  			belongto = s_y_map[j];
  			{
  				fprintf( stderr, "belongto = %d (%dinori)\n", belongto, y_o_map[belongto]+1 ); 
  			}
  			if( belongto == -1 && closeh[j] != -1 )
  #endif
  #if 0
  			if( closeh[j] != -1 )
  			{
  				belongto = s_y_map[closeh[j]];
  //				if( belongto != -1 )
  //					fprintf( stderr, "known, %d(%dinori)->%d(%dinori)\n", j, scores[j].numinseq+1, belongto, y_o_map[belongto]+1 );
  			}
  			else
  //			if( belongto == -1 )
  #else
  			belongto = s_y_map[j];
  			if( belongto == -1 )
  #endif
  			{
  				belongto = 0;  // default ha horyu
  				minscore = dfromc[0][j];
  				for( i=0; i<nyuko; i++ )
  				{
  //					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] );
  					if( scores[j].shimon == scores[yukos[i]].shimon && !strcmp( seq[scores[j].numinseq], seq[y_o_map[i]] ) ) 
  					{
  //						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 );
  						belongto = i;
  						break;
  					}
  					if( dfromc[i][j] < minscore )
  //					if( dfromc[i][j] < minscore && minscore-dfromc[i][j] > ( minscoreinpick[yukos[i]] + minscoreinpick[j] ) * 1.0 )
  //					if( rnd() < 0.5 ) // CHUUI !!!!!
  					{
  //						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] );
  						minscore = dfromc[i][j];
  						belongto = i;
  					}
  				}
  			}
  #if 0
  			if( dfromc[belongto][j] > minscoreinpick[belongto] )
  			{
  				fprintf( stderr, "dame, %f > %f\n", dfromc[belongto][j], minscoreinpick[belongto] );
  				belongto = npick;
  			}
  			else
  				fprintf( stderr, "ok, %f < %f\n", dfromc[belongto][j], minscoreinpick[belongto] );
  #endif
  //			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] );
  //			fprintf( stderr, "numin = %d\n", numin[belongto] );
  			outs[belongto] = realloc( outs[belongto], sizeof( Scores ) * ( numin[belongto] + 1 ) );
  			outs[belongto][numin[belongto]] = scores[j];
  			numin[belongto]++;
  
  		}
  		free( dfromcp );
  		FreeDoubleMtx( dfromc );
  
  //		fprintf( stderr, "##### npick = %d\n", npick );
  //		fprintf( stderr, "##### nyuko = %d\n", nyuko );
  
  
  		if( nyuko > 2 )
  		{
  			fprintf( stderr, "upgma       " );
  //			veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
  			fixed_musclesupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len, NULL );
  			fprintf( stderr, "\r                      \r" );
  		}
  		else
  		{
  			topol[0][0] = (int *)realloc( topol[0][0], 2 * sizeof( int ) );
  			topol[0][1] = (int *)realloc( topol[0][1], 2 * sizeof( int ) );
  			topol[0][0][0] = 0;
  			topol[0][0][1] = -1;
  			topol[0][1][0] = 1;
  			topol[0][1][1] = -1;
  		}
  		FreeFloatHalfMtx( yukomtx, npick );
  
  #if 0
  		ii = nyuko-1;
  		fprintf( stderr, "nyuko = %d, topol[][] = \n", nyuko );
  		for( j=0; j<nyuko-1; j++ )
  		{
  			fprintf( stderr, "STEP%d \n", j );
  			for( i=0; ; i++ )
  			{
  				fprintf( stderr, "%d ", ( topol[j][0][i] )+0 );
  				if( topol[j][0][i] == -1 ) break;
  			}
  			fprintf( stderr, "\n" );
  			for( i=0; ; i++ )
  			{
  				fprintf( stderr, "%d ", ( topol[j][1][i] )+0 );
  				if( topol[j][1][i] == -1 ) break;
  			}
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "\n" );
  		}
  		exit( 1 );
  #endif
  		jptr = treeorder; 
  		iptr = topol[nyuko-2][0]; while( *iptr != -1 ) *jptr++ = *iptr++;
  		iptr = topol[nyuko-2][1]; while( *iptr != -1 ) *jptr++ = *iptr++;
  		*jptr++ = -1;
  		for( j=0; j<nyuko; j++ )
  		{
  //			fprintf( stderr, "treeorder[%d] = %d\n", j, treeorder[j] );
  			if( treeorder[j] == -1 ) break;
  		}
  	}
  
  	aligned = 1;
  	for( i=0; i<nyuko; i++ )
  	{
  		ii = treeorder[i];
  #if 0
  		if( numin[ii] > 1 )
  		{
  			fprintf( stderr, "\ncalling a child, pick%d (%d inori): # of mem=%d\n", i, p_o_map[ii]+1, numin[ii] );
  			for( j=0; j<numin[ii]; j++ )
  			{
  				fprintf( stderr, "%d ", outs[ii][j].numinseq+1 );
  			}
  			fprintf( stderr, "\n" );
  		}
  #endif
  		aligned *= splitseq_mq( outs[ii], numin[ii], nlen, seq, orialn, name, inputfile, uniform, children+ii, alloclen, order, whichgroup, weight, depthpt, scores->numinseq );
  	}
  
  
  	for( i=0; i<nyuko; i++ )
  	{
  		if( !numin[i] )
  		{
  			fprintf( stderr, "i=%d/%d, ERROR!\n", i, nyuko );
  			for( j=0; j<nyuko; j++ )
  				fprintf( stderr, "numin[%d] = %d (rep=%d inori)\n", j, numin[j], y_o_map[j] );
  			exit( 1 );
  		}
  	}
  
  #if TREE
  	if( treeout )
  	{
  		treelen = 0;
  		for( i=0; i<nyuko; i++ )
  			treelen += strlen( children[i] );
  		*tree = calloc( treelen + nin * 3, sizeof ( char ) );
  	}
  #endif
  
  
  	if( nin >= classsize || !aligned )
  		val = 0;
  	else
  		val = 1;
  
  	if( nyuko > 1 )
  	{
  		int *mem1p, *mem2p;
  		int mem1size, mem2size;
  		int v1 = 0, v2 = 0, v3 = 0;
  		int nlim;
  		int l;
  		static int *mem1 = NULL;
  		static int *mem2 = NULL;
  		char **parttree = NULL; // by Mathog
  
  #if TREE
  		if( treeout )
  		{
  			parttree = (char **)calloc( nyuko, sizeof( char * ) );
  			for( i=0; i<nyuko; i++ )
  			{
  //				fprintf( stderr, "allocating parttree, size = %d\n", treelen + nin * 5 );
  				parttree[i] = calloc( treelen + nin * 5, sizeof ( char ) );
  				strcpy( parttree[i], children[i] );
  				free( children[i] );
  			}
  			free( children );
  		}
  #endif
  		if( mem1 == NULL )
  		{
  			mem1 = AllocateIntVec( njob+1 );
  			mem2 = AllocateIntVec( njob+1 );
  		}
  
  //		veryfastsupg_float_realloc_nobk_halfmtx( nyuko, yukomtx, topol, len );
  	
  //		counteff_simple_float( nyuko, topol, len, eff );
  
  
  		nlim = nyuko-1;
  		for( l=0; l<nlim; l++ )
  		{
  			mem1p = topol[l][0];
  			mptr = mem1;
  			mem1size = 0;
  			while( *mem1p != -1 )
  			{
  //				fprintf( stderr, "*mem1p = %d (%d inori), numin[]=%d\n", *mem1p, p_o_map[*mem1p], numin[*mem1p] );
  				i = numin[*mem1p]; ptr = outs[*(mem1p++)];
  				mem1size += i;
  				while( i-- )
  				{
  					*mptr++ = (ptr++)->numinseq;
  				}
  			}
  			*mptr = -1;
  
  			mem2p = topol[l][1];
  			mptr = mem2;
  			mem2size = 0;
  			while( *mem2p != -1 )
  			{
  //				fprintf( stderr, "*mem2p = %d (%d inori), numin[]=%d\n", *mem2p, p_o_map[*mem2p], numin[*mem2p] );
  				i = numin[*mem2p]; ptr = outs[*(mem2p++)];
  				mem2size += i;
  				while( i-- )
  				{
  					*mptr++ = (ptr++)->numinseq;
  				}
  			}
  			*mptr = -1;
  
  			qsort( mem1, mem1size, sizeof( int ), (int (*)())intcompare );
  			qsort( mem2, mem2size, sizeof( int ), (int (*)())intcompare );
  //			selhead( mem1, numin[0] );
  //			selhead( mem2, numin[1] );
  
  
  #if 0
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "mem1 (nin=%d) = \n", nin );
  			for( i=0; ; i++ )
  			{
  				fprintf( stderr, "%d ", mem1[i]+1 );
  				if( mem1[i] == -1 ) break;
  			}
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "mem2 (nin=%d) = \n", nin );
  			for( i=0; ; i++ )
  			{
  				fprintf( stderr, "%d ", mem2[i]+1 );
  				if( mem2[i] == -1 ) break;
  			}
  			fprintf( stderr, "\n" );
  #endif
  
  #if 0
  			fprintf( stderr, "before pairalign, l = %d, nyuko=%d, mem1size=%d, mem2size=%d\n", l, nyuko, mem1size, mem2size );
  			fprintf( stderr, "before alignment\n" );
  			for( j=0; j<mem1size; j++ )
  				fprintf( stderr, "%s\n", seq[mem1[j]] );
  			fprintf( stderr, "----\n" );
  			for( j=0; j<mem2size; j++ )
  				fprintf( stderr, "%s\n", seq[mem2[j]] );
  			fprintf( stderr, "----\n\n" );
  #endif
  
  			if( val )
  			{
  				fprintf( stderr, "\r  Alignment %d-%d                                 \r", mem1size, mem2size );
  				if( *mem1 < *mem2 )
  					pairalign( njob, nlen, seq, mem1, mem2, weight, alloclen );
  				else
  					pairalign( njob, nlen, seq, mem2, mem1, weight, alloclen );
  			}
  
  #if TREE
  			if( treeout )
  			{
  				v1 = topol[l][0][0];
  				v2 = topol[l][1][0];
  	
  //				fprintf( stderr, "nyuko=%d, v1=%d, v2=%d\n", nyuko, v1, v2 );
  				if( v1 > v2 )
  				{
  					v3 = v1;
  					v1 = v2;
  					v2 = v3;
  				}
  //				fprintf( stderr, "nyuko=%d, v1=%d, v2=%d after sort\n", nyuko, v1, v2 );
  //				fprintf( stderr, "nyuko=%d, v1=%d, v2=%d\n", nyuko, v1, v2 );
  //				fprintf( stderr, "v1=%d, v2=%d, parttree[v1]=%s, parttree[v2]=%s\n", v1, v2, parttree[v1], parttree[v2] );
  				sprintf( *tree, "(%s,%s)", parttree[v1], parttree[v2] );
  				strcpy( parttree[v1], *tree );
  //				fprintf( stderr, "parttree[%d] = %s\n", v1, parttree[v1] );
  //				fprintf( stderr, "*tree = %s\n", *tree );
  				free( parttree[v2] ); parttree[v2] = NULL;
  			}
  #endif
  
  #if 0
  			fprintf( stderr, "after alignment\n" );
  			for( j=0; j<mem1size; j++ )
  				fprintf( stderr, "%s\n", seq[mem1[j]] );
  			fprintf( stderr, "----\n" );
  			for( j=0; j<mem2size; j++ )
  				fprintf( stderr, "%s\n", seq[mem2[j]] );
  			fprintf( stderr, "----\n\n" );
  #endif
  		}
  #if TREE
  		if( treeout )
  		{
  			free( parttree[v1] ); parttree[v1] = NULL;
  //			fprintf( stderr, "*tree = %s\n", *tree );
  //			FreeCharMtx( parttree );
  			free( parttree ); parttree = NULL;
  		}
  #endif
  
  #if 0
  		fprintf( stderr, "after alignment\n" );
  		for( i=0; i<nyuko; i++ )
  		{
  			for( j=0; j<numin[i]; j++ )
  				fprintf( stderr, "%s\n", seq[outs[i][j].numinseq] );
  		}
  #endif
  
  		if( val )
  		{
  			groupid++;
  			mptr = mem1; while( *mptr != -1 ) 
  			{
  #if 0
  				fprintf( stdout, "==g1-%d \n", *mptr+1 );
  				fprintf( stdout, "%s \n", seq[*mptr] );
  #endif
  				whichgroup[*mptr] = groupid;
  				weight[*mptr++] *= 0.5;
  			}
  	
  			mptr = mem2; while( *mptr != -1 ) 
  			{
  #if 0
  				fprintf( stdout, "=g2-%d ", *mptr+1 );
  				fprintf( stdout, "%s \n", seq[*mptr] );
  #endif
  				whichgroup[*mptr] = groupid;
  				weight[*mptr++] *= 0.5;
  			}
  	
  			if( numin[1] == 0 )
  			{
  				mptr = mem1; while( *mptr != -1 ) 
  				{
  					whichgroup[*mptr] = groupid;
  					weight[*mptr++] *= (double)2.0/numin[0];
  				}
  			}
  		}
  		{
  			if( *depthpt > maxdepth ) maxdepth = *depthpt;
  			(*depthpt)++;
  		}
  	}
  	else
  	{
  #if TREE
  		if( treeout )
  		{
  			sprintf( *tree, "%s", children[0] );
  			free( children[0] );
  			free( children );
  		}
  #endif
  	}
  	for( i=0; i<npick; i++ ) free( (void *)outs[i] );
  //	FreeFloatHalfMtx( pickmtx, npick );
  //	FreeFloatHalfMtx( yukomtx, npick );
  	FreeFloatMtx( len );
  	FreeIntCub( topol );
  	FreeIntVec( treeorder );
  	free( outs );
  	free( numin );
  	free( picks );
  	free( yukos );
  	free( s_p_map );
  	free( s_y_map );
  	free( p_o_map );
  	free( y_o_map );
  	free( hanni );
  	free( closeh );
  	free( pickkouho );
  	free( tsukau );
  //	free( minscoreinpick );
  	return val;
  }
  #endif
  
  static void alignparaphiles( int nseq, int *nlen, double *weight, char **seq, int nmem, int *members, int *alloclen )
  {
  	int i, ilim;
  	int *mem1 = AllocateIntVec( nmem );
  	int *mem2 = AllocateIntVec( 2 );
  
  	mem2[1] = -1;
  	ilim = nmem-1;
  	for( i=0; i<ilim; i++ )
  	{
  		mem1[i] = members[i];
  		mem1[i+1] = -1;
  		mem2[0] = members[i+1];
  		pairalign( nseq, nlen, seq, mem1, mem2, weight, alloclen );
  	}
  	free( mem1 );
  	free( mem2 );
  }
  
  
  
  
  
  
  
  
  
  
  
  int main( int argc, char *argv[] )
  {
  	static char **name, **seq, **orialn;
  	static int *grpseq;
  	static char *tmpseq;
  	static int  **pointt;
  	static int *nlen;
  	int i, st, en;
  	FILE *infp;
  	FILE *treefp;
  	char *treefile = NULL; //by Mathog
  	char c;
  	int alloclen;
  	static int *order;
  	static int *whichgroup;
  	static double *weight;
  	static char tmpname[B+100];
  	int groupnum;
  	int groupid;
  	int pos;
  	int *paramem;
  	int npara;
  	int completed;
  	int orilen;
  	int pscore;
  	char *pt;
  	int **tmpaminodis;
  	static char com[1000];
  	int depth;
  	int aan;
  
  	static Scores *scores;
  	static short *table1;
  	static char **tree;
  
  
  
  	arguments( argc, argv );
  
  	if( inputfile )
  	{
  		infp = fopen( inputfile, "r" );
  		if( !infp )
  		{
  			fprintf( stderr, "Cannot open %s\n", inputfile );
  			exit( 1 );
  		}
  	}
  	else
  		infp = stdin;
  
  	getnumlen( infp );
  	rewind( infp );
  
  	if( njob < 2 )
  	{
  		fprintf( stderr, "At least 2 sequences should be input!\n"
  						 "Only %d sequence found.\n", njob ); 
  		exit( 1 );
  	}
  
  
  	if( picksize == NOTSPECIFIED || picksize < 2 )
  		picksize = PICKSIZE;
  
  	if( classsize == NOTSPECIFIED || classsize < 0 )
  	{
  		classsize = njob + 1;
  	}
  	else
  	{
  //		picksize = MIN( picksize, (int)sqrt( classsize ) + 1);
  	}
  
  	if( tokyoripara == NOTSPECIFIED )
  	{
  		if( doalign ) tokyoripara = TOKYORIPARA_A;
  		else tokyoripara = TOKYORIPARA;
  	}
  	if( picksize > njob )
  		tokyoripara = 0.0;
  
  
  	alloclen = nlenmax * 2;
  	name = AllocateCharMtx( njob, B+1 );
  
  	if( classsize == 1 )
  		seq = AllocateCharMtx( njob, 0 );
  	else
  		seq = AllocateCharMtx( njob, alloclen+1 );
  
  
  	nlen = AllocateIntVec( njob ); 
  	tmpseq = calloc( nlenmax+1, sizeof( char )  );
  	pointt = AllocateIntMtx( njob, 0 );
  	grpseq = AllocateIntVec( nlenmax + 1 );
  	order = (int *)calloc( njob + 1, sizeof( int ) );
  	whichgroup = (int *)calloc( njob, sizeof( int ) );
  	weight = (double *)calloc( njob, sizeof( double ) );
  
  	fprintf( stderr, "alloclen = %d in main\n", alloclen );
  
  	for( i=0; i<njob; i++ ) whichgroup[i] = 0;
  	for( i=0; i<njob; i++ ) weight[i] = 1.0;
  	for( i=0; i<njob; i++ ) order[i] = -1;
  
  	if( classsize == 1 )
  		readData_varlen( infp, name, nlen, seq );
  	else
  		readData_pointer( infp, name, nlen, seq );
  
  	fclose( infp );
  
  	if( fromaln ) doalign = 1;
  
  	if( fromaln )
  	{
  		orialn = AllocateCharMtx( njob, alloclen+1 );
  		for( i=0; i<njob; i++ )
  		{
  			if( strlen( seq[i] ) != nlenmax )
  			{
  				fprintf( stderr, "Input sequences must be aligned\n" );
  				exit( 1 );
  			}
  			strcpy( orialn[i], seq[i] );
  		}
  	}
  
  	constants( njob, seq );
  
  	if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
  	else              tsize = (int)pow( 6, 6 );
  
  	if( dorp == 'd' )
  	{
  		lenfaca = DLENFACA;
  		lenfacb = DLENFACB;
  		lenfacc = DLENFACC;
  		lenfacd = DLENFACD;
  	}
  	else
  	{
  		lenfaca = PLENFACA;
  		lenfacb = PLENFACB;
  		lenfacc = PLENFACC;
  		lenfacd = PLENFACD;
  	}
  
  	maxl = 0;
  	for( i=0; i<njob; i++ ) 
  	{
  		gappick0( tmpseq, seq[i] );
  		nlen[i] = strlen( tmpseq );
  		strcpy( seq[i], tmpseq );
  		pointt[i] = AllocateIntVec( nlen[i]+1 ); // ??
  
  		if( nlen[i] < 6 )
  		{
  			fprintf( stderr, "Seq %d, too short, %d characters\n", i+1, nlen[i] );
  			fprintf( stderr, "name = %s\n", name[i] );
  			fprintf( stderr, "seq = %s\n", seq[i] );
  			exit( 1 );
  //			continue;
  		}
  		if( nlen[i] > maxl ) maxl = nlen[i];
  		if( dorp == 'd' ) /* nuc */
  		{
  			if( seq_grp_nuc( grpseq, tmpseq ) < 6 )
  			{
  				fprintf( stderr, "Seq %d, too short.\n", i+1 );
  				fprintf( stderr, "name = %s\n", name[i] );
  				fprintf( stderr, "seq = %s\n", seq[i] );
  				exit( 1 );
  //				continue;
  			}
  			makepointtable_nuc( pointt[i], grpseq );
  		}
  		else                 /* amino */
  		{
  			if( seq_grp( grpseq, tmpseq ) < 6 )
  			{
  				fprintf( stderr, "Seq %d, too short.\n", i+1 );
  				fprintf( stderr, "name = %s\n", name[i] );
  				fprintf( stderr, "seq = %s\n", seq[i] );
  				exit( 1 );
  //				continue;
  			}
  			makepointtable( pointt[i], grpseq );
  		}
  //		fprintf( stdout, ">%s\n", name[i] );
  //		fprintf( stdout, "%s\n", seq[i] );
  	}
  //	exit( 1 );
  
  #if 0
  	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
  #endif
  
  	initSignalSM();
  
  	initFiles();
  
  	WriteOptions( trap_g );
  
  	c = seqcheck( seq );
  	if( c )
  	{
  		fprintf( stderr, "Illeagal character %c\n", c );
  		exit( 1 );
  	}
  
  	pid = (int)getpid();
  	sprintf( datafile, "/tmp/data-%d", pid );
  	sprintf( queryfile, "/tmp/query-%d", pid );
  	sprintf( resultfile, "/tmp/fasta-%d", pid );
  
  	scores = (Scores *)calloc( njob, sizeof( Scores ) );
  
  //	fprintf( stderr, "\nCalculating i-i scores ... \n" );
  	for( i=0; i<njob; i++ ) 
  	{
  		orilen = strlen( seq[i] );
  		scores[i].numinseq = i; // irimasu
  		scores[i].orilen = orilen;
  	}
  
  	if( doalign == 'f' )
  	{
  		fastapath = getenv( "FASTA_4_MAFFT" );
  		if( !fastapath ) fastapath = "fasta34";
  		fprintf( stderr, "fastapath=%s\n", fastapath );
  		tmpaminodis = AllocateIntMtx( 0x80, 0x80 );
  		getfastascoremtx( tmpaminodis );
  	}
  	else
  		tmpaminodis = NULL;
  	
  	for( i=0; i<njob; i++ )
  	{
  		scores[i].pointt = pointt[i];
  		scores[i].shimon = (int)pointt[i][0] + (int)pointt[i][1] + (int)pointt[i][scores[i].orilen-6];
  		scores[i].name = name[i];
  		if( doalign )
  		{
  			fprintf( stderr, "\r %05d/%05d   ", i, njob );
  			free( scores[i].pointt );
  			if( doalign == 'f' )
  			{
  #if 0
  #define KIZAMI 100
  				int ipos = (int)( i / KIZAMI ) * KIZAMI;
  				int iposamari = i % KIZAMI;
  
  				fprintf( stderr, "%d / %d\r", i, njob );
  //				fprintf( stderr, "ipos = %d\n", ipos );
  //				fprintf( stderr, "iposamari = %d\n", iposamari );
  
  //				fprintf( stderr, " calling blast, i=%d\n", i );
  //				blastresults = callfasta( seq, scores+i, 1, 0, 1 );
  				blastresults = callfasta( seq, scores+ipos, MIN(KIZAMI,njob-ipos), NULL, iposamari, (iposamari==0)  );
  //				fprintf( stderr, "done., i=%d\n\n", i );
  				scores[i].selfscore = (int)blastresults[iposamari]; 
  #if 0
  				for( j=0; j<100; j++ )
  				{
  					fprintf( stderr, "res[%d] = %f\n", j, blastresults[j] );
  				}
  #endif
  //				fprintf( stderr, "%d->selfscore = %d\n", i, scores[i].selfscore );
  				free( blastresults );
  #else
  				pscore = 0;
  				if( scoremtx == -1 )
  				{
  					st = 1;
  					en = 0;
  					for( pt=seq[i]; *pt; pt++ )
  					{
  						if( *pt == 'u' ) *pt = 't';
  						aan = amino_n[(int)*pt];
  						if( aan<0 || aan >= 4 ) *pt = 'n';
  
  						if( *pt == 'n' ) 
  						{
  							en++;
  							if( st ) continue;
  							else pscore += tmpaminodis[(int)*pt][(int)*pt];
  						}
  						else
  						{
  							st = 0;
  							en = 0;
  							pscore += tmpaminodis[(int)*pt][(int)*pt];
  						}
  					}
  					scores[i].selfscore = pscore - en * tmpaminodis['n']['n']; 
  				}
  				else
  				{
  					st = 1;
  					en = 0;
  					for( pt=seq[i]; *pt; pt++ )
  					{
  						aan = amino_n[(int)*pt];
  						if( aan<0 || aan >= 20 ) *pt = 'X';
  						if( *pt == 'X' ) 
  						{
  							en++;
  							if( st ) continue;
  							else pscore += tmpaminodis[(int)*pt][(int)*pt];
  						}
  						else
  						{
  							st = 0;
  							en = 0;
  							pscore += tmpaminodis[(int)*pt][(int)*pt];
  						}
  					}
  					scores[i].selfscore = pscore - en * tmpaminodis['X']['X']; 
  				}
  #endif
  			}
  			else
  			{
  				pscore = 0;
  				for( pt=seq[i]; *pt; pt++ )
  				{
  					pscore += amino_dis[(int)*pt][(int)*pt];
  				}
  				scores[i].selfscore = pscore; 
  			}
  //			fprintf( stderr, "selfscore[%d] = %d\n", i+1, scores[i].selfscore );
  		}
  		else
  		{
  			table1 = (short *)calloc( tsize, sizeof( short ) );
  			if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
  			makecompositiontable_p( table1, pointt[i] );
  			scores[i].selfscore = commonsextet_p( table1, pointt[i] );
  			free( table1 );
  		}
  	}
  	if( tmpaminodis ) FreeIntMtx( tmpaminodis );
  
  	depth = 0;
  #if TREE
  	if( treeout )
  	{
  		tree = (char **)calloc( 1, sizeof( char *) );
  		*tree = NULL;
  //		splitseq_bin( scores, njob, nlen, seq, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight );
  		completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
  		treefile = (char *)calloc( strlen( inputfile ) + 10, sizeof( char ) );
  		if( inputfile )
  			sprintf( treefile, "%s.tree", inputfile );
  		else
  			sprintf( treefile, "splittbfast.tree" );
  		treefp = fopen( treefile, "w" );
  		fprintf( treefp, "%s\n", *tree );
  		fclose( treefp );
  	}
  	else
  		completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
  #else
  	completed = splitseq_mq( scores, njob, nlen, seq, orialn, name, inputfile, 0, tree, &alloclen, order, whichgroup, weight, &depth, -1 );
  #endif
  
  	fprintf( stderr, "\nDone.\n\n" );
  
  #if 1
  	groupnum = 0;
  	groupid = -1;
  	paramem = NULL;
  	npara = 0;
  	for( i=0; i<njob; i++ )
  	{
  		pos = order[i];
  		if( whichgroup[pos] != groupid )
  		{
  			groupnum++;
  			groupid = whichgroup[pos];
  		}
  		if( whichgroup[pos] )
  		{
  			if( paramem )
  			{
  				paramem[npara] = -1;
  				if( npara > 1 && classsize > 2 ) 
  				{
  					qsort( paramem, npara, sizeof( int ), (int (*)(const void *, const void*))intcompare );
  //					selhead( paramem, npara );
  					alignparaphiles( njob, nlen, weight, seq, npara, paramem, &alloclen );
  				}
  				free( paramem ); paramem = NULL; npara = 0;
  			}
  			sprintf( tmpname, "Group-%d %s", groupnum, name[pos]+1 );
  		}
  		else
  		{
  			paramem = realloc( paramem, sizeof( int) * ( npara + 2 ) );
  			paramem[npara++] = pos;
  			sprintf( tmpname, "Group-para %s", name[pos]+1 );
  		}
  		tmpname[B-1] = 0;
  		if( classsize > 1 && classsize <= njob )
  			strcpy( name[pos]+1, tmpname );
  	}
  	if( paramem )
  	{
  		paramem[npara] = -1;
  		if( npara > 1 && classsize > 2 ) 
  		{
  			qsort( paramem, npara, sizeof( int ), (int (*)(const void *, const void*))intcompare );
  //			selhead( paramem, npara );
  			alignparaphiles( njob, nlen, weight, seq, npara, paramem, &alloclen );
  		}
  		free( paramem ); paramem = NULL; npara = 0;
  	}
  #else
  	for( i=0; i<njob; i++ )
  	{
  		sprintf( tmpname, "Group-%d %s", whichgroup[i], name[i]+1 );
  		strcpy( name[i]+1, tmpname );
  	}
  #endif
  
  
  //	maketanni( name, seq,  njob, nlenmax, nlen );
  
  	fclose( trap_g );
  
  #if DEBUG
  	fprintf( stderr, "writing alignment to stdout\n" );
  #endif
  	if( reorder )
  		writeData_reorder_pointer( stdout, njob, name, nlen, seq, order );
  	else
  		writeData_pointer( stdout, njob, name, nlen, seq );
  #if IODEBUG
  	fprintf( stderr, "OSHIMAI\n" );
  #endif
  	if( classsize == 1 )
  	{
  		fprintf( stderr, "\n\n" );
  		fprintf( stderr, "----------------------------------------------------------------------------\n" );
  		fprintf( stderr, "\n" );
  		fprintf( stderr, "nseq = %d\n", njob );
  		fprintf( stderr, "groupsize = %d, picksize=%d\n", classsize, picksize );
  		fprintf( stderr, "The input sequences have been sorted so that similar sequences are close.\n" );
  		if( reorder )
  			fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
  #if TREE
  		if( treeout )
  		{
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
  		}
  //		else
  //		{
  //			fprintf( stderr, "To output guide tree,\n" );
  //			fprintf( stderr, "%% %s -t  -i %s\n", progName( argv[0] ), "inputfile" );
  //		}
  #endif
  		if( !doalign )
  		{
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
  			fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
  		}
  		fprintf( stderr, "\n" );
  		fprintf( stderr, "----------------------------------------------------------------------------\n" );
  	}
  	else if( groupnum > 1 )
  	{
  		fprintf( stderr, "\n\n" );
  		fprintf( stderr, "----------------------------------------------------------------------------\n" );
  		fprintf( stderr, "\n" );
  		fprintf( stderr, "groupsize = %d, picksize=%d\n", classsize, picksize );
  		fprintf( stderr, "The input sequences have been classified into %d groups + some paraphyletic groups\n", groupnum );
  		fprintf( stderr, "Note that the alignment is not completed.\n" );
  		if( reorder )
  			fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
  #if TREE
  		if( treeout )
  		{
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
  		}
  //		else
  //		{
  //			fprintf( stderr, "To output guide tree,\n" );
  //			fprintf( stderr, "%% %s -t  -i %s\n", progName( argv[0] ), "inputfile" );
  //		}
  #endif
  		if( !doalign )
  		{
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
  			fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
  		}
  		fprintf( stderr, "\n" );
  		fprintf( stderr, "----------------------------------------------------------------------------\n" );
  	}			
  	else
  	{
  		fprintf( stderr, "\n\n" );
  		fprintf( stderr, "----------------------------------------------------------------------------\n" );
  		fprintf( stderr, "\n" );
  		fprintf( stderr, "nseq = %d\n", njob );
  		fprintf( stderr, "groupsize = %d, partsize=%d\n", classsize, picksize );
  //		fprintf( stderr, "A single alignment containing all the input sequences has been computed.\n" );
  //		fprintf( stderr, "If the sequences are highly diverged and you feel there are too many gaps,\n" );
  //		fprintf( stderr, "please try \n" );
  //		fprintf( stderr, "%% mafft --parttree --groupsize 100 inputfile\n" );
  //		fprintf( stderr, "which classifies the sequences into several groups with <~ 100 sequences\n" );
  //		fprintf( stderr, "and performs only intra-group alignments.\n" );
  		if( reorder )
  			fprintf( stderr, "The order of sequences has been changed according to estimated similarity.\n" );
  #if TREE
  		if( treeout )
  		{
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "A guide tree is in the '%s' file.\n", treefile );
  		}
  //		else
  //		{
  //			fprintf( stderr, "To output guide tree,\n" );
  //			fprintf( stderr, "%% %s -t  -i %s\n", progName( argv[0] ), "inputfile" );
  //		}
  #endif
  		if( !doalign || fromaln )
  		{
  			fprintf( stderr, "\n" );
  			fprintf( stderr, "mafft --dpparttree might give a better result, although slow.\n" );
  			fprintf( stderr, "mafft --fastaparttree is also available if you have fasta34.\n" );
  		}
  		fprintf( stderr, "\n" );
  		fprintf( stderr, "----------------------------------------------------------------------------\n" );
  	}
  #if TREE
  	if( treeout ) free( treefile );
  #endif
  
  #if 0
  	fprintf( stdout, "weight =\n" );
  	for( i=0; i<njob; i++ )
  		fprintf( stdout, "%d: %f\n", i+1, weight[i] );
  #endif
  
  	if( doalign == 'f' )
  	{
  		strcpy( com, "rm -f" );
  		strcat( com, " " );
  		strcat( com, datafile );
  		strcat( com, "*  " );
  		strcat( com, queryfile );
  		strcat( com, " " );
  		strcat( com, resultfile );
  		fprintf( stderr, "%s\n", com );
  		system( com );
  	}
  
  	SHOWVERSION;
  
  	return( 0 );
  }
  