/********* Sequence input routines for CLUSTAL W *******************/
/* DES was here.  FEB. 1994 */
/* Now reads PILEUP/MSF and CLUSTAL alignment files */

#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include "clustalw.h"	

#define MIN(a,b) ((a)<(b)?(a):(b))



/*
*	Prototypes
*/
static sint get_seq_index(SEQ *seqs,sint nseqs,sint ii);

/*
 *	Global variables
 */
static sint debug;


Boolean open_alignment_output(char *filename,ALNOUT_OPTPTR alnout_opt)
{

	if(!alnout_opt->output_clustal && !alnout_opt->output_nbrf && !alnout_opt->output_gcg &&
		 !alnout_opt->output_phylip && !alnout_opt->output_gde) {
                error("You must select an alignment output format");
                return FALSE;
        }

	if(alnout_opt->output_clustal) 
		if((alnout_opt->clustal_outfile=open_output_file(filename, "CLUSTAL", "aln", alnout_opt->clustal_outname))==NULL) return FALSE;

	if(alnout_opt->output_nbrf) 
		if((alnout_opt->nbrf_outfile=open_output_file(filename, "NBRF", "pir", alnout_opt->nbrf_outname))==NULL) return FALSE;

	if(alnout_opt->output_gcg) 
		if((alnout_opt->gcg_outfile=open_output_file(filename, "GCG", "msf", alnout_opt->gcg_outname))==NULL) return FALSE;

	if(alnout_opt->output_phylip) 
		if((alnout_opt->phylip_outfile=open_output_file(filename, "PHYLIP", "phy", alnout_opt->phylip_outname))==NULL) return FALSE;

	if(alnout_opt->output_gde) 
		if((alnout_opt->gde_outfile=open_output_file(filename, "GDE", "gde", alnout_opt->gde_outname))==NULL) return FALSE;

	return TRUE;
}

/* gets a filename from the user and opens the file 
returns the file handle (the filename is written into char *file_name) */

FILE *  open_output_file(char *in_name, char *prompt, char *file_extension, char *out_name)
{
	char temp[FILENAMELEN+1];
	char path[FILENAMELEN+1];
	char local_prompt[MAXLINE];
	FILE * file_handle;
	Boolean usemenu;

/* if the output filename is already specified, just open the file and return */
	if (out_name[0]!=EOS) {
		file_handle = open_explicit_file(out_name);
		return file_handle;
	}

	usemenu=get_usemenu();

/* otherwise, suggest a name to the user and prompt if this is ok */
	get_path(in_name,path);
	strcpy(out_name,path);
	strcat(out_name,file_extension);

	if(strcmp(out_name,in_name)==0) {
		warning("Output file name is the same as input file.");
		if (usemenu) {
			strcpy(local_prompt,"\n\nEnter new name to avoid overwriting ");
			strcat(local_prompt," [%s]: ");          
			fprintf(stdout,local_prompt,out_name);
			gets(temp);
			if(*temp != EOS) strcpy(out_name,temp);
		}
	}
	else if (usemenu) {
		strcpy(local_prompt,"\nEnter a name for the ");
		strcat(local_prompt,prompt);
		strcat(local_prompt," output file ");
		strcat(local_prompt," [%s]: ");          
		fprintf(stdout,local_prompt,out_name);
		gets(temp);
		if(*temp != EOS) strcpy(out_name,temp);
	}

#ifdef VMS
	if((file_handle=fopen(out_name,"w","rat=cr","rfm=var"))==NULL) {
#else
	if((file_handle=fopen(out_name,"w"))==NULL) {
#endif
		error("Cannot open output file [%s]",out_name);
		return NULL;
	}
	return file_handle;
}



FILE *  open_explicit_file(char *file_name)
{ 
	FILE * file_handle;

	if (*file_name == EOS) {
		error("Bad output file [%s]",file_name);
		return NULL;
	}
#ifdef VMS
	if((file_handle=fopen(file_name,"w","rat=cr","rfm=var"))==NULL) {
#else
	if((file_handle=fopen(file_name,"w"))==NULL) {
#endif
		error("Cannot open output file [%s]",file_name);
		return NULL;
	}
	return file_handle;
}

void create_alignment_output(ALN mult_aln,sint fseq, sint nseq,ALNOUT_OPT alnout_opt,SS_OPT ss_opt)
{
	sint i,length;
	Boolean usemenu;

	length=0;
        for (i=fseq;i<nseq;i++)
                if (length < mult_aln.seqs[i].len)
                        length = mult_aln.seqs[i].len;

	usemenu=get_usemenu();

	if (usemenu) info("Consensus length = %d",(pint)length);
	if(alnout_opt.output_clustal) {
		clustal_out(alnout_opt.clustal_outfile, 0, length, fseq, nseq,alnout_opt.seq_numbers,mult_aln,ss_opt);
		fclose(alnout_opt.clustal_outfile);
		info("CLUSTAL-Alignment file created  [%s]",alnout_opt.clustal_outname);
	}
	if(alnout_opt.output_nbrf)  {
		nbrf_out(alnout_opt.nbrf_outfile, 0, length, fseq, nseq,mult_aln);
		fclose(alnout_opt.nbrf_outfile);
		info("NBRF/PIR-Alignment file created [%s]",alnout_opt.nbrf_outname);
	}
	if(alnout_opt.output_gcg)  {
		gcg_out(alnout_opt.gcg_outfile, 0, length, fseq, nseq,mult_aln);
		fclose(alnout_opt.gcg_outfile);
		info("GCG-Alignment file created      [%s]",alnout_opt.gcg_outname);
	}
	if(alnout_opt.output_phylip)  {
		phylip_out(alnout_opt.phylip_outfile, 0, length, fseq, nseq,mult_aln);
		fclose(alnout_opt.phylip_outfile);
		info("PHYLIP-Alignment file created   [%s]",alnout_opt.phylip_outname);
	}
	if(alnout_opt.output_gde)  {
		gde_out(alnout_opt.gde_outfile, 0, length, fseq, nseq,alnout_opt.lowercase,mult_aln,ss_opt);
		fclose(alnout_opt.gde_outfile);
		info("GDE-Alignment file created      [%s]",alnout_opt.gde_outname);
	}
}


void clustal_out(FILE *clusout, sint fres, sint len, sint fseq, sint nseq,Boolean seq_numbers,ALN mult_aln,SS_OPT ss_opt)
{
	static char *seq1;
	static sint *seq_no;
	static sint *print_seq_no;
	char *revision_level;
	char *ss_mask1, *ss_mask2;
	char 	    temp[MAXLINE];
	char c;
	char val;
	sint ii,lv1,catident1[NUMRES],catident2[NUMRES],ident,chunks;
	sint i,j,k,l;
	sint pos,ptr;
	sint line_length;
	sint max_aln_length,max_names;

/*These are all the positively scoring groups that occur in the Gonnet Pam250
matrix. There are strong and weak groups, defined as strong score >0.5 and
weak score =<0.5. Strong matching columns to be assigned ':' and weak matches
assigned '.' in the clustal output format.
*/

char *res_cat1[] = {
                "sta",
                "neqk",
                "nhqk",
                "ndeq",
                "qhrk",
                "milv",
                "milf",
                "hy",
                "fyw",
                NULL };

char *res_cat2[] = {
                "csa",
                "atv",
                "sag",
                "stnk",
                "stpa",
                "sgnd",
                "sndeqk",
                "ndeqhk",
                "neqhrk",
                "fvlim",
                "hfy",
                NULL };


	seq_no = (sint *)ckalloc((nseq+1) * sizeof(sint));
	print_seq_no = (sint *)ckalloc((nseq+1) * sizeof(sint));
	max_aln_length=max_names=0;
	for (i=fseq;i<nseq;i++)
        {
		if(mult_aln.seqs[i].len>max_aln_length) max_aln_length=mult_aln.seqs[i].len;
		if(strlen(mult_aln.seqs[i].name)>max_names) max_names=strlen(mult_aln.seqs[i].name);
		print_seq_no[i] = seq_no[i] = 0;
                for(j=0;j<fres;j++) {
			val = mult_aln.seqs[i].data[j];
			if(isalpha(val)) seq_no[i]++;
		}
        }


	seq1 = (char *)ckalloc((max_aln_length+1) * sizeof(char));

	if (mult_aln.prf1.ss.struct_penalties == SECST && ss_opt.use_ss1 == TRUE) {
		ss_mask1=calc_sec_struct_mask(mult_aln.seqs[0].len,mult_aln.prf1.ss.sec_struct_mask,ss_opt);
	}
	if (mult_aln.prf2.ss.struct_penalties == SECST && ss_opt.use_ss2 == TRUE) {
		ss_mask2=calc_sec_struct_mask(mult_aln.seqs[mult_aln.prf1.nseqs].len,mult_aln.prf2.ss.sec_struct_mask,ss_opt);
	}
	
	revision_level=get_revision_level();
	fprintf(clusout,"CLUSTAL %s multiple sequence alignment\n\n",
				 revision_level);

/* decide the line length for this alignment - maximum is LINELENGTH */
	line_length=PAGEWIDTH-max_names;
	line_length=line_length-line_length % 10; /* round to a multiple of 10*/
	if (line_length > LINELENGTH) line_length=LINELENGTH;

	chunks = len/line_length;
	if(len % line_length != 0)
		++chunks;
		
	for(lv1=0;lv1<chunks;lv1++) {
		pos = lv1*line_length;
		ptr = (len<pos+line_length) ? len : pos+line_length;

		fprintf(clusout,"\n");
		
		if (ss_opt.output_struct_penalties == 0 || ss_opt.output_struct_penalties == 2) {
			if (mult_aln.prf1.ss.struct_penalties == SECST && ss_opt.use_ss1 == TRUE) {
				for(i=pos;i<ptr;i++) {
					val=ss_mask1[i+fres];
					if (!isalpha(val))
						temp[i-pos]='-';
					else
						temp[i-pos]=val;
				}
				temp[ptr-pos]=EOS;
				fprintf(clusout,"!SS_%-*s  %s\n",max_names,mult_aln.prf1.ss.name,temp);
				}
		}
		if (ss_opt.output_struct_penalties == 1 || ss_opt.output_struct_penalties == 2) {
			if (mult_aln.prf1.ss.struct_penalties != NONE && ss_opt.use_ss1 == TRUE) {
				for(i=pos;i<ptr;i++) {
					val=mult_aln.prf1.ss.gap_penalty_mask[i+fres-2];
					if (!isalpha(val))
						temp[i-pos]='-';
					else
						temp[i-pos]=val;
				}
				temp[ptr-pos]=EOS;
				fprintf(clusout,"!GM_%-*s  %s\n",max_names,mult_aln.prf1.ss.name,temp);
				}
		}
		if (ss_opt.output_struct_penalties == 0 || ss_opt.output_struct_penalties == 2) {
			if (mult_aln.prf2.ss.struct_penalties == SECST && ss_opt.use_ss2 == TRUE) {
				for(i=pos;i<ptr;i++) {
					val=ss_mask2[i+fres];
					if (!isalpha(val))
						temp[i-pos]='-';
					else
						temp[i-pos]=val;
				}
				temp[ptr-pos]=EOS;
				fprintf(clusout,"!SS_%-*s  %s\n",max_names,mult_aln.prf2.ss.name,temp);
			}
		}
		if (ss_opt.output_struct_penalties == 1 || ss_opt.output_struct_penalties == 2) {
			if (mult_aln.prf2.ss.struct_penalties != NONE && ss_opt.use_ss2 == TRUE) {
				for(i=pos;i<ptr;i++) {
					val=mult_aln.prf2.ss.gap_penalty_mask[i+fres];
					if (!isalpha(val))
						temp[i-pos]='-';
					else
						temp[i-pos]=val;
				}
				temp[ptr-pos]=EOS;
				fprintf(clusout,"!GM_%-*s  %s\n",max_names,mult_aln.prf2.ss.name,temp);
			}
		}
		
		for(ii=fseq;ii<nseq;ii++) {
			i=get_seq_index(mult_aln.seqs,nseq,ii);
			print_seq_no[i] = 0;
			for(j=pos;j<ptr;j++) {
				if (j+fres<mult_aln.seqs[i].len)
					val = mult_aln.seqs[i].data[j+fres];
				else val = -3;
				if(val == EOS) break;
				else if(!isalpha(val))
                                        seq1[j]='-';
				else {
					seq1[j]=toupper(val);
					seq_no[i]++;
					print_seq_no[i]=1;
				}
			}
			for(;j<ptr;j++) seq1[j]='-';
			strncpy(temp,&seq1[pos],ptr-pos);
			temp[ptr-pos]=EOS;
			fprintf(clusout,"%-*s %s",max_names+5,mult_aln.seqs[i].name,temp);
			if (seq_numbers && print_seq_no[i])
					fprintf(clusout," %d",seq_no[i]);
			fprintf(clusout,"\n");
		}
	
		for(i=pos;i<ptr;i++) {
			seq1[i]=' ';
			ident=0;
			for(j=0;res_cat1[j]!=NULL;j++) catident1[j] = 0;
			for(j=0;res_cat2[j]!=NULL;j++) catident2[j] = 0;
			for(j=fseq;j<nseq;j++) {
				if(isalpha(mult_aln.seqs[fseq].data[i])) {
					if(mult_aln.seqs[fseq].data[i] == mult_aln.seqs[j].data[i])
					++ident;
					for(k=0;res_cat1[k]!=NULL;k++) {
					        for(l=0;(c=res_cat1[k][l]);l++) {
							if (mult_aln.seqs[j].data[i]==c)
							{
								catident1[k]++;
								break;
							}
						}
					}
					for(k=0;res_cat2[k]!=NULL;k++) {
					        for(l=0;(c=res_cat2[k][l]);l++) {
							if (mult_aln.seqs[j].data[i]==c)
							{
								catident2[k]++;
								break;
							}
						}
					}
				}
			}
			if(ident==nseq-fseq)
				seq1[i]='*';
			else if (!mult_aln.dnaflag) {
				for(k=0;res_cat1[k]!=NULL;k++) {
					if (catident1[k]==nseq-fseq) {
						seq1[i]=':';
						break;
					}
				}
				if(seq1[i]==' ')
					for(k=0;res_cat2[k]!=NULL;k++) {
						if (catident2[k]==nseq-fseq) {
							seq1[i]='.';
							break;
						}
					}
			}
		}
		strncpy(temp,&seq1[pos],ptr-pos);
		temp[ptr-pos]=EOS;
		for(k=0;k<max_names+6;k++) fprintf(clusout," ");
		fprintf(clusout,"%s\n",temp);
	}
		

	seq1=ckfree((void *)seq1);
	if (mult_aln.prf1.ss.struct_penalties == SECST && ss_opt.use_ss1 == TRUE) ckfree(ss_mask1);
	if (mult_aln.prf2.ss.struct_penalties == SECST && ss_opt.use_ss2 == TRUE) ckfree(ss_mask2);
	
}


static sint get_seq_index(SEQ *seqs,sint nseqs,sint ii)
{
	sint i;

	for(i=0;i<nseqs;i++) {
		if(seqs[i].output_index==ii) break;
	}

	return i;
}



void gcg_out(FILE *gcgout, sint fres, sint len, sint fseq, sint nseq,ALN mult_aln)
{
	char *seq, residue;
	char val;
	sint *all_checks;
	sint i,ii,chunks,block;
	sint j,pos1,pos2;	
	sint max_aln_length,max_names;
	long grand_checksum;

	max_aln_length=max_names=0;
	for(i=fseq; i<nseq; i++) {
		if(mult_aln.seqs[i].len>max_aln_length) max_aln_length=mult_aln.seqs[i].len;
		if(strlen(mult_aln.seqs[i].name)>max_names) max_names=strlen(mult_aln.seqs[i].name);
	}

	seq = (char *)ckalloc((max_aln_length+1) * sizeof(char));
	all_checks = (sint *)ckalloc((nseq+1) * sizeof(sint));
	
	for(i=fseq; i<nseq; i++) {
		for(j=fres; j<fres+len; j++) {
			val = mult_aln.seqs[i].data[j];
			if(val == EOS) break;
			else if(!isalpha(val))
				residue = '.';
                        else {
				residue = val;
			}
			seq[j-fres] = toupper(residue);
		}
/* pad any short sequences with gaps, to make all sequences the same length */
		for(; j<fres+len; j++) 
			seq[j-fres] = '.';
		all_checks[i] = SeqGCGCheckSum(seq, (int)len);
	}	

	grand_checksum = 0;
	for(i=0; i<nseq; i++) grand_checksum += all_checks[i];
	grand_checksum = grand_checksum % 10000;
        fprintf(gcgout,"PileUp\n\n");
	fprintf(gcgout,"\n\n   MSF:%5d  Type: ",(pint)len);
	if(mult_aln.dnaflag)
		fprintf(gcgout,"N");
	else
		fprintf(gcgout,"P");
	fprintf(gcgout,"    Check:%6ld   .. \n\n", (long)grand_checksum);
	for(ii=fseq; ii<nseq; ii++)  {
		i = get_seq_index(mult_aln.seqs,nseq,ii);
		fprintf(gcgout,
			" Name: %s oo  Len:%5d  Check:%6ld  Weight:  %.1f\n",
			mult_aln.seqs[i].name,(pint)len,(long)all_checks[i],(float)mult_aln.seqs[i].weight*100.0/(float)INT_SCALE_FACTOR);
        }
	fprintf(gcgout,"\n//\n");  

	chunks = len/GCG_LINELENGTH;
	if(len % GCG_LINELENGTH != 0) ++chunks;

	for(block=0; block<chunks; block++) {
		fprintf(gcgout,"\n\n");
		pos1 = (block * GCG_LINELENGTH);
		pos2 = (len<pos1+GCG_LINELENGTH)? len : pos1+GCG_LINELENGTH;
		for(ii=fseq; ii<nseq; ii++) {
			i = get_seq_index(mult_aln.seqs,nseq,ii);
			fprintf(gcgout,"\n%-*s ",max_names+5,mult_aln.seqs[i].name);
			for(j=pos1; j<pos2; j++) {
/*
    JULIE -
    check for sint sequences - pad out with '.' characters to end of alignment
*/
				if (j+fres<mult_aln.seqs[i].len)
					val = mult_aln.seqs[i].data[j+fres];
				else val = -3;
				if(val == EOS)
					residue = '.';
				else if(!isalpha(val))
					residue = '.';
				else {
					residue = toupper(val);
				}
				fprintf(gcgout,"%c",residue);
				if((j+1) % 10 == 0) fprintf(gcgout," ");
			}
		}
	}

	seq=ckfree((void *)seq);
	all_checks=ckfree((void *)all_checks);
	
	fprintf(gcgout,"\n\n");
}



void phylip_out(FILE *phyout, sint fres, sint len, sint fseq, sint nseq,ALN mult_aln)
{
	char residue;
	char val;
	sint i,ii,chunks,block;	
	sint j,pos1,pos2;	
	sint name_len;
	Boolean warn;
	char **snames;
	
	snames=(char **)ckalloc((nseq-fseq+2)*sizeof(char *));
	name_len=0;
	for(i=fseq; i<nseq; i++)  {
		snames[i]=(char *)ckalloc((11)*sizeof(char));
		ii=strlen(mult_aln.seqs[i].name);
		strncpy(snames[i],mult_aln.seqs[i].name,10);
		if(name_len<ii) name_len=ii;
	}
	if(name_len>10) {
		warn=FALSE;
		for(i=fseq; i<nseq; i++)  {
                	for(j=i+1;j<nseq;j++) {
                        	if (strcmp(snames[i],snames[j]) == 0) 
					warn=TRUE;
                	}
        	}
		if(warn)
			warning("Truncating sequence names to 10 characters for PHYLIP output.\n"
			"Names in the PHYLIP format file are NOT unambiguous.");
		else
			warning("Truncating sequence names to 10 characters for PHYLIP output.");
	}


	chunks = len/GCG_LINELENGTH;
	if(len % GCG_LINELENGTH != 0) ++chunks;

	fprintf(phyout,"%6d %6d",(pint)nseq-fseq,(pint)len);

	for(block=0; block<chunks; block++) {
		pos1 = (block * GCG_LINELENGTH);
		pos2 = (len<pos1+GCG_LINELENGTH)? len : pos1+GCG_LINELENGTH;
		for(ii=fseq; ii<nseq; ii++)  {
			i = get_seq_index(mult_aln.seqs,nseq,ii);
			if(block == 0)  {
				fprintf(phyout,"\n%-10s ",snames[i]);
			}
			else
				fprintf(phyout,"\n           ");
			for(j=pos1; j<pos2; j++) {
				if (j+fres<mult_aln.seqs[i].len)
					val = mult_aln.seqs[i].data[j+fres];
				else val = -3;
				if(val == EOS)
					break;
				else if(!isalpha(val))
					residue = '-';
				else {
					residue = toupper(val);
				}
				fprintf(phyout,"%c",residue);
				if((j+1) % 10 == 0) fprintf(phyout," ");
			}
		}
		fprintf(phyout,"\n");
	}

	for(i=fseq;i<nseq;i++)
		ckfree(snames[i]);
	ckfree(snames);
}





void nbrf_out(FILE *nbout, sint fres, sint len, sint fseq, sint nseq,ALN mult_aln)
{
/*      static char *aacids = "XCSTPAGNDEQHRKMILVFYW";*/
/*		static char *nbases = "XACGT";	*/
	char *seq, residue;
	char val;
	sint i,ii;
	sint j,slen;	
	sint line_length;
	sint max_aln_length,max_names;

	max_aln_length=max_names=0;
	for(i=fseq; i<nseq; i++)  {
		if(mult_aln.seqs[i].len>max_aln_length) max_aln_length=mult_aln.seqs[i].len;
		if(strlen(mult_aln.seqs[i].name)>max_names) max_names=strlen(mult_aln.seqs[i].name);
	}

	seq = (char *)ckalloc((max_aln_length+1) * sizeof(char));
	
/* decide the line length for this alignment - maximum is LINELENGTH */
	line_length=PAGEWIDTH-max_names;
	line_length=line_length-line_length % 10; /* round to a multiple of 10*/
	if (line_length > LINELENGTH) line_length=LINELENGTH;

	for(ii=fseq; ii<nseq; ii++) {
		i = get_seq_index(mult_aln.seqs,nseq,ii);
		fprintf(nbout, mult_aln.dnaflag ? ">DL;" : ">P1;");
		fprintf(nbout, "%s\n%s\n", mult_aln.seqs[i].name, mult_aln.seqs[i].title);
		slen = 0;
		for(j=fres; j<fres+len; j++) {
			val = mult_aln.seqs[i].data[j];
			if(val == EOS)
				break;
			else if(!isalpha(val))
				residue = '-';
			else {
				residue = val;
			}
			seq[j-fres] = toupper(residue);
			slen++;
		}
		for(j=0; j<slen; j++) {
			fprintf(nbout,"%c",seq[j]);
			if(((j+1) % line_length == 0) || (j == slen-1)) 
				fprintf(nbout,"\n");
		}
		fprintf(nbout,"*\n");
	}	

	seq=ckfree((void *)seq);
}


void gde_out(FILE *gdeout, sint fres, sint len, sint fseq, sint nseq,Boolean lowercase,ALN mult_aln,SS_OPT ss_opt)
{
	char *seq, residue;
	char val;
	char *ss_mask1, *ss_mask2;
	sint i,ii;
	sint j,slen;	
	sint line_length;
	sint max_aln_length,max_names;

	max_aln_length=max_names=0;
	for(i=fseq; i<nseq; i++)  {
		if(mult_aln.seqs[i].len>max_aln_length) max_aln_length=mult_aln.seqs[i].len;
		if(strlen(mult_aln.seqs[i].name)>max_names) max_names=strlen(mult_aln.seqs[i].name);
	}

	seq = (char *)ckalloc((max_aln_length+1) * sizeof(char));

/* decide the line length for this alignment - maximum is LINELENGTH */
	line_length=PAGEWIDTH-max_names;
	line_length=line_length-line_length % 10; /* round to a multiple of 10*/
	if (line_length > LINELENGTH) line_length=LINELENGTH;

	if (mult_aln.prf1.ss.struct_penalties == SECST && ss_opt.use_ss1 == TRUE) {
		ss_mask1=calc_sec_struct_mask(mult_aln.seqs[0].len,mult_aln.prf1.ss.sec_struct_mask,ss_opt);
	}
	if (mult_aln.prf2.ss.struct_penalties == SECST && ss_opt.use_ss2 == TRUE) {
                ss_mask2=calc_sec_struct_mask(mult_aln.seqs[mult_aln.prf1.nseqs].len,mult_aln.prf2.ss.sec_struct_mask,ss_opt);


			}

	
	for(ii=fseq; ii<nseq; ii++) {
		i = get_seq_index(mult_aln.seqs,nseq,ii);
		fprintf(gdeout, mult_aln.dnaflag ? "#" : "%%");
		fprintf(gdeout, "%s\n", mult_aln.seqs[i].name);
		slen = 0;
		for(j=fres; j<fres+len; j++) {
			val = mult_aln.seqs[i].data[j];
			if(val==EOS)
				break;
			else if(!isalpha(val))
				residue = '-';
			else {
				residue = val;
			}
			if (lowercase)
				seq[j-fres] = (char)tolower((int)residue);
			else
				seq[j-fres] = (char)toupper((int)residue);
			slen++;
		}
		for(j=0; j<slen; j++) {
			fprintf(gdeout,"%c",seq[j]);
			if(( j!=0 && (j+1) % line_length == 0) || (j == slen-1)) 
				fprintf(gdeout,"\n");
		}
	}	

	if (ss_opt.output_struct_penalties == 0 || ss_opt.output_struct_penalties == 2) {
		if (mult_aln.prf1.ss.struct_penalties == SECST && ss_opt.use_ss1 == TRUE) {
			fprintf(gdeout,"\"SS_%-*s\n",max_names,mult_aln.prf1.ss.name);
			for(i=fres; i<fres+len; i++) {
				val=ss_mask1[i];
				if (!isalpha(val))
					seq[i-fres]='-';
				else
					seq[i-fres]=val;
			}
			seq[i-fres]=EOS;
			for(i=0; i<len; i++) {
				fprintf(gdeout,"%c",seq[i]);
				if((i % line_length == 0) || (i == len)) 
					fprintf(gdeout,"\n");
			}
		}
		
		if (mult_aln.prf2.ss.struct_penalties == SECST && ss_opt.use_ss2 == TRUE) {
			fprintf(gdeout,"\"SS_%-*s\n",max_names,mult_aln.prf2.ss.name);
			for(i=fres; i<fres+len; i++) {
				val=ss_mask2[i];
				if (!isalpha(val))
					seq[i-fres]='-';
				else
					seq[i-fres]=val;
			}
			seq[i]=EOS;
			for(i=0; i<len; i++) {
				fprintf(gdeout,"%c",seq[i-1]);
				if((i % line_length == 0) || (i == len)) 
					fprintf(gdeout,"\n");
			}
		}
	}
	if (ss_opt.output_struct_penalties == 1 || ss_opt.output_struct_penalties == 2) {
		if (mult_aln.prf1.ss.struct_penalties != NONE && ss_opt.use_ss1 == TRUE) {
			fprintf(gdeout,"\"GM_%-*s\n",max_names,mult_aln.prf1.ss.name);
			for(i=fres; i<fres+len; i++) {
				val=mult_aln.prf1.ss.gap_penalty_mask[i];
				if (!isalpha(val))
					seq[i-fres]='-';
				else
					seq[i-fres]=val;
			}
			seq[i]=EOS;
			for(i=0; i<len; i++) {
				fprintf(gdeout,"%c",seq[i]);
				if((i % line_length == 0) || (i == len)) 
					fprintf(gdeout,"\n");
			}
		}
		
		if (mult_aln.prf2.ss.struct_penalties != NONE && ss_opt.use_ss2 == TRUE) {
			fprintf(gdeout,"\"GM_%-*s\n",max_names,mult_aln.prf2.ss.name);
			for(i=fres; i<fres+len; i++) {
				val=mult_aln.prf2.ss.gap_penalty_mask[i];
				if (!isalpha(val))
					seq[i-fres]='-';
				else
					seq[i-fres]=val;
			}
			seq[i]=EOS;
			for(i=0; i<len; i++) {
				fprintf(gdeout,"%c",seq[i]);
				if((i % line_length == 0) || (i == len)) 
					fprintf(gdeout,"\n");
			}
		}
	}
		
	if (mult_aln.prf1.ss.struct_penalties == SECST && ss_opt.use_ss1 == TRUE) ckfree(ss_mask1);
	if (mult_aln.prf2.ss.struct_penalties == SECST && ss_opt.use_ss2 == TRUE) ckfree(ss_mask2);

	seq=ckfree((void *)seq);
}

/*
    calculate a sec structure mask from the secondary structures in the input file
    this is for display purposes only - the sec structure mask is not used during
    the alignment calculation. For the alignment, use calc_gap_penalty_mask().
*/
char * calc_sec_struct_mask(int prf_length, char *mask,SS_OPT ss_opt)
{
	int i,j;
	char *tmp_mask;

	tmp_mask = (char *)ckalloc((prf_length+10) * sizeof(char));
	for (i=0;i<prf_length;i++)
		tmp_mask[i] = mask[i];
	i=0;
	while (i<prf_length) {
		if (tolower(mask[i]) == 'a' || mask[i] == '$') {
			for (j = 0; j<ss_opt.helix_end_minus; j++) {
				if (i+j>=prf_length || (tolower(mask[i+j]) != 'a'
				                    && mask[i+j] != '$')) break;
				tmp_mask[i+j] = 'a';
			}
			i += j;
			while (tolower(mask[i]) == 'a'
				                    || mask[i] == '$') {
				if (i>=prf_length) break;
				if (mask[i] == '$') {
					tmp_mask[i] = 'A';
					i++;
					break;
				}
				else tmp_mask[i] = mask[i];
				i++;
			}
			for (j = 0; j<ss_opt.helix_end_minus; j++) {
				if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'a'
				                    || mask[i-j-1] == '$'))
					tmp_mask[i-j-1] = 'a';
			}
		}
	 	else if (tolower(mask[i]) == 'b' || mask[i] == '%') {
			for (j = 0; j<ss_opt.strand_end_minus; j++) {
				if (i+j>=prf_length || (tolower(mask[i+j]) != 'b'
				                    && mask[i+j] != '%')) break;
				tmp_mask[i+j] = 'b';
			}
			i += j;
			while (tolower(mask[i]) == 'b'
				                    || mask[i] == '%') {
				if (i>=prf_length) break;
				if (mask[i] == '%') {
					tmp_mask[i] = 'B';
					i++;
					break;
				}
				else tmp_mask[i] = mask[i];
				i++;
			}
			for (j = 0; j<ss_opt.strand_end_minus; j++) {
				if ((i-j-1>=0) && (tolower(mask[i-j-1]) == 'b'
				                    || mask[i-j-1] == '%'))
				tmp_mask[i-j-1] = 'b';
			}
		}
	else i++;
	}

	return tmp_mask;
}


/*
    calculate a gap penalty mask from the secondary structures in the input file
*/
void calc_gap_penalty_mask(int prf_length, char *ss_mask, char *gap_mask, SS_OPT ss_opt)
{
	int i;
	char *struct_mask;

/*
    first calculate a temp secondary structure mask 
*/
	struct_mask = calc_sec_struct_mask(prf_length,ss_mask,ss_opt);

/*
    then the gap penalty mask
*/
	for(i=0;i<prf_length;i++) {
		switch (struct_mask[i]) {
			case 'A':
				gap_mask[i] = ss_opt.helix_penalty+'0';
				break;
			case 'a':
				gap_mask[i] = ss_opt.helix_end_penalty+'0';
				break;
			case 'B':
				gap_mask[i] = ss_opt.strand_penalty+'0';
				break;
			case 'b':
				gap_mask[i] = ss_opt.strand_end_penalty+'0';
				break;
			default:
				gap_mask[i] = ss_opt.loop_penalty+'0';
				break;
		}
	}

        struct_mask=ckfree(struct_mask);

}