//lastest version, noise in considered.

#include	<stdio.h>
#include	<stdlib.h>
#include	<math.h>
#include	<string.h>
#include	<assert.h>

#define	MASS_H 	1.0079
#define	MASS_C 	12.011
#define	MASS_O 	15.9994

#define	ERROR_OF_PRECURSOR_CHARGE1 15
#define	ERROR_OF_PRECURSOR_CHARGE2 30
#define	ERROR_OF_PRECURSOR_CHARGE3 45

#define	MIN_NUM_OF_PEAKS	30
#define PEAK_FACTOR	6
#define	MAX_FRAGMENT_LENGTH 30
#define MIN_FRAGMENT_LENGTH 5
#define COEFFICIENT	2.0

#define	SIGMA	9.0

#define	THRESHOLD 4.0
//#define	DEBUG 1

//#define	PRINT_REAL_IDEAL

#define	SCALE_REGION 65536
double	scale_factor[ 2*SCALE_REGION + 1 ];

typedef struct  _DIGEST_SITE_{
        int     offset;
        char    AA;
        int     yes_or_no;
}_DIGEST_SITE;

typedef	struct	_PROBABILITY_ {
	double	offset;
	double	pr;
}_PROBABILITY_;

typedef	struct	_DTA_ {
	double	precursor;
	int	charge;
	double	sum_intensity;
	long	min_ion;
	long	max_ion;
	long	num_of_peaks;
	double 	*real_spectrum;
	char	*ID;
	int	index;
        double  best_score;
        long int        index_of_best_matched_fragment;
        int     num_of_matched_fragments;
}_DTA_;

typedef	struct	_FRAGMENT_ {
	char	*sequence;
	int	length;
	double	*ideal_spectrum_charge1;
	double	*ideal_spectrum_charge2;
	long	sum_of_pr1;
	long	sum_of_pr2;
	int	num_of_ions1;
	int	num_of_ions2;
	double	mass;
        double  best_score;
        long int    index_of_best_matched_dta;
//	double	*score;
//      long int  *index_of_matched_dta;
        int num_of_matched_dtas;
}_FRAGMENT_;

typedef	struct	_PROTEIN_ {
	int	index;
	char	*protein_name;
	char	*protein_sequence;
	long	num_of_fragments;
	struct	_FRAGMENT_ **fragment_list;
}_PROTEIN_;

double	AA_Table[256];
struct  _DIGEST_SITE_   *digest_site;
int     num_of_digest_site;
int     num_of_N1, num_of_C1, num_of_sN2, num_of_dN2, num_of_sC2, num_of_dC2;
struct  _PROBABILITY_   *Pr_N1, *Pr_C1, *Pr_sN2, *Pr_dN2, *Pr_sC2, *Pr_dC2;
double  noise_1, noise_2;

struct	_DTA_	**dta_list;
long int	num_of_dta;
struct	_PROTEIN_ **protein_list;
int	num_of_proteins;
struct	_FRAGMENT_ **fragment_library;
long int	total_num_of_fragments;

void	init_AA_table()
{
	int	i;
	
	for( i = 0; i < 256; i++ )    {
		AA_Table[i] = -1;
        }
        AA_Table['G'] = 570.21;
        AA_Table['A'] = 710.37;
        AA_Table['S'] = 870.23;
        AA_Table['P'] = 970.53;
        AA_Table['V'] = 990.68;
        AA_Table['T'] = 1010.48;
        AA_Table['C'] = 1030.09;
        AA_Table['I'] = 1130.84;
        AA_Table['L'] = 1130.84;
        AA_Table['X'] = 1130.84;
        AA_Table['N'] = 1140.43;
        AA_Table['O'] = 1140.79;
        AA_Table['B'] = 1145.35;
        AA_Table['D'] = 1150.27;
        AA_Table['Q'] = 1280.59;
        AA_Table['K'] = 1280.90;
        AA_Table['Z'] = 1285.51;
        AA_Table['E'] = 1290.43;
        AA_Table['M'] = 1310.40;
        AA_Table['H'] = 1370.59;
        AA_Table['F'] = 1470.68;
        AA_Table['R'] = 1561.00;
        AA_Table['Y'] = 1630.63;
        AA_Table['W'] = 1860.79;

        return;
}

int    initiate_probability_table( char *filename_of_probability_list )
{
        FILE    *fp;
        char    flag[256];
        double  offset, probability;
        int     result;
        int     i;

        num_of_N1 = 0;
        num_of_C1 = 0;
        num_of_sN2 = 0;
        num_of_dN2 = 0;
        num_of_sC2 = 0;
        num_of_dC2 = 0;
        
        fp = fopen( filename_of_probability_list, "r" );
        if( fp == NULL ){
		printf("Can't open file %s!\n", filename_of_probability_list );
		exit(-1);
	}

        while( !feof(fp) )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( result < 3 )
                        break;

                if( strcmp( flag, "N1" ) == 0 )
                        num_of_N1++;
                if( strcmp( flag, "C1" ) == 0 )
                        num_of_C1++;
                if( strcmp( flag, "sN2" ) == 0 )
                        num_of_sN2++;
                if( strcmp( flag, "dN2" ) == 0 )
                        num_of_dN2++;
                if( strcmp( flag, "sC2" ) == 0 )
                        num_of_sC2++;
                if( strcmp( flag, "dC2" ) == 0 )
                        num_of_dC2++;
        }
            
        Pr_N1 = malloc( num_of_N1 * sizeof( struct _PROBABILITY_ ) );
        Pr_C1 = malloc( num_of_C1 * sizeof( struct _PROBABILITY_ ) );
        Pr_sN2 = malloc( num_of_sN2 * sizeof( struct _PROBABILITY_ ) );
        Pr_dN2 = malloc( num_of_dN2 * sizeof( struct _PROBABILITY_ ) );
        Pr_sC2 = malloc( num_of_sC2 * sizeof( struct _PROBABILITY_ ) );
        Pr_dC2 = malloc( num_of_dC2 * sizeof( struct _PROBABILITY_ ) );
        if( (Pr_N1== NULL) || (Pr_C1== NULL) || (Pr_sN2== NULL) ||
            (Pr_dN2== NULL) || (Pr_sC2== NULL) || (Pr_dC2== NULL) )
        {
                printf( "Memory request failed!\n" );
                exit(-1);
        }
        
        rewind(fp);

        memset( flag, '\0', 255 );
        result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
        if( strcmp( flag, "noise1" ) == 0 )
        {
                noise_1 = probability;
        }
        else
        {
                fclose(fp);
                return -1;
        }

        for( i=0; i < num_of_N1; i++ )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( strcmp( flag, "N1" ) != 0 )
                {
                        fclose( fp );;
                        return -1;;
                }
                offset = lrint( offset*10.0 );
                Pr_N1[i].offset = offset;
                Pr_N1[i].pr = probability;
        }

        for( i=0; i < num_of_C1; i++ )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( strcmp( flag, "C1" ) != 0 )
                {
                        fclose(fp);
                        return -1;
                }
                offset = lrint( offset*10.0 );
                Pr_C1[i].offset = offset;
                Pr_C1[i].pr = probability;
        }

        memset( flag, '\0', 255 );
        result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
        if( strcmp( flag, "noise2" ) == 0 )
        {
                noise_2 = probability;
        }
        else
        {
                fclose(fp);
                return -1;
        }
        
        for( i=0; i < num_of_sN2; i++ )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( strcmp( flag, "sN2" ) != 0 )
                {
                        fclose(fp);
                        return -1;
                }
                offset = lrint( offset*10.0 );
                Pr_sN2[i].offset = offset;
                Pr_sN2[i].pr = probability;
        }

        for( i=0; i < num_of_dN2; i++ )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( strcmp( flag, "dN2" ) != 0 )
                {
                        fclose(fp);
                        return -1;
                }
                offset = lrint( offset*10.0 );
                Pr_dN2[i].offset = offset;
                Pr_dN2[i].pr = probability;
        }

        for( i=0; i < num_of_sC2; i++ )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( strcmp( flag, "sC2" ) != 0 )
                {
                        fclose(fp);
                        return -1;
                }
                offset = lrint( offset*10.0 );
                Pr_sC2[i].offset = offset;
                Pr_sC2[i].pr = probability;
        }

        for( i=0; i < num_of_dC2; i++ )
        {
                memset( flag, '\0', 255 );
                result = fscanf( fp, "%s%lf%lf", flag, &offset, &probability );
                if( strcmp( flag, "dC2" ) != 0 )
                {
                        fclose(fp);
                        return -1;
                }
                offset = lrint( offset*10.0 );
                Pr_dC2[i].offset = offset;
                Pr_dC2[i].pr = probability;
        }

        fclose(fp);
        return 1;
}

void	readin_enzyme_site_file( char *digest_file_name )
{
	FILE	*fp;
	int	count;
	int	offset;
	char	AA;

	fp = fopen( digest_file_name, "r" );
	if( fp == NULL ){
		printf("Can't open file %s!\n", digest_file_name );
		exit(-1);
	}
	count = 0;
	while( !feof(fp) ){
		int	result;
		result = fscanf( fp, "%d%c", &offset, &AA );
		if( result < 2 )
			break;
		count++;
	}
	if( count == 0 ){
		printf("Enzyme site configuration error!\n");
		exit(-1);
	}

	num_of_digest_site = count;
	digest_site = malloc( count * sizeof(struct _DIGEST_SITE_ ) );
	if( digest_site == NULL ){
		printf("Memory request failed! For digest site\n");
		exit(-1);
	}

	rewind( fp );
	count = 0;
	while( !feof(fp )){
		int	result;
		result = fscanf( fp, "%d%c", &offset, &AA );
		if( result < 2 )
			break;
		digest_site[count].offset = offset;
		if( AA >= 'a' && AA <= 'z' ){
			digest_site[count].AA = 'A' + AA - 'a';
			digest_site[count].yes_or_no = -1;
		}else{
			digest_site[count].AA = AA;
			digest_site[count].yes_or_no = 1;
		}
		count++;
	}
	fclose(fp);

    return;
}

int	readin_dta( char *dta_filename, _DTA_ *p )
{
	FILE	*fp;
	double	precursor, sum_intensity, intensity;
	double	max_ion, min_ion, ion;
	long	i, result, loc, count;
	int	charge;
	double	first_max, second_max;
	char	ch, buffer[1024];

	max_ion = min_ion = -1;
	count = 0;

	// get the max ion
	fp = fopen( dta_filename, "r" );
	if( fp == NULL ){
		printf("Can't open file %s!\n", dta_filename );
		exit (-1);
	}

	i = 0;
	while( 1 ){
		ch = fgetc( fp );
		if( ch == '\n' )
			break;
		if( ch == 0x0d ){
			ch = fgetc(fp);
			assert( ch == 0x0a );
			break;
		}
		buffer[ i++ ] = ch;
	}
	buffer[ i ] = '\0';
	sscanf( buffer, "%lf%d", &precursor, &charge );

	first_max = second_max = -1;

	while( !feof( fp ) ){
		result = fscanf( fp, "%lf%lf", &ion, &intensity );
		if( result < 2 ){
			break;
		}
		count++;
		if( ion > max_ion || max_ion < 0 )
			max_ion = ion;
		if( ion < min_ion || min_ion < 0 )
			min_ion = ion;

		if( intensity > second_max ){
			if( intensity > first_max ){
				second_max = first_max;
				first_max = intensity;
			}else{
				second_max = intensity;
			}
		}
	}

	if( count < MIN_NUM_OF_PEAKS ){
    		fclose(fp);
		return	-1;
	}
	if( first_max > PEAK_FACTOR * second_max ){
    		fclose(fp);
		return -1;
	}

	// malloc memory
	count = lrint( max_ion * 10.0 + 10 );
	p->num_of_peaks = count;

	p->real_spectrum = (double *)malloc( count * sizeof(double) );
	p->ID = (char *)malloc( strlen( dta_filename ) + 1 );
	if( p->real_spectrum == NULL || p->ID == NULL ){
		printf("Memory request failed! For real_spectrum and ID of a dta file.\n");
    		fclose(fp);
		exit(-1);
	}
	strcpy( p->ID, dta_filename );

	for( i = 0; i < count; i++ )
		p->real_spectrum[i] = 0;

	// reread the dta file
	sum_intensity = 0;

        rewind( fp );

	i = 0;
	while( 1 ){
		ch = fgetc( fp );
		if( ch == '\n' || ch == 0x0d || ch == 0x0a )
			break;
		buffer[ i++ ] = ch;
	}
	buffer[ i ] = '\0';
	sscanf( buffer, "%lf%d", &precursor, &charge );

	p->precursor = precursor * 10.0;
	p->charge = charge;
	while( !feof( fp ) ){
		long	delta;

		result = fscanf( fp, "%lf%lf", &ion, &intensity );
		if( result < 2 )
			break;
		loc = lrint( ion * 10.0 - 0.5 );
		delta = loc - lrint(p->precursor);
		if( delta <=  10 && delta >= -10 )
			continue;

		sum_intensity += intensity;
		p->real_spectrum[ loc ] = intensity;
	}
	p->sum_intensity = sum_intensity;
        //+4 or -4 is for the entension of each peak;
	p->min_ion = lrint( min_ion * 10.0 - 0.5 );
	p->max_ion = lrint( max_ion * 10.0 - 0.5 );
        p->best_score = -1;
        p->index_of_best_matched_fragment = -1;
        // initial value of "num_of_matched_fragments" must be set as 0 but not -1;
        // see "compare_ideal_real_spectrum" for details.
        p->num_of_matched_fragments = 0;

	fclose(fp);

	return	1;
}

void    sort_peak( double *peak, long int *peak_index, int left, int right )
{
        int     temp;
        double  	average;
        int 	i, j;

	if( left >= right )
		return;
	i = left;
	j = right;
	average = peak[ peak_index[(left+right)/2] ];

	do{
		while( i <= right && peak[ peak_index[i] ] < average )
			i++;
		while( j >= left  && peak[ peak_index[j] ] > average )
			j--;
		if( i <= j ){
			temp = peak_index[i];
			peak_index[i] = peak_index[j];
			peak_index[j] = temp;
			i++;
			j--;
		}
	}while( i <= j );

	if( i < right ) sort_peak( peak, peak_index, i, right );
	if( j > left ) sort_peak( peak, peak_index, left, j );

	return;
}

int	process_real_spectrum( _DTA_ *p )
{
    	long int  	*peak_index;
    	double  	threshold;
    	double  	intensity;
    	double 		sum;
    	int 		count;
        double  	average;
	int		i;

    	if ( !( p->sum_intensity > 0 ) ) {
        	printf( "%s: summary of intensity is less than 0!\n", p->ID );
        	return -1;
    	}

    	peak_index = (long int *) malloc( p->num_of_peaks * sizeof( long int ) );
	if( peak_index == NULL ){
		printf("Memory request failed!\n");
		exit(-1);
	}

	for( i = 0; i < p->num_of_peaks; i++ )
		peak_index[i] = i;

	if( p->charge == 1 )
		threshold = (noise_1/100.0) * p->sum_intensity;
	else
		threshold = (noise_2/100.0) * p->sum_intensity;

//    	sort_peak( p->real_spectrum, peak_index, 0, p->num_of_peaks-1 );

//    assert( threshold > 0 );
/****
	sum = 0.0;
    	for( i=0; i < p->num_of_peaks; i++ )
        {
        	if( sum > threshold )
            		break;
        	else
                {
            		sum += p->real_spectrum[ peak_index[i] ];
            		p->real_spectrum[ peak_index[i] ] = 0;
        	}
    	}
***/

    	for( i=p->min_ion; i <= p->max_ion; i++ )
        {
        	if( p->real_spectrum[i] > 0 )
                {
            		intensity = p->real_spectrum[i];
            		intensity /= 9;
			p->real_spectrum[i] = intensity;
			p->real_spectrum[ i+1 ] += intensity;
    			p->real_spectrum[ i+2 ] += intensity;
    			p->real_spectrum[ i+3 ] += intensity;
    			p->real_spectrum[ i+4 ] += intensity;
			p->real_spectrum[ i-1 ] += intensity;
    			p->real_spectrum[ i-2 ] += intensity;
    			p->real_spectrum[ i-3 ] += intensity;
    			p->real_spectrum[ i-4 ] += intensity;
        	}
    	}
	p->min_ion -= 4;
	p->max_ion += 4;

//	average = sum / (double)(p->max_ion - p->min_ion );
	average = threshold / (double)(p->max_ion - p->min_ion );

    	for( i= p->min_ion; i <= p->max_ion; i++ )
		p->real_spectrum[i] += average;

	for( i = p->min_ion; i <= p->max_ion; i++ )
		p->real_spectrum[i] /= p->sum_intensity;

    	free(peak_index);

	return 1;
}

int	readin_dta_list( char *dta_list_filename )
{
	FILE	*fp;
	long int	count;
    	int i, result;
	char	buff[1024];

	fp = fopen ( dta_list_filename, "r" );
	if( fp == NULL )
    	{
		printf("File %s open error!\n", dta_list_filename );
		exit(-1);
	}

	count = 0;
	while( !feof( fp ) )
    	{
		result = fscanf( fp, "%s", buff );
		if( result < 1 )
			break;
		count++;
	}

	num_of_dta = count;
	dta_list = malloc( count * sizeof( struct _DTA_ *) );
	if( dta_list == NULL )
    	{
		printf("Memory request failed! For list of dta files\n");
		exit(-1);
	}

	for( i = 0; i < count; i++ ) {
		dta_list[i] = (struct _DTA_ *)malloc( sizeof(struct _DTA_) );
		if( dta_list[i] == NULL ) {
			printf("Memory request failed! for a dta file\n");
			exit(-1);
		}
		dta_list[i]->precursor = -1;
        	dta_list[i]->charge = -1;
        	dta_list[i]->sum_intensity = -1;
        	dta_list[i]->min_ion = -1;
        	dta_list[i]->max_ion = -1;
        	dta_list[i]->num_of_peaks = -1;
		dta_list[i]->index = i;
		dta_list[i]->ID = NULL;
		dta_list[i]->real_spectrum = NULL;
        	dta_list[i]->best_score = -1;
        	dta_list[i]->index_of_best_matched_fragment = -1;
	}

	rewind( fp );

	for( i = 0; i < count; i++ ) {
		result = fscanf( fp, "%s", buff );

		if( result < 1 )
			break;
		if( readin_dta( buff, dta_list[i] ) < 0 ) {
			free( dta_list[i] );
			dta_list[i] = NULL;
			continue;
		}

		if( process_real_spectrum( dta_list[i] ) < 0 ) {
	            	if (dta_list[i]->ID != NULL) {
                		free( dta_list[i]->ID );
		                dta_list[i]->ID = NULL;
        	    	}

		   	if (dta_list[i]->real_spectrum != NULL ) {
                		free( dta_list[i]->real_spectrum );
                		dta_list[i]->real_spectrum = NULL;
            		}

			free( dta_list[i] );
            		dta_list[i] = NULL;
		}
	}

	fclose(fp);
	return 1;
}

int	readin_protein_sequence( FILE *fp, char **protein_name, char **protein_sequence )
{
	char	ch;
	int	new_line_flag = 1;
	char	*name, *sequence;
	int	limit,count;

	new_line_flag = 1;
 	// read until a CR/LF followd by '>'
	*protein_name = NULL;
	*protein_sequence = NULL;

	while( 1 ){
		ch = fgetc( fp );
		if( ch == EOF )
			return -1;
		if( new_line_flag && ch == '>' )
			break;
		if( ch == '\n' || ch == 0x0a || ch == 0x0d )
			new_line_flag = 1;
		else
			new_line_flag = 0;
	}

	// read the head line till a CR/LF
	limit = 0;
	name = NULL;
	count = 0;

	while( 1 ){
		ch = fgetc( fp );
		if( ch == '\n' || ch == 0x0a || ch == 0x0d )
			break;
		if( ch == EOF )
			return -1;

		if( count == limit ){
			limit += 1024;
			name = realloc( name, limit*sizeof(char) );
			if( name == NULL ){
				printf("Memory request failed! For protein name.\n");
				exit(-1);
			}
		}
		name[ count ] = ch;
		count++;
	}
	if( name ){
		name[ count ] = '\0';
	}

	count = 0;
	limit = 0;
	sequence = NULL;
	// read in sequence
	new_line_flag = 1;
	while( 1 ){
		ch = fgetc( fp );
		if( ch == 0x0a || ch == 0x0d || ch == '\n' ){
			new_line_flag = 1;
			continue;
		}
		if( new_line_flag && ch == '>' ){
			fseek( fp, -1L, SEEK_CUR );
			break;
		}
		new_line_flag = 0;
		if( ch == EOF )
			break;
		if( count > ( limit - 10 )){
			limit += 1024;
			sequence = realloc( sequence, limit*sizeof(char));
			if( sequence == NULL ){
				printf("Memory request failed! For protein sequence.\n");
				exit(-1);
			}
		}
		sequence[ count ] = ch;
		count++;
	}

	if( sequence ){
		sequence[ count ] = '\0';
	}
	*protein_name = name;
	*protein_sequence = sequence;
	return	1;
}

struct	_FRAGMENT_ * new_a_fragment( char *buff, int count )
{
	long	i, length;
	struct	_FRAGMENT_ *f;
	double	mass;

	length = count + 1;

	f = malloc( sizeof( struct _FRAGMENT_ ) );
	if( f == NULL ){
		printf("Memory request failed! For a new fragment.\n");
		exit(-1);
	}
	f->sequence = malloc( length * sizeof(char) );
	if( f->sequence == NULL ){
		printf("Memory request failed! For seuqence of a new fragment.\n");
		exit(-1);
	}
	memset( f->sequence, '\0', length );
	memcpy( f->sequence, buff, count );
	f->length = count;

	mass = 0;
	for( i = 0; i < count; i++ ){
		if( AA_Table[ buff[i] ] < 0 ){
			printf("Unknow AA id: %c. Check it!\n", buff[i] );
			exit(-1);
		}
		mass += AA_Table[ buff[i] ];
	}
	f->mass = mass;

	f->ideal_spectrum_charge1 = NULL;
	f->ideal_spectrum_charge2 = NULL;
	f->sum_of_pr1 = -1;
	f->sum_of_pr2 = -1;
	f->num_of_ions1 = -1;
	f->num_of_ions2 = -1;
        f->best_score = -1.0;
        f->index_of_best_matched_dta = -1;
        f->num_of_matched_dtas = -1;

        return	f;
}

int	is_a_site( char *p, char *buff )
{
	int	i;
	int	flag[4];
	int	result = 0;

	for( i = 0; i < 4; i++ )
		flag[i] = 0;

	if( p == buff )
		result++;
	if( *(p+1) == '\0' )
		return 1;
	for( i = 0; i < num_of_digest_site; i++ ){
		if( (digest_site[i].AA == '*') ||
	  	    (digest_site[i].AA == *(p + digest_site[i].offset ) )) {
			if( digest_site[i].yes_or_no == -1 )
				flag[ digest_site[i].offset ] = -1;
			else{
				if( flag[ digest_site[i].offset ] != -1 )
					flag[ digest_site[i].offset ] = digest_site[i].yes_or_no;
			}
		}
	}
	if( flag[0] == 1 && flag[1] == 1 && flag[2] == 1 )
		result++;
	return	result;
}

int	get_fragment_num( char *buff )
{
	char	*p, *end_ptr;
	int	num = 0;

	if( buff == NULL )
		return 0;
	for( p = buff; *p; p++ ){
		int	result;
		result = is_a_site( p, buff );
		if( result == 0 )
			continue;

		for( end_ptr = p + 1; *end_ptr; end_ptr++ ){
			if( !is_a_site( end_ptr, buff ) )
				continue;

			if( (end_ptr - p) > MAX_FRAGMENT_LENGTH )
				break;
			if( p == buff ){
				if( end_ptr < p + MIN_FRAGMENT_LENGTH-1 )
					continue;
			}else{
				if( end_ptr < p + MIN_FRAGMENT_LENGTH )
					continue;
			}
			num++;

			if( p == buff && result == 2 && end_ptr >= p+MIN_FRAGMENT_LENGTH )
				num++;
		}
	}
	return num;
}

int	digest_to_fragments( struct _PROTEIN_ *pro )
{
	char	*p, *end_ptr, *buff;
	int	i, n;

	n = get_fragment_num( pro->protein_sequence );
	if( n == 0 )
		return -1;
	pro->num_of_fragments = n;

	pro->fragment_list = malloc( n * sizeof( struct _FRAGMENT_ *) );
	if( pro->fragment_list == NULL ){
		printf("Memory request failed! For fragment list of a protein\n");
		return -1;
	}

	buff = pro->protein_sequence;
	i = 0;
	for( p = buff; *p; p++ ){
		int	result;

		result = is_a_site( p, buff );
		if( result == 0 )
			continue;

		for( end_ptr = p + 1; *end_ptr; end_ptr++ ){
			if( !is_a_site( end_ptr, buff ) )
				continue;
			if( (end_ptr - p) > MAX_FRAGMENT_LENGTH )
				break;
			if( p == buff ) {
				if( end_ptr < p + MIN_FRAGMENT_LENGTH-1 )
					continue;
			}else{
				if( end_ptr < p + MIN_FRAGMENT_LENGTH )
					continue;
			}

			// get a fragment
			if( p == buff ){
				pro->fragment_list[ i ] = new_a_fragment( p, (end_ptr-p)+1);
				i++;
			}else{
				pro->fragment_list[ i ] = new_a_fragment( p+1, (end_ptr-p));
				i++;
			}
			if( p == buff && result == 2 && end_ptr >= (p+MIN_FRAGMENT_LENGTH) ){
				pro->fragment_list[ i ] = new_a_fragment( p+1, (end_ptr-p));
				i++;
			}
		}
	}
	return	1;
}

int	read_and_digest_proteins( char *protein_list_filename )
{
	FILE	*fp;
	int	i, count, result;
	char	*protein_name, *protein_sequence;

	fp = fopen( protein_list_filename, "r" );
	if( fp == NULL ){
		printf("Can't open file %s\n", protein_list_filename );
		exit(-1);
	}

	count = 0;
	while( 1 ){
		result = readin_protein_sequence( fp, &protein_name, &protein_sequence );

		if( result < 1 ){
			if( protein_name != NULL )
         		       free( protein_name );
			if( protein_sequence != NULL )
	                	free( protein_sequence );
			break;
		}

		if( protein_sequence != NULL )
			count++;
		if( protein_name != NULL )
	            	free( protein_name );
		if( protein_sequence != NULL )
            		free( protein_sequence );
		protein_name = protein_sequence = NULL;
	}

	if( count <= 0 )
		return -1;

	num_of_proteins = count;

	protein_list = malloc( count * sizeof( struct _PROTEIN_ *) );
	if( protein_list == NULL ){
		printf("Memory request failed! For a list of protein.\n");
		exit(-1);
	}
	for( i = 0; i < count; i++ ){
		protein_list[i] = malloc( sizeof(struct _PROTEIN_ ) );
		if( protein_list[i] == NULL ){
			printf("Memory request failed! For a protein.\n");
			exit(-1);
		}
	}

	rewind( fp );

	for( i = 0; i < count; i++ ){
		result = readin_protein_sequence( fp, &protein_name, &protein_sequence );
		if( result < 0 )
			break;
		if( protein_sequence == NULL )
			continue;
		protein_list[i]->index = i;
		protein_list[i]->protein_name = protein_name;
		protein_list[i]->protein_sequence = protein_sequence;

		if( digest_to_fragments( protein_list[i] ) < 0 )
			return -1;
   	}

   	fclose( fp );

   	return 1;
}

int	produce_ideal_spectrum_charge1( struct _FRAGMENT_ *f, char *fragment, int count )
{
	int	i, j, index;
	double	mass;
	long	ion, mass_int;
    	double  delta_b, delta_y;
	double	pr, pr2, pr3, pr4, pr5;
	double	*pr_tmp;
	double	sum_of_pr;
    	double  sum_of_noise;
    	double  average;
	double	middle_point;
	double	delta, r;

    	delta_b  = 10.0 * MASS_H;
    	delta_y  = 10.0 * ( 3*MASS_H + MASS_O );

	mass_int = lrint( f->mass + 100 );
	middle_point = mass_int / 2.0;

	pr_tmp = malloc( mass_int * sizeof(double) );
	if( pr_tmp == NULL ){
		printf("Malloc error!\n");
		exit(-1);
	}

	for( i = 0; i < mass_int; i++ )
		pr_tmp[i] = 0;

	mass = 0.0;
	for( i = 0; i < count; i++ ){
		if( AA_Table[ fragment[i] ] < 0 ){
			printf("Unknow AA id: %c. Check it!\n", fragment[i] );
			exit(-1);
		}
		mass += AA_Table[ fragment[i] ];

		for( j = 0; j < num_of_N1; j++ ){
			ion = lrint( mass + delta_b + Pr_N1[j].offset - 0.5 );
			if( ion <= 0 || ion >= mass_int )
				continue;
			pr = Pr_N1[j].pr / 25.0;
			pr2 = pr + pr;
			pr3 = pr2 + pr;
			pr4 = pr3 + pr;
			pr5 = pr4 + pr;

			pr_tmp[ ion ] += pr5;
			pr_tmp[ ion+1 ] += pr4;
			pr_tmp[ ion-1 ] += pr4;
			pr_tmp[ ion+2 ] += pr3;
			pr_tmp[ ion-2 ] += pr3;
			pr_tmp[ ion+3 ] += pr2;
			pr_tmp[ ion-3 ] += pr2;
			pr_tmp[ ion+4 ] += pr;
			pr_tmp[ ion-4 ] += pr;
		}
	}

	mass = 0.0;
	for( i = count-1; i > 0; i-- ){
		if( AA_Table[ fragment[i] ] < 0 ){
			printf("Unknow AA id: %c. Check it!\n", fragment[i] );
			exit(-1);
		}
		mass += AA_Table[ fragment[i] ];
		for( j = 0; j < num_of_C1; j++ ){
			ion = lrint( mass + delta_y + Pr_C1[j].offset - 0.5 );
			if( ion <= 0 || ion >= mass_int )
				continue;

			pr = Pr_C1[j].pr / 25.0;
			pr2 = pr + pr;
			pr3 = pr2 + pr;
			pr4 = pr3 + pr;
			pr5 = pr4 + pr;

			pr_tmp[ ion ] += pr5;
			pr_tmp[ ion+1 ] += pr4;
			pr_tmp[ ion-1 ] += pr4;
			pr_tmp[ ion+2 ] += pr3;
			pr_tmp[ ion-2 ] += pr3;
			pr_tmp[ ion+3 ] += pr2;
			pr_tmp[ ion-3 ] += pr2;
			pr_tmp[ ion+4 ] += pr;
			pr_tmp[ ion-4 ] += pr;
		}
	}
	
	// screen out decarboxy form of precursor ion
	ion = lrint( f->mass + 10.0 * ( 3*MASS_H -MASS_C - MASS_O ) - 0.5 );
	pr_tmp[ ion ] = 0;
	pr_tmp[ ion+1 ] = 0;
	pr_tmp[ ion-1 ] = 0;
	pr_tmp[ ion+2 ] = 0;
	pr_tmp[ ion-2 ] = 0;
	pr_tmp[ ion+3 ] = 0;
	pr_tmp[ ion-3 ] = 0;
	pr_tmp[ ion+4 ] = 0;

	delta = 2.0 * SCALE_REGION / (double) mass_int;
	r = 0;
	sum_of_pr = 0;
	for( i = 0; i < mass_int; i++ )
	{
		int	num;
	
		if( ! (pr_tmp[i] >0 ) ){
			r += delta;
			continue;
		}

		num = lrint( r );

		pr_tmp[i] *= scale_factor[ num ];
		sum_of_pr += pr_tmp[i];

		r += delta;
	}

        sum_of_noise = ( noise_1/(100.0-noise_1) )*sum_of_pr;
        average = sum_of_noise / mass_int;

	for( i = 0; i < mass_int; i++ )
            	pr_tmp[i] += average;

	f->num_of_ions1 = mass_int;
	f->ideal_spectrum_charge1 = pr_tmp;
	f->sum_of_pr1 = sum_of_pr;

    	return;
}

int	produce_ideal_spectrum_charge2( struct _FRAGMENT_ *f, char *fragment, int count )
{
	int	i, j, index;
	double	mass;
	double	score;
	long	ion, mass_int;
        double  delta_b, delta_y;
	double	pr, pr2, pr3, pr4, pr5;
	double	*pr_tmp;
        double  sum_of_pr;
        double  sum_of_noise;
        double  average;
	double	middle_point;
	double	r, delta;

        delta_b  = 10.0 * MASS_H;
        delta_y  = 10.0 * ( 3*MASS_H + MASS_O );

	mass_int = lrint( f->mass + 100 );
	middle_point = mass_int / 2.0;

	pr_tmp = malloc( mass_int * sizeof(double) );
	if( pr_tmp == NULL ){
		printf("Malloc error!\n");
		exit(-1);
	}

	for( i = 0; i < mass_int; i++ )
		pr_tmp[i] = 0;

	mass = 0.0;
	for( i = 0; i < count; i++ ){

		if( AA_Table[ fragment[i] ] < 0 ){
			printf("Unknow AA id: %c. Check it!\n", fragment[i] );
			exit(-1);
		}

		mass += AA_Table[ fragment[i] ];
		for( j = 0; j < num_of_sN2; j++ ){
			ion = lrint( mass + delta_b + Pr_sN2[j].offset - 0.5 );
			if( ion <= 0 || ion >= mass_int )
				continue;
			pr = Pr_sN2[j].pr / 25.0;
			pr2 = pr + pr;
			pr3 = pr2 + pr;
			pr4 = pr3 + pr;
			pr5 = pr4 + pr;

			pr_tmp[ ion ] += pr5;
			pr_tmp[ ion+1 ] += pr4;
			pr_tmp[ ion-1 ] += pr4;
			pr_tmp[ ion+2 ] += pr3;
			pr_tmp[ ion-2 ] += pr3;
			pr_tmp[ ion+3 ] += pr2;
			pr_tmp[ ion-3 ] += pr2;
			pr_tmp[ ion+4 ] += pr;
			pr_tmp[ ion-4 ] += pr;
		}

		// deal charge 2
		for( j = 0; j < num_of_dN2; j++ ){
			ion = lrint( ( mass + delta_b )/2.0 + Pr_dN2[j].offset - 0.5 );
			if( ion <= 0 || ion >= mass_int )
				continue;
			pr = Pr_dN2[j].pr / 25.0;
			pr2 = pr + pr;
			pr3 = pr2 + pr;
			pr4 = pr3 + pr;
			pr5 = pr4 + pr;

			pr_tmp[ ion ] += pr5;
			pr_tmp[ ion+1 ] += pr4;
			pr_tmp[ ion-1 ] += pr4;
			pr_tmp[ ion+2 ] += pr3;
			pr_tmp[ ion-2 ] += pr3;
			pr_tmp[ ion+3 ] += pr2;
			pr_tmp[ ion-3 ] += pr2;
			pr_tmp[ ion+4 ] += pr;
			pr_tmp[ ion-4 ] += pr;
		}
	}

	mass = 0.0;
	for( i = count-1; i > 0; i-- ){
		if( AA_Table[ fragment[i] ] < 0 ){
			printf("Unknow AA id: %c. Check it!\n", fragment[i] );
			exit(-1);
		}
		mass += AA_Table[ fragment[i] ];
                
		for( j = 0; j < num_of_sC2; j++ ){
			ion = lrint( mass + delta_y + Pr_sC2[j].offset - 0.5 );
			if( ion <= 0 || ion >= mass_int )
				continue;
//			assert( ion < p->num_of_peaks );
//			assert( ion > 4 );
			pr = Pr_sC2[j].pr / 25.0;
			pr2 = pr + pr;
			pr3 = pr2 + pr;
			pr4 = pr3 + pr;
			pr5 = pr4 + pr;

			pr_tmp[ ion ] += pr5;
			pr_tmp[ ion+1 ] += pr4;
			pr_tmp[ ion-1 ] += pr4;
			pr_tmp[ ion+2 ] += pr3;
			pr_tmp[ ion-2 ] += pr3;
			pr_tmp[ ion+3 ] += pr2;
			pr_tmp[ ion-3 ] += pr2;
			pr_tmp[ ion+4 ] += pr;
			pr_tmp[ ion-4 ] += pr;
		}

                // deal charge 2
		for( j = 0; j < num_of_dC2; j++ ){
			ion = lrint( ( mass + delta_y )/2.0 + Pr_dC2[j].offset - 0.5 );
			if( ion <= 0 || ion >= mass_int )
				continue;
			pr = Pr_dC2[j].pr / 25.0;
			pr2 = pr + pr;
			pr3 = pr2 + pr;
			pr4 = pr3 + pr;
			pr5 = pr4 + pr;

			pr_tmp[ ion ] += pr5;
			pr_tmp[ ion+1 ] += pr4;
			pr_tmp[ ion-1 ] += pr4;
			pr_tmp[ ion+2 ] += pr3;
			pr_tmp[ ion-2 ] += pr3;
			pr_tmp[ ion+3 ] += pr2;
			pr_tmp[ ion-3 ] += pr2;
			pr_tmp[ ion+4 ] += pr;
			pr_tmp[ ion-4 ] += pr;
		}
	}

        // screen out precursor ion
        ion = lrint( ( f->mass + 10.0 * ( 4*MASS_H + MASS_O ) )/2.0 - 0.5 );
        
	pr_tmp[ ion ] = 0;
	pr_tmp[ ion+1 ] = 0;
	pr_tmp[ ion-1 ] = 0;
	pr_tmp[ ion+2 ] = 0;
	pr_tmp[ ion-2 ] = 0;
	pr_tmp[ ion+3 ] = 0;
	pr_tmp[ ion-3 ] = 0;
	pr_tmp[ ion+4 ] = 0;
	pr_tmp[ ion-4 ] = 0;

       // screen out dehydrated form of precursor ion
	ion = lrint( ( f->mass + 10.0 * ( 2*MASS_H ) )/2.0 - 0.5 );
	pr_tmp[ ion ] = 0;
	pr_tmp[ ion+1 ] = 0;
	pr_tmp[ ion-1 ] = 0;
	pr_tmp[ ion+2 ] = 0;
	pr_tmp[ ion-2 ] = 0;
	pr_tmp[ ion+3 ] = 0;
	pr_tmp[ ion-3 ] = 0;
	pr_tmp[ ion+4 ] = 0;
       
	// screen out decarboxy form of precursor ion
	ion = lrint( ( f->mass + 10.0 * ( 3*MASS_H -MASS_C - MASS_O ) )/2.0 - 0.5 );
	pr_tmp[ ion ] = 0;
	pr_tmp[ ion+1 ] = 0;
	pr_tmp[ ion-1 ] = 0;
	pr_tmp[ ion+2 ] = 0;
	pr_tmp[ ion-2 ] = 0;
	pr_tmp[ ion+3 ] = 0;
	pr_tmp[ ion-3 ] = 0;
	pr_tmp[ ion+4 ] = 0;
	pr_tmp[ ion-4 ] = 0;

	delta = 2.0 * SCALE_REGION / (double) mass_int;
	r = 0;
	sum_of_pr = 0;
	for( i = 0; i < mass_int; i++ )
	{
		int	num;
	
		if( ! (pr_tmp[i] >0 ) ){
			r += delta;
			continue;
		}

		num = lrint( r );

		pr_tmp[i] *= scale_factor[ num ];
		sum_of_pr += pr_tmp[i];

		r += delta;
	}
        sum_of_noise = ( noise_2/(100.0-noise_2) )*sum_of_pr;
        average = sum_of_noise / mass_int;

	for( i = 0; i < mass_int; i++ )
            	pr_tmp[i] += average;
        
        f->num_of_ions2 = mass_int;
	f->ideal_spectrum_charge2 = pr_tmp;
	f->sum_of_pr2 = sum_of_pr;
        
        return;
}

double	compare_ideal_real_spectrum( _DTA_ *p, _FRAGMENT_ *f, long int *ion_count )
{
	double	score;
        long int        min_ion, max_ion;
	long int	count;
	long int	num_of_ions;
	double	*ideal_spectrum;
	double	pi, qi;
	double	sum_of_pr;
	long int        i;

#ifdef	PRINT_REAL_IDEAL
fprintf( stderr, "===== %s\n", p->ID );
#endif
	
	if( p == NULL )
		return	-1;

	if( p->charge == 1 && 
	    fabs( f->mass + 10 * (2*MASS_H + MASS_O + MASS_H) - p->precursor  ) > ERROR_OF_PRECURSOR_CHARGE1 ){
		return	-1;
	}
	if( p->charge == 2 && 
	    fabs( f->mass + 10 * ( 2*MASS_H + MASS_O + 2*MASS_H ) - p->precursor ) > ERROR_OF_PRECURSOR_CHARGE2 ) {
		return	-1;
	}

	// work out each points pr
	if( p->charge == 1 ){
		num_of_ions = f->num_of_ions1;
		ideal_spectrum = f->ideal_spectrum_charge1;
	}
	if( p->charge == 2 ){
		num_of_ions = f->num_of_ions2;
		ideal_spectrum = f->ideal_spectrum_charge2;
	}

        min_ion = p->min_ion;
        if( p->max_ion > num_of_ions )
                max_ion = num_of_ions;
        else
                max_ion = p->max_ion;

	sum_of_pr = 0;
	for( i = min_ion; i <= max_ion; i++ )
		sum_of_pr += ideal_spectrum[i];

	score = 0.0;
	for( i = min_ion; i < max_ion; i++ ){

		count++;
//		if( i >= p->num_of_peaks )
//			break;
	
		pi = p->real_spectrum[i] ;
		qi = ideal_spectrum[i] / sum_of_pr;

		if( pi > 0 && qi > 0 ) {
			double	delta;
			delta = ((log(pi) - log(qi)) * ( pi - qi ) );
			score += delta;
#ifdef	PRINT_REAL_IDEAL
	fprintf( stderr, "%.1lf\t%.10lf\t%.10lf\t%.10lf\t%d\t%d\t%d\n", (double) i / 10.0, pi, qi, delta, min_ion, max_ion, p->num_of_peaks );
#endif
		}
	}

	*ion_count = count;
	if( count == 0 )
		return	-1;
    
    	return	score;
}

void	compute_score_for_dta_list()
{
	struct	_FRAGMENT_ *f;
    	long int     i, j, ion_count;
    	double  score, best_score;
        long    dta_index;
        int     count;

	for (i = 0; i < total_num_of_fragments; i++) {

                count = 0;
        	f = fragment_library[i];
        	produce_ideal_spectrum_charge1( f, f->sequence, f->length );
	    	produce_ideal_spectrum_charge2( f, f->sequence, f->length );
               
 		for( j = 0; j < num_of_dta; j++ ) {
    	    		if( dta_list[j] == NULL )
	            		continue;
                        
                        best_score = -1;
                        
 			score = compare_ideal_real_spectrum( dta_list[j], f, &ion_count );

            		if( score < 0.0 
// || score > THRESHOLD
			)
                		continue;

            		if( dta_list[j]->best_score < 0 || dta_list[j]->best_score > score )
                        {
                		dta_list[j]->best_score = score;
                		dta_list[j]->index_of_best_matched_fragment = i;
                                // initial value of num_of_matched_fragments is 0 but not -1;
                                dta_list[j]->num_of_matched_fragments ++;
	            	}

                        if( best_score < 0 || best_score > score )
                        {
                                best_score = score;
                                dta_index = j;
                                count ++;
                        }
        	}

                f->best_score = best_score;
                f->index_of_best_matched_dta = dta_index;
                f->num_of_matched_dtas = count;                
                
        	if (f->ideal_spectrum_charge1 != NULL) {
            		free( f->ideal_spectrum_charge1 );
            		f->ideal_spectrum_charge1 = NULL;
        	}
   		if (f->ideal_spectrum_charge2 != NULL) {
            		free( f->ideal_spectrum_charge2 );
            		f->ideal_spectrum_charge2 = NULL;
        	}
    	}
     
    	return;
}
                 		                                    		
void process_dta_list()
{
    	long int    i, j;
    	long int count;

    	total_num_of_fragments = 0;
    	for (i = 0; i < num_of_proteins; i++) 
        	total_num_of_fragments += protein_list[i]->num_of_fragments;
   
    	fragment_library = malloc( total_num_of_fragments * sizeof(struct _FRAGMENT_ *) );

    	if( fragment_library == NULL ) {
        	printf("Memory request failed! For the fragment library\n");
        	exit(-1);
    	}
    
    	count = 0;
    	for (i=0; i<num_of_proteins; i++) {
        	for (j=0; j<protein_list[i]->num_of_fragments; j++ ) {
            		fragment_library[count] = protein_list[i]->fragment_list[j];
            		count ++;
        	} 
    	}

    	compute_score_for_dta_list();

	return;
}

void    output( char *result_filename )
{
        FILE    *fp;
        long int    index;
        long int i;
    
	fp = fopen( result_filename, "a" );
	if( fp == NULL ){
		printf("Can't open file %s!\n", result_filename );
		exit(-1);
	}

        for( i=0; i<num_of_dta; i++ )
        {
       	        if( dta_list[i] == NULL )
       	                continue;            
//              if( dta_list[i]->score < 0 || dta_list[i]->score > THRESHOLD )
                if( dta_list[i]->best_score < 0 )
                        continue;
                index = dta_list[i]->index_of_best_matched_fragment;    
                fprintf( fp, "%s\t%s\t%5.2lf\n", fragment_library[index]->sequence, dta_list[i]->ID, dta_list[i]->best_score ); 
        }
        fclose( fp );  
        
        return;
}

void free_memory()
{
        long int    i;

        free( digest_site );

        free( Pr_N1 );
        free( Pr_C1 );
        free( Pr_sN2 );
        free( Pr_dN2 );
        free( Pr_sC2 );
        free( Pr_dC2 );
        
	for (i=0; i<num_of_dta; i++ ){
		if( dta_list[i] == NULL )
			continue;
            
		if( dta_list[i]->real_spectrum != NULL )
                {
                        free( dta_list[i]->real_spectrum );
                        dta_list[i]->real_spectrum = NULL;
                }
		if( dta_list[i]->ID != NULL )
                {
                        free( dta_list[i]->ID );
                        dta_list[i]->ID = NULL;
                }
                free( dta_list[i] );
                dta_list[i] = NULL;
        }
	free( dta_list );    
	dta_list = NULL;

//content of "protein_list[i]->fragment_list[j]" will be deleted in the next step.    
	for (i=0; i<num_of_proteins; i++)
        {
                if (protein_list[i]->fragment_list != NULL)
                {
                        free( protein_list[i]->fragment_list );	
                        protein_list[i]->fragment_list = NULL;
                }
                if (protein_list[i]->protein_sequence != NULL)
                {
                        free( protein_list[i]->protein_sequence );	
                        protein_list[i]->protein_sequence = NULL;
                }
                if (protein_list[i]->protein_name != NULL)
                {
                        free( protein_list[i]->protein_name );	
                        protein_list[i]->protein_name = NULL;
                }
        }
        free( protein_list );
        protein_list = NULL;
	
        for (i=0; i<total_num_of_fragments; i++)
        {
                if (fragment_library[i]->sequence != NULL)
                {
                        free( fragment_library[i]->sequence );
                        fragment_library[i]->sequence = NULL;
                }
            
                if (fragment_library[i]->ideal_spectrum_charge1 != NULL)
                {
                        free( fragment_library[i]->ideal_spectrum_charge1 );
                        fragment_library[i]->ideal_spectrum_charge1 = NULL;
                }
            
                if (fragment_library[i]->ideal_spectrum_charge2 != NULL)
                {
                        free( fragment_library[i]->ideal_spectrum_charge2 );
                        fragment_library[i]->ideal_spectrum_charge2 = NULL;
                }
            
                free( fragment_library[i] );
                fragment_library[i] = NULL;
        }
        free( fragment_library );
        fragment_library = NULL;
    
        return;
}

void	init_scale_factor()
{
	int	i;
	double	e;

	for( i = 0; i <= 2*SCALE_REGION ; i++ ) {
		e = ( i - SCALE_REGION ) / (double) SCALE_REGION / 2.0;
		
		e = e * e;
		
		scale_factor[i] = exp( - e * (double)SIGMA );
		//printf("bdb	%d	%.10lf\n", i, scale_factor[i] );
	}
}

int main( int argc, char **argv )
{

	if( argc < 5 ){
		printf("Usage: %s dta.list protein.list enzyme.cfg probability.cfg match.out\n", argv[0]);
		exit(-1);
	}

	init_AA_table();

	init_scale_factor();

	readin_enzyme_site_file( argv[3] );

	if( initiate_probability_table( argv[4] ) < 0 )
        {
                printf( "It fails to read in %s!\n", argv[4] );
		exit(-1);
	}
    
	if( readin_dta_list( argv[1] ) < 0 ){
		printf("Read dta list file %s error!\n", argv[1]);
		exit(-1);
	}
        
        if (read_and_digest_proteins( argv[2] ) < 0) {
		printf("Read  and digest protein database %s error!\n", argv[2]);
		exit(-1);
	}
    
        process_dta_list();

        output( argv[5] );
    
    	free_memory();

	return 1;
}