#include "mltaln.h"

#define DEBUG 0
#define IODEBUG 0
#define SCOREOUT 0

#define GLOBAL 0

#define END_OF_VEC -1

//int nadd;
double thresholdtorev;
int dodp;
int addfragment;
int mode = '2';
int reflim = 1000;
int contrastsort = 1;

typedef struct _thread_arg
{
	int iend; 
	char **seq;
	int *map;
	char *tmpseq;
	int *res;
	int **spointt;
	short *table1;
	int iq;
#ifdef enablemultithread
	int *jshare;
	int thread_no;
	pthread_mutex_t *mutex_counter;
#endif
} thread_arg_t;

typedef struct _selfdpthread_arg
{
	int iend;
	char **seq;
	double *res;
#ifdef enablemultithread
	int *jshare;
	int thread_no;
	pthread_mutex_t *mutex_counter;
#endif
} selfdpthread_arg_t;

typedef struct _contrast
{
	int pos; 
	double dif;
} contrastarr;

static void	*selfdpthread( void *arg )
{
	selfdpthread_arg_t *targ = (selfdpthread_arg_t *)arg;
	int iend = targ->iend;
	char **seq = targ->seq;
	double *res = targ->res;
#ifdef enablemultithread
	int thread_no = targ->thread_no;
	int *jshare = targ->jshare; 
#endif
	int j;
	char **revseq;

	revseq = AllocateCharMtx( 1, nlenmax+1 );

	j = -1;
	while( 1 )
	{
#ifdef enablemultithread
		if( nthread )
		{
			pthread_mutex_lock( targ->mutex_counter );
			j = *jshare;
			if( j%100 == 0 ) reporterr( "%d / %d (thread %d)   \r", j, iend, thread_no );
			if( j == iend )
			{
				pthread_mutex_unlock( targ->mutex_counter );
				break;
			}
			++(*jshare);
			pthread_mutex_unlock( targ->mutex_counter );
		}
		else
#endif
		{
			j++;
			if( j%100 == 0 ) reporterr( "%d / %d      \r", j, iend );
			if( j == iend ) 
			{
				break;
			}
		}

		sreverse( revseq[0], seq[j] );
#if GLOBAL
		res[j] =  G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, seq+j, seq+j, 0 );
		res[j] -= G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, seq+j, revseq, 0 );
#else
		res[j] =  L__align11_noalign( n_dis_consweight_multi, seq+j, seq+j );
		res[j] -= L__align11_noalign( n_dis_consweight_multi, seq+j, revseq );
#endif
	}

	creverse( 0 );
	FreeCharMtx( revseq );
#if GLOBAL
	G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
#else
	L__align11_noalign( NULL, NULL, NULL );
#endif
	return( NULL );
}

#if 0
static void partshuffle( int size, int outsize, int *ary )
{
	int i;

//	reporterr( "ary before shuffle = \n" );
 //   for(i=0;i<size;i++) reporterr( "%d ", ary[i] );
//	reporterr( "\n" );

    for(i=0;i<outsize;i++)
    {
        int j = rand()%size;
        int t = ary[i];
        ary[i] = ary[j];
        ary[j] = t;
    }

//	reporterr( "ary after shuffle = \n" );
 //   for(i=0;i<outsize;i++) reporterr( "%d ", ary[i] );
//	reporterr( "|" );
 //   for(i=outsize;i<size;i++) reporterr( "%d ", ary[i] );
//	reporterr( "\n" );
}
#endif

void arguments( int argc, char *argv[] )
{
    int c;

	nthread = 1;
	inputfile = NULL;
	nadd = 0;
	dodp = 0;
	alg = 'a';
	alg = 'm';
	dorp = NOTSPECIFIED;
	fmodel = 0;
//	ppenalty = (int)( -2.0 * 1000 - 0.5 );
//	ppenalty_ex = (int)( -0.1 * 1000 - 0.5 );
//	poffset = (int)( 0.1 * 1000 - 0.5 ); 
	ppenalty = NOTSPECIFIED;
	ppenalty_ex = NOTSPECIFIED;
	poffset = NOTSPECIFIED;
	kimuraR = 2;
	pamN = 200;
	thresholdtorev = 0.0;
	addfragment = 0;


    while( --argc > 0 && (*++argv)[0] == '-' )
	{
        while ( (c = *++argv[0]) )
		{
            switch( c )
            {
				case 'i':
					inputfile = *++argv;
					fprintf( stderr, "inputfile = %s\n", inputfile );
					--argc;
					goto nextoption;
				case 'I':
					nadd = myatoi( *++argv );
					fprintf( stderr, "nadd = %d\n", nadd );
					--argc;
					goto nextoption;
				case 'C':
					nthread = myatoi( *++argv );
					fprintf( stderr, "nthread = %d\n", nthread );
					--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 = myatoi( *++argv );
					fprintf( stderr, "kappa = %d\n", kimuraR );
					--argc;
					goto nextoption;
				case 'j':
					pamN = myatoi( *++argv );
					scoremtx = 0;
					TMorJTT = JTT;
					fprintf( stderr, "jtt/kimura %d\n", pamN );
					--argc;
					goto nextoption;
				case 't':
					thresholdtorev = atof( *++argv );
					fprintf( stderr, "thresholdtorev = %f\n", thresholdtorev );
					--argc; 
					goto nextoption;
				case 'o':
					mode = *(*++argv);
					fprintf( stderr, "mode = %c\n", mode );
					--argc; 
					goto nextoption;
				case 'r':
					reflim = myatoi(*++argv);
					fprintf( stderr, "reflim = %d\n", reflim );
					--argc; 
					goto nextoption;
				case 'c':
					contrastsort = 0;
					break;
				case 'd':
					dodp = 1;
					break;
				case 'F':
					addfragment = 1;
					break;
#if 1
				case 'a':
					fmodel = 1;
					break;
#endif
				case 'S':
					alg = 'S';
					break;
				case 'M':
					alg = 'M';
					break;
				case 'm':
					alg = 'm';
					break;
				case 'G':
					alg = 'G';
					break;
				case 'D':
					dorp = 'd';
					break;
				case 'P':
					dorp = 'p';
					break;
                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 );
	}
}





void 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;
	if( grp - grpbk < 6 )
	{
//		fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
//		exit( 1 );
		*grpbk = -1;
	}
}

void 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;
	if( grp - grpbk < 6 )
	{
//		fprintf( stderr, "\n\nWARNING: Too short.\nPlease also consider use mafft-ginsi, mafft-linsi or mafft-ginsi.\n\n\n" );
//		exit( 1 );
		*grpbk = -1;
	}
}

void makecompositiontable_p( short *table, int *pointt )
{
	int point;

	while( ( point = *pointt++ ) != END_OF_VEC )
		table[point]++;
}


void makepointtable_nuc( int *pointt, int *n )
{
	int point;
	register int *p;

	if( *n == -1 )
	{
		*pointt = -1;
		return;
	}

	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;

	if( *n == -1 )
	{
		*pointt = -1;
		return;
	}

	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;
}

static int localcommonsextet_p2( short *table, int *pointt )
{
	int value = 0;
	short tmp;
	int point;
	short *memo;
	int *ct;
	int *cp;

	if( *pointt == -1 )
		return( 0 );

	memo = (short *)calloc( tsize, sizeof( short ) );
	if( !memo ) ErrorExit( "Cannot allocate memo\n" );
	ct = (int *)calloc( MIN( maxl, tsize )+1, sizeof( int ) ); // chuui!!
	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;

	free( memo );
	free( ct );
	return( value );
}

static int compfunc( const void *a, const void *b )
{
	return ((contrastarr *)b)->dif - ((contrastarr *)a)->dif; // correct
//	return ((contrastarr *)a)->dif - ((contrastarr *)b)->dif; // incorrect!
} 

static void makecontrastorder6mer( int *order, int **pointt, int **pointt_rev, char **seq, int iend, int shift )
{
	int i;
	double *res;
	contrastarr *arr;
	short *table1, *table1_rev;


	arr = calloc( iend, sizeof( contrastarr ) );
	res = calloc( iend, sizeof( double ) );

	for( i=0; i<iend; i++ )
	{
		if( i % 100 == 1 ) reporterr( "%d   \r", i );
		table1 = (short *)calloc( tsize, sizeof( short ) );
		if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
		makecompositiontable_p( table1, pointt[i] );
		res[i] = localcommonsextet_p2( table1, pointt[i] );
		free( table1 );

		table1_rev = (short *)calloc( tsize, sizeof( short ) );
		if( !table1_rev ) ErrorExit( "Cannot allocate table1\n" );
		makecompositiontable_p( table1_rev, pointt_rev[i] );
		res[i] -= localcommonsextet_p2( table1_rev, pointt[i] );
		free( table1_rev );

	}

	for( i=0; i<iend; i++ )
	{
		arr[i].pos = i;
		arr[i].dif = res[i];
	}

	qsort( arr, iend, sizeof( contrastarr ), compfunc );

	for( i=0; i<iend; i++ )
		order[i] = arr[i].pos + shift;

//	for( i=0; i<iend; i++ ) reporterr( "%f\n", arr[i].dif );
//	reporterr( "highest contrast, %s\n", seq[order[0]] );
//	reporterr( "lowest contrast, %s\n", seq[order[iend-1]] );

	free( arr );
	free( res );

}
static void makecontrastorder( int *order, char **seq, int iend, int shift )
{
	int i;
	double *res;
	contrastarr *arr;

	arr = calloc( iend, sizeof( contrastarr ) );
	res = calloc( iend, sizeof( double ) );

#ifdef enablemultithread
	if( nthread )
	{
		int j;
		pthread_t *handle;
		pthread_mutex_t mutex_counter;
		selfdpthread_arg_t *targ;
		int *jsharept;
		
		targ = calloc( nthread, sizeof( selfdpthread_arg_t ) );
		handle = calloc( nthread, sizeof( pthread_t ) );
		pthread_mutex_init( &mutex_counter, NULL );
		jsharept = calloc( 1, sizeof(int) );
		*jsharept = 0;
		
		for( j=0; j<nthread; j++ )
		{
			targ[j].iend = iend;
			targ[j].seq = seq;
			targ[j].res = res; 
			targ[j].jshare = jsharept;
			targ[j].mutex_counter = &mutex_counter;
			targ[j].thread_no = j;
			pthread_create( handle+j, NULL, selfdpthread, (void *)(targ+j) );
		}
		for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
		pthread_mutex_destroy( &mutex_counter );
		free( handle );
		free( targ );
		free( jsharept );
	}
	else
#endif
	{
		selfdpthread_arg_t *targ;
		targ = calloc( 1, sizeof( selfdpthread_arg_t ) );
		targ[0].iend = iend;
		targ[0].seq = seq;
		targ[0].res = res; 
		selfdpthread( targ );
		free( targ );
	}

	for( i=0; i<iend; i++ )
	{
		arr[i].pos = i;
		arr[i].dif = res[i];
	}

	qsort( arr, iend, sizeof( contrastarr ), compfunc );

	for( i=0; i<iend; i++ )
		order[i] = arr[i].pos + shift;

//	for( i=0; i<iend; i++ ) reporterr( "%f\n", arr[i].dif );
//	reporterr( "highest contrast, %s\n", seq[order[0]] );
//	reporterr( "lowest contrast, %s\n", seq[order[iend-1]] );

	free( arr );
	free( res );

}


static void	*directionthread( void *arg )
{
	thread_arg_t *targ = (thread_arg_t *)arg;
	int iend = targ->iend;
	char **seq = targ->seq;
	int *map = targ->map;
	char *tmpseq = targ->tmpseq;
	int *res = targ->res;
	int **spointt = targ->spointt;
	short *table1 = targ->table1;
//	int iq = targ->iq;
#ifdef enablemultithread
//	int thread_no = targ->thread_no;
	int *jshare = targ->jshare; 
#endif
	int j;
	char **mseq1, **mseq2;


	if( dodp ) // nakuserukamo
	{
		mseq1 = AllocateCharMtx( 1, 0 );
		mseq2 = AllocateCharMtx( 1, 0 );
	}

	j = -1;
	while( 1 )
	{
#ifdef enablemultithread
		if( nthread )
		{
			pthread_mutex_lock( targ->mutex_counter );
			j = *jshare;
			if( j == iend )
			{
				pthread_mutex_unlock( targ->mutex_counter );
				break;
			}
			++(*jshare);
			pthread_mutex_unlock( targ->mutex_counter );
		}
		else
#endif
		{
			j++;
			if( j == iend ) 
			{
//				if( iq%100==1 ) fprintf( stderr, "\r %d / %d  \r", iq, njob );
				break;
			}
		}


		if( dodp )
		{
//			strcpy( mseq1[0], tmpseq );
//			strcpy( mseq2[0], seq[j] );
			mseq1[0] = tmpseq;
			mseq2[0] = seq[map[j]];
#if GLOBAL
			res[j] = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1, mseq2, 0 );
#else
			res[j] = L__align11_noalign( n_dis_consweight_multi, mseq1, mseq2 );
#endif
		}
		else
		{
//			reporterr( "\n\nj=%d, map[j]=%d\n\n", j, map[j] );
			res[j] = localcommonsextet_p2( table1, spointt[map[j]] );
		}
	}
	if( dodp ) // nakuserukamo
	{
		free( mseq1 );
		free( mseq2 );
#if GLOBAL
		G__align11_noalign( NULL, 0, 0, NULL, NULL, 0 );
#else
		L__align11_noalign( NULL, NULL, NULL );
#endif
	}
//	else
//		if( nthread )  // inthread == 0 no toki free suru to, error. nazeda
//			localcommonsextet_p( NULL, NULL );
	return( NULL );
}

int main( int argc, char *argv[] )
{
	static int  *nlen;	
	static int  *nogaplen;	
	static char **name, **seq;
	int i, j, istart, iend, ic;
	FILE *infp;
//	FILE *adfp;
	char c;

	int *grpseq;
	char *tmpseq, *revseq;
	int  **pointt, **pointt_rev, **spointt;
	double res_forward, res_reverse, res_max;
	int ires, mres, mres2;
	int *res, *resr, *resf;
	int *map;
	static short *table1, *table1_rev;
	static char **mseq1f, **mseq1r, **mseq2;
	int *contrastorder;

	arguments( argc, argv );
#ifndef enablemultithread
	nthread = 0;
#endif

	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( alg == 'a' )
	{
		if( nlenmax < 10000 )
			alg = 'G';
		else
			alg = 'S';
	}

	seq = AllocateCharMtx( njob, nlenmax*1+1 );

#if 0
	Read( name, nlen, seq );
	readData( infp, name, nlen, seq );
#else
    name = AllocateCharMtx( njob, B+1 );
    nlen = AllocateIntVec( njob ); 
    nogaplen = AllocateIntVec( njob ); 
	readData_pointer( infp, name, nlen, seq );
	fclose( infp );

	if( dorp != 'd' )
	{
		fprintf( stderr, "Not necessary!\n" );
		for( i=0; i<njob; i++ ) 
			fprintf( stdout, "_F_%-10.10s\n", name[i]+1 );
		exit( 1 );
	}
#endif

	constants( njob, seq );


#if 0
	fprintf( stderr, "params = %d, %d, %d\n", penalty, penalty_ex, offset );
#endif

	initSignalSM();

	initFiles();

	c = seqcheck( seq );
	if( c )
	{
		fprintf( stderr, "Illegal character %c\n", c );
		exit( 1 );
	}

	fprintf( stderr, "\n" );
	if( alg == 'G' ) // dp to the first sequence
	{
		mseq1f = AllocateCharMtx( 1, nlenmax+nlenmax );
		mseq1r = AllocateCharMtx( 1, nlenmax+nlenmax );
		mseq2 = AllocateCharMtx( 1, nlenmax+nlenmax );
	    tmpseq = AllocateCharVec( MAX( nlenmax, B ) +1 );

		gappick0( mseq1f[0], seq[0] );
		sreverse( mseq1r[0], mseq1f[0] );
		strcpy( seq[0], mseq1f[0] );

		if( nadd )
			istart = njob - nadd;
		else
			istart = 1;

		fprintf( stderr, "\n" );

		for( i=0; i<istart; i++ )
		{
			gappick0( tmpseq, seq[i] );
			strcpy( seq[i], tmpseq );
			strcpy( tmpseq, name[i] );
			strcpy( name[i], "_F_" );
			strncpy( name[i]+3, tmpseq+1, 10 );
			name[i][13] = 0;
		}
		for( i=istart; i<njob; i++ ) 
		{
			gappick0( mseq2[0], seq[i] );

#if GLOBAL
			res_forward = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1f, mseq2, 0 );
			res_reverse = G__align11_noalign( n_dis_consweight_multi, penalty, penalty_ex, mseq1r, mseq2, 0 );
#else
			res_forward = L__align11_noalign( n_dis_consweight_multi, mseq1f, mseq2 );
			res_reverse = L__align11_noalign( n_dis_consweight_multi, mseq1r, mseq2 );
#endif
#if 0

			strcpy( mseq2[0], seq[i] );
			strcpy( mseq1f[0], seq[0] );
			res_forward = G__align11( n_dis_consweight_multi, mseq1f, mseq2, nlenmax*2, 0, 0 );
			fprintf( stdout, "%s\n", mseq1f[0] );
			fprintf( stdout, "%s\n", mseq2[0] );

			strcpy( mseq2[0], seq[i] );
			sreverse( mseq1r[0], seq[0] );
			res_reverse = G__align11( n_dis_consweight_multi, mseq1r, mseq2, nlenmax*2, 0, 0 );
			fprintf( stdout, "%s\n", mseq1r[0] );
			fprintf( stdout, "%s\n", mseq2[0] );
#endif

//			fprintf( stdout, "\nscore_for(%d,%d) = %f\n", 0, i, res_forward );
//			fprintf( stdout, "score_rev(%d,%d) = %f\n", 0, i, res_reverse );
			res_max = MAX(res_reverse,res_forward);
			if( (res_reverse-res_forward)/res_max > thresholdtorev ) // tekitou
			{
//				fprintf( stderr, "REVERSE!!!\n" );
				sreverse( seq[i], mseq2[0] );

				strcpy( tmpseq, name[i] );
				strcpy( name[i], "_R_" );
				strncpy( name[i]+3, tmpseq+1, 10 );
				name[i][13] = 0;
			}
			else
			{
				strcpy( seq[i], mseq2[0] );

				strcpy( tmpseq, name[i] );
				strcpy( name[i], "_F_" );
				strncpy( name[i]+3, tmpseq+1, 10 );
				name[i][13] = 0;
			}
		}
		FreeCharMtx( mseq1f );
		FreeCharMtx( mseq1r );
		FreeCharMtx( mseq2 );
		free( tmpseq );
	}
	else if( alg == 'm' )
	{

		if( dodp ) // nakuserukamo
		{
			mseq1f = AllocateCharMtx( 1, nlenmax+1);
			mseq1r = AllocateCharMtx( 1, nlenmax+1 );
			mseq2 = AllocateCharMtx( 1, nlenmax+1 );
		}
		else
		{
//			nthread = 0; // heiretsu keisan no kouritsu ha warui node
			spointt = AllocateIntMtx( njob, 0 ); 
			pointt = AllocateIntMtx( njob, nlenmax+1 );
			pointt_rev = AllocateIntMtx( njob, nlenmax+1 );
		}
	    tmpseq = AllocateCharVec( MAX( nlenmax, B ) +1 );
	    revseq = AllocateCharVec( nlenmax+1 );
		grpseq = AllocateIntVec( nlenmax+1 );
		res = AllocateIntVec( njob );
		resr = AllocateIntVec( njob );
		resf = AllocateIntVec( njob );
		map = AllocateIntVec( njob );
		contrastorder = AllocateIntVec( njob );
		if( dorp == 'd' ) tsize = (int)pow( 4, 6 );
		else              tsize = (int)pow( 6, 6 ); // iranai

		maxl = 0;
		for( i=0; i<njob; i++ )
		{
			gappick0( tmpseq, seq[i] );
			nogaplen[i] = strlen( tmpseq );
			if( nogaplen[i] > maxl ) maxl = nogaplen[i];
		}

		reporterr( "Step 1/2\n" );


		if( !dodp )
		{
			if( nadd )
				iend = njob - nadd;
			else
				iend = 0; // keisan shinai
	
			for( i=0; i<iend; i++ )
			{
				gappick0( tmpseq, seq[i] );
				strcpy( seq[i], tmpseq );
				seq_grp_nuc( grpseq, tmpseq );
				makepointtable_nuc( pointt[i], grpseq );
				spointt[i] = pointt[i];
			}
	
			if( nadd )
				istart = njob - nadd;
			else
				istart = 0;
			for( i=istart; i<njob; i++ ) 
			{
				gappick0( tmpseq, seq[i] );
				strcpy( seq[i], tmpseq );
				sreverse( revseq, tmpseq );

				seq_grp_nuc( grpseq, tmpseq );
				makepointtable_nuc( pointt[i], grpseq );
//				makecompositiontable_p( table1, pointt[i] ); -> moto no basho ni modosu
				seq_grp_nuc( grpseq, revseq );
				makepointtable_nuc( pointt_rev[i], grpseq );
//				makecompositiontable_p( table1_rev, pointt_rev[i] ); -> moto no basho ni modosu
				spointt[i] = pointt[i];


//				reporterr( "pointt[i] = %p\n", pointt[i] );
//				reporterr( "pointt[i][0] = %p\n", pointt[i][0] );

			}	
		}


		if( contrastsort ) // sukoshi chuui
		{


			if( nadd )
			{
				iend = njob-nadd;
				for( i=0; i<iend; i++ ) contrastorder[i] = i;
				istart = njob-nadd;
				iend = nadd;
			}
			else
			{
				istart = 0;
				iend = njob;
			}

			if( dodp )
				makecontrastorder( contrastorder+istart, seq+istart, iend, istart );
			else
				makecontrastorder6mer( contrastorder+istart, pointt+istart, pointt_rev+istart, seq+istart, iend, istart );
		}
		else
		{
			for( i=0; i<njob; i++ ) contrastorder[i] = i;
		}


//		reporterr( "contrastorder = \n" );
//		for( i=0; i<njob; i++ )
//			reporterr( "%d ", contrastorder[i] );
//		reporterr( "\n" );



		if( nadd )
			iend = njob - nadd;
		else
			iend = 1;
		for( i=0; i<iend; i++ )
		{
			ic = contrastorder[i];
//			fprintf( stdout, "%d, SKIP\n", i );
			gappick0( tmpseq, seq[ic] );
			strcpy( seq[ic], tmpseq );
//			if( !nadd ) strcpy( seq[i], tmpseq ); // seq ha tsukawanaikara ii.

#if 0 // -> makecontrastorder() no mae ni idou
			if( !dodp )
			{
				seq_grp_nuc( grpseq, tmpseq );
				makepointtable_nuc( pointt[ic], grpseq );
				spointt[ic] = pointt[ic];
			}
#endif

			strcpy( tmpseq, name[ic] );
			strcpy( name[ic], "_F_" );
			strncpy( name[ic]+3, tmpseq+1, 10 );
			name[ic][13] = 0;
		}

		reporterr( "\n\nStep 2/2\n" );

		if( nadd )
			istart = njob - nadd;
		else
			istart = 1;
		for( i=istart; i<njob; i++ ) 
		{
//			fprintf( stderr, "\r %d / %d ", i, njob );
			ic = contrastorder[i];
			gappick0( tmpseq, seq[ic] );
			strcpy( seq[ic], tmpseq );
			sreverse( revseq, tmpseq );

#if 0 // -> makecontrastorder() no mae ni idou
			if( !dodp )
			{
				table1 = (short *)calloc( tsize, sizeof( short ) );
				if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
				table1_rev = (short *)calloc( tsize, sizeof( short ) );
				if( !table1_rev ) ErrorExit( "Cannot allocate table1_rev\n" );
				seq_grp_nuc( grpseq, tmpseq );
				makepointtable_nuc( pointt[ic], grpseq );
				makecompositiontable_p( table1, pointt[ic] );
				seq_grp_nuc( grpseq, revseq );
				makepointtable_nuc( pointt_rev[ic], grpseq );
				makecompositiontable_p( table1_rev, pointt_rev[ic] );
			}
#else
			if( !dodp )
			{
				table1 = (short *)calloc( tsize, sizeof( short ) );
				if( !table1 ) ErrorExit( "Cannot allocate table1\n" );
				table1_rev = (short *)calloc( tsize, sizeof( short ) );
				if( !table1_rev ) ErrorExit( "Cannot allocate table1_rev\n" );
				makecompositiontable_p( table1, pointt[ic] );
				makecompositiontable_p( table1_rev, pointt_rev[ic] );
			}
#endif

			if( nadd && addfragment )
				iend = njob-nadd;
			else
				iend = i;


			if( iend > reflim ) 
			{
//				reporterr( "iend = %d -> %d\n", iend, reflim );
#if 0
				for( j=0; j<iend; j++ ) map[j] = j;
				partshuffle( iend, reflim, map );
#else
				for( j=0; j<iend; j++ ) map[j] = contrastorder[j];
#endif
				iend = reflim; // approximation
			}
			else
			{
#if 0
				for( j=0; j<iend; j++ ) map[j] = j;
#else
				for( j=0; j<iend; j++ ) map[j] = contrastorder[j];
#endif
			}

//			reporterr( "reflim = %d, seq[%d] = %s\n", reflim, contrastorder[0], seq[contrastorder[0]] );

#ifdef enablemultithread
			if( nthread )
			{
				pthread_t *handle;
				pthread_mutex_t mutex_counter;
				thread_arg_t *targ;
				int *jsharept;
		
				targ = calloc( nthread, sizeof( thread_arg_t ) );
				handle = calloc( nthread, sizeof( pthread_t ) );
				pthread_mutex_init( &mutex_counter, NULL );
				jsharept = calloc( 1, sizeof(int) );
				*jsharept = 0;
		
				if( i%100==1 ) fprintf( stderr, " %d / %d (%d threads)   \r", i, njob, nthread );
				for( j=0; j<nthread; j++ )
				{
					targ[j].iend = iend;
					targ[j].map = map;
					targ[j].seq = seq;
					targ[j].tmpseq = tmpseq; 
					targ[j].res = resf; 
					targ[j].spointt = spointt; 
					targ[j].table1 = table1; 
					targ[j].jshare = jsharept;
					targ[j].iq = i; // iranai
					targ[j].mutex_counter = &mutex_counter;
					targ[j].thread_no = j;
					pthread_create( handle+j, NULL, directionthread, (void *)(targ+j) );
				}
				for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
				pthread_mutex_destroy( &mutex_counter );
				free( handle );
				free( targ );
				free( jsharept );
			}
			else
#endif
			{
				thread_arg_t *targ;

				if( i%100==1 ) fprintf( stderr, " %d / %d   \r", i, njob );
				targ = calloc( 1, sizeof( thread_arg_t ) );
				targ[0].iend = iend;
				targ[0].map = map;
				targ[0].seq = seq;
				targ[0].tmpseq = tmpseq; 
				targ[0].res = resf; 
				targ[0].spointt = spointt; 
				targ[0].table1 = table1; 
				targ[0].iq = i;  // iranai
				directionthread( targ );
				free( targ );
			}



#ifdef enablemultithread
			if( nthread )
			{
				pthread_t *handle;
				pthread_mutex_t mutex_counter;
				thread_arg_t *targ;
				int *jsharept;
		
				targ = calloc( nthread, sizeof( thread_arg_t ) );
				handle = calloc( nthread, sizeof( pthread_t ) );
				pthread_mutex_init( &mutex_counter, NULL );
				jsharept = calloc( 1, sizeof(int) );
				*jsharept = 0;
		
				for( j=0; j<nthread; j++ )
				{
					targ[j].iend = iend;
					targ[j].seq = seq;
					targ[j].map = map;
					targ[j].tmpseq = revseq; 
					targ[j].res = resr; 
					targ[j].spointt = spointt; 
					targ[j].table1 = table1_rev; 
					targ[j].jshare = jsharept;
					targ[j].iq = i; // iranai
					targ[j].mutex_counter = &mutex_counter;
					targ[j].thread_no = j;
					pthread_create( handle+j, NULL, directionthread, (void *)(targ+j) );
				}
				for( j=0; j<nthread; j++ ) pthread_join( handle[j], NULL );
				pthread_mutex_destroy( &mutex_counter );
				free( handle );
				free( targ );
				free( jsharept );
			}
			else
#endif
			{
				thread_arg_t *targ;
				targ = calloc( 1, sizeof( thread_arg_t ) );
				targ[0].iend = iend;
				targ[0].seq = seq;
				targ[0].map = map;
				targ[0].tmpseq = revseq; 
				targ[0].res = resr; 
				targ[0].spointt = spointt;
				targ[0].table1 = table1_rev; 
				targ[0].iq = i;  // iranai
				directionthread( targ );
				free( targ );
			}

			if( mode == '2' )
			{
				mres = mres2 = 0;
				for( j=0; j<iend; j++ )
				{
					ires = resf[j];
//					fprintf( stdout, "ires (%d,%d) = %d\n", i, j, ires );
//					fflush( stdout );
					if( ires>mres2 ) 
					{
						if( ires>mres ) 
						{
							mres2 = mres;
							mres = ires;
						}
						else
							mres2 = ires;
					}
				}
				res_forward = (double)( mres + mres2 ) / 2;
				mres = mres2 = 0;
				for( j=0; j<iend; j++ )
				{
					ires = resr[j];
					if( ires>mres2 )
					{
						if( ires>mres ) 
						{
							mres2 = mres;
							mres = ires;
						}
						else
							mres2 = ires;
					}
				}
				res_reverse = (double)( mres + mres2 ) / 2;
				res_max = MAX(res_reverse,res_forward);
			}
//			reporterr( "i=%d, res_reverse = %f\n", i, res_reverse );
			else if( mode == '1' )
			{
				res_reverse = 0.0;
				for( j=0; j<iend; j++ ) if( res_reverse < (double)resr[j] ) res_reverse = (double)resr[j];
				res_forward = 0.0;
				for( j=0; j<iend; j++ ) if( res_forward < (double)resf[j] ) res_forward = (double)resf[j];
				res_max = 1.0;
			}

			else if( mode == 'd' )
			{
				res_reverse = 0.0;
				for( j=0; j<iend; j++ ) if( res_reverse < (double)(resr[j]-resf[j]) ) res_reverse = (double)(resr[j]-resf[j]);
				res_forward = 0.0;
				for( j=0; j<iend; j++ ) if( res_forward < (double)(resf[j]-resr[j]) ) res_forward = (double)(resf[j]-resr[j]);
				res_max = 1.0;
			}

			else if( mode == 'a' )
			{
				res_reverse = 0.0;
				for( j=0; j<iend; j++ ) res_reverse += (double)resr[j];
				res_reverse /= (double)iend;
				res_forward = 0.0;
				for( j=0; j<iend; j++ ) res_forward += (double)resf[j];
				res_forward /= (double)iend;
				res_max = 1.0;
			}
			else
			{
				reporterr( "Unknown mode!\n" );
				exit( 1 );
			}


			if( (res_reverse>res_forward) ) // tekitou
//			if( (res_reverse-res_forward)/res_max > thresholdtorev ) // tekitou
			{
				strcpy( seq[ic], revseq );

				strcpy( tmpseq, name[ic] );
				strcpy( name[ic], "_R_" );
				strncpy( name[ic]+3, tmpseq+1, 10 );
				name[ic][13] = 0;
				if( !dodp ) spointt[ic] = pointt_rev[ic];
			}
			else
			{
				strcpy( tmpseq, name[ic] );
				strcpy( name[ic], "_F_" );
				strncpy( name[ic]+3, tmpseq+1, 10 );
				name[ic][13] = 0;
				if( !dodp ) spointt[ic] = pointt[ic];
			}

			if( !dodp )
			{
				free( table1 );
				free( table1_rev );
			}
		}

		if( name[0][1] == 'R' )
		{
			for( j=0; j<njob; j++ ) 
			{
				if( name[j][1] == 'R' ) 
					name[j][1] = 'F';
				else
					name[j][1] = 'R';
			}
		}

		creverse( 0 );
		free( tmpseq );
		free( revseq );
		free( grpseq );
		free( res );
		free( resr );
		free( resf );
		free( map );
		free( nlen );
		free( nogaplen );
		free( contrastorder );
		if( dodp )
		{
			FreeCharMtx( mseq1f );
			FreeCharMtx( mseq1r );
			FreeCharMtx( mseq2 );
		}
		else
		{
			FreeIntMtx( pointt );
			FreeIntMtx( pointt_rev );
			free( spointt );
		}
	}
	else
	{
		fprintf( stderr, "Unknown alg %c\n", alg );
		exit( 1 );
	}
//	writeData_pointer( stdout, njob, name, nlen, seq );
	for( i=0; i<njob; i++ ) 
	{
//		fprintf( stdout, ">%s\n", name[i] );
//		fprintf( stdout, "%s\n", seq[i] );
		fprintf( stdout, "%s\n", name[i] );
	}

	FreeCharMtx( seq );
	FreeCharMtx( name );
	freeconstants();
	closeFiles();

	fprintf( stderr, "\n" );
	SHOWVERSION;
	return( 0 );
}