palign.c

在任务级并行平台P2HP上

2005.3.2
***/
void writetofile(Boolean wfflag,char *filename,int si,int sj,int a,int b)
{
	FILE *to;
	if (wfflag) {
		if ((to=fopen(filename,"wb+"))==NULL)
		{info("Could not open result  file (%s) ",filename);
		exit(1) ; 
		}
	}
	else
	if ((to=fopen(filename,"ab+"))==NULL)
	{info("Could not open result  file (%s) ",filename);
	exit(1) ; 
	}
	//printf("a %d, b %d, si %d , sj %d",a,b,si,sj);
	//fputc(maxsore,to);
	/***fputc('#',to);
	fputc(a,to);
	fputc('#',to);
	fputc(b,to);
	fputc('#',to);
	fputc(tmat[si+1][sj+1],to);
	fputc('#',to);
	fputc(tmat[sj+1][si+1],to);
	fputc('#',to);***/
	fwrite(&a,sizeof(int),1,to);
	fwrite(&b,sizeof(int),1,to);
	if(fwrite(&tmat[si+1][sj+1],sizeof(double),1,to)!=1)
	{fprintf(stderr,"write error!")	;
	}
	if(fwrite(&tmat[sj+1][si+1],sizeof(double),1,to)!=1)
		fprintf(stderr,"write error");
	fclose(to);
}

/***
added by jf.
2005.3.11
***/
static void  read_paramfile(char * paramfilename)
{
	FILE *fp;
	//int lineno;
	char linebuf[30];
	//char fName[]="d:\\hpctest\\rtemp.txt";
	int temp=0;
	int lineno;
	char * str;
	
	if((fp=fopen(paramfilename,"r"))==NULL){
		fprintf(stderr,"Can't open %s\n",paramfilename);
		exit(-1);
	}
	for(lineno=1;fgets(linebuf,30,fp)!=NULL&&lineno<=4;++lineno)
	{//printf("%5d %s %d",lineno,linebuf,strlen(linebuf));
		int j;
		int num;
		for(j=0,num=0;j<strlen(linebuf)-1;j++)
		{
			temp=linebuf[j]-'0';
			num=num*10;
			num+=temp;
		}
		/****11.20**
		if(lineno==1) Asequence=num;
		if(lineno==2) Bsequence=num;
		if(lineno==3) Esequence=num;****/
		
		if(lineno==1) Asequence=num;
		if(lineno==2) Astart=num;
		if(lineno==3) Bsequence=num;
		if (lineno==4) Bend=num;
	}
	//**2005.11.20  fgets(linebuf,30,fp);
	fclose(fp);
	/*printf("%s",linebuf);
	for(int j=0;j<strlen(linebuf);j++)
	{
	temp=linebuf[j]-'0';
	num=num*10;
	num+=temp;
	}
	sum+=num;*/
	/****2005.11.20**
	if(linebuf[0]=='T') align_model=1;
	else align_model=0;****/
	
	
	//printf("Asequence: %d ,Bsequence: %d ,Esequence: %d ,align_model: %d \n",Asequence,Bsequence,Esequence,align_model);
	
}
/*
static void  read_paramfile(char * paramfilename)
{
FILE *fp;
//int lineno;
char linebuf[30];
//char fName[]="d:\\hpctest\\rtemp.txt";
int temp=0;

		if((fp=fopen(paramfilename,"r"))==NULL){
		fprintf(stderr,"Can't open %s\n",paramfilename);
		exit(-1);
		}
		for(int lineno=1;fgets(linebuf,30,fp)!=NULL&&lineno<=3;++lineno)
		{//printf("%5d %s",lineno,linebuf);
		for(int j=0,int num=0;j<strlen(linebuf)-1;j++)
		{
		temp=linebuf[j]-'0';
		num=num*10;
		num+=temp;
		}
		if (lineno==1) Asequence=num;
		if(lineno==2) Bsequence=num;
		if(lineno==3) Esequence=num;
		}
		fgets(linebuf,30,fp);
		fclose(fp);
		/*printf("%s",linebuf);
		for(int j=0;j<strlen(linebuf);j++)
		{
		temp=linebuf[j]-'0';
		num=num*10;
		num+=temp;
		}
sum+=num;*/
/*if(linebuf[0]=="T") align_model=1;
else align_model=0;

}*/



static void add(sint v)
{
	
	if(last_print<0) {
		displ[print_ptr-1] = v;
		displ[print_ptr++] = last_print;
	}
	else
		last_print = displ[print_ptr++] = v;
}

static sint calc_score(sint iat,sint jat,sint v1,sint v2)
{
	sint ipos,jpos;
	sint ret;
	
	ipos = v1 + iat;
	jpos = v2 + jat;
	
	ret=matrix[(int)seq_array[seq1][ipos]][(int)seq_array[seq2][jpos]];
	
	return(ret);
}


static float tracepath(sint tsb1,sint tsb2)
{
	char c1,c2;
    sint  i1,i2,r;
    sint i,k,pos,to_do;
	sint count;
	float score;
	char s1[600], s2[600];
	
	to_do=print_ptr-1;
	i1 = tsb1;
	i2 = tsb2;
	
	pos = 0;
	count = 0;
	for(i=1;i<=to_do;++i) {
		
		if (debug>1) fprintf(stdout,"%d ",(pint)displ[i]);
		if(displ[i]==0) {
			c1 = seq_array[seq1][i1];
			c2 = seq_array[seq2][i2];
			
			if (debug>0)
			{
				if (c1>max_aa) s1[pos] = '-';
				else s1[pos]=amino_acid_codes[c1];
				if (c2>max_aa) s2[pos] = '-';
				else s2[pos]=amino_acid_codes[c2];
			}
			
			if ((c1!=gap_pos1) && (c1 != gap_pos2) &&
				(c1 == c2)) count++;
			++i1;
			++i2;
			++pos;
		}
		else {
			if((k=displ[i])>0) {
				
				if (debug>0)
					for (r=0;r<k;r++)
					{
						s1[pos+r]='-';
						if (seq_array[seq2][i2+r]>max_aa) s2[pos+r] = '-';
						else s2[pos+r]=amino_acid_codes[seq_array[seq2][i2+r]];
					}
					
					i2 += k;
					pos += k;
			}
			else {
				
				if (debug>0)
					for (r=0;r<(-k);r++)
					{
						s2[pos+r]='-';
						if (seq_array[seq1][i1+r]>max_aa) s1[pos+r] = '-';
						else s1[pos+r]=amino_acid_codes[seq_array[seq1][i1+r]];
					}
					
					i1 -= k;
					pos -= k;
			}
		}
	}
	if (debug>0) fprintf(stdout,"\n");
	if (debug>0) 
	{
		for (i=0;i<pos;i++) fprintf(stdout,"%c",s1[i]);
		fprintf(stdout,"\n");
		for (i=0;i<pos;i++) fprintf(stdout,"%c",s2[i]);
		fprintf(stdout,"\n");
	}
	/*
	if (count <= 0) count = 1;
	*/
	score = 100.0 * (float)count;
	return(score);
}


static void forward_pass(char *ia, char *ib, sint n, sint m)
{
	
	sint i,j;
	pwint f,hh,p,t;
	
	maxscore = 0;
	se1 = se2 = 0;
	for (i=0;i<=m;i++)
    {
		HH[i] = 0;
		DD[i] = -g;
    }
	
	for (i=1;i<=n;i++)
	{
        hh = p = 0;
		f = -g;
		
        for (j=1;j<=m;j++)
		{
			
			f -= gh; 
			t = hh - g - gh;
			if (f<t) f = t;
			
			DD[j] -= gh;
			t = HH[j] - g - gh;
			if (DD[j]<t) DD[j] = t;
			
			hh = p + matrix[(int)ia[i]][(int)ib[j]];
			if (hh<f) hh = f;
			if (hh<DD[j]) hh = DD[j];
			if (hh<0) hh = 0;
			
			p = HH[j];
			HH[j] = hh;
			
			if (hh > maxscore)
			{
				maxscore = hh;
				se1 = i;
				se2 = j;
			}
		}
	}
	
}


static void reverse_pass(char *ia, char *ib)
{
	
	sint i,j;
	pwint f,hh,p,t;
	pwint cost;
	
	cost = 0;
	sb1 = sb2 = 1;
	for (i=se2;i>0;i--)
    {
		HH[i] = -1;
		DD[i] = -1;
    }
	
	for (i=se1;i>0;i--)
	{
        hh = f = -1;
        if (i == se1) p = 0;
        else p = -1;
		
        for (j=se2;j>0;j--)
		{
			
			f -= gh; 
			t = hh - g - gh;
			if (f<t) f = t;
			
			DD[j] -= gh;
			t = HH[j] - g - gh;
			if (DD[j]<t) DD[j] = t;
			
			hh = p + matrix[(int)ia[i]][(int)ib[j]];
			if (hh<f) hh = f;
			if (hh<DD[j]) hh = DD[j];
			
			p = HH[j];
			HH[j] = hh;
			
			if (hh > cost)
			{
				cost = hh;
				sb1 = i;
				sb2 = j;
				if (cost >= maxscore) break;
			}
		}
        if (cost >= maxscore) break;
	}
	
}

static int diff(sint A,sint B,sint M,sint N,sint tb,sint te)
{
	sint type;
	sint midi,midj,i,j;
	int midh;
	static pwint f, hh, e, s, t;
	
	if(N<=0)  {
		if(M>0) {
			del(M);
		}
		
		return(-(int)tbgap(M));
	}
	
	if(M<=1) {
		if(M<=0) {
			add(N);
			return(-(int)tbgap(N));
		}
		
		midh = -(tb+gh) - tegap(N);
		hh = -(te+gh) - tbgap(N);
		if (hh>midh) midh = hh;
		midj = 0;
		for(j=1;j<=N;j++) {
			hh = calc_score(1,j,A,B)
				- tegap(N-j) - tbgap(j-1);
			if(hh>midh) {
				midh = hh;
				midj = j;
			}
		}
		
		if(midj==0) {
			del(1);
			add(N);
		}
		else {
			if(midj>1)
				add(midj-1);
			displ[print_ptr++] = last_print = 0;
			if(midj<N)
				add(N-midj);
		}
		return midh;
	}
	
	/* Divide: Find optimum midpoint (midi,midj) of cost midh */
	
	midi = M / 2;
	HH[0] = 0.0;
	t = -tb;
	for(j=1;j<=N;j++) {
		HH[j] = t = t-gh;
		DD[j] = t-g;
	}
	
	t = -tb;
	for(i=1;i<=midi;i++) {
		s=HH[0];
		HH[0] = hh = t = t-gh;
		f = t-g;
		for(j=1;j<=N;j++) {
			if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
			if ((hh=HH[j]-g-gh) > (e=DD[j]-gh)) e=hh;
			hh = s + calc_score(i,j,A,B);
			if (f>hh) hh = f;
			if (e>hh) hh = e;
			
			s = HH[j];
			HH[j] = hh;
			DD[j] = e;
		}
	}
	
	DD[0]=HH[0];
	
	RR[N]=0;
	t = -te;
	for(j=N-1;j>=0;j--) {
		RR[j] = t = t-gh;
		SS[j] = t-g;
	}
	
	t = -te;
	for(i=M-1;i>=midi;i--) {
		s = RR[N];
		RR[N] = hh = t = t-gh;
		f = t-g;
		
		for(j=N-1;j>=0;j--) {
			
			if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
			if ((hh=RR[j]-g-gh) > (e=SS[j]-gh)) e=hh;
			hh = s + calc_score(i+1,j+1,A,B);
			if (f>hh) hh = f;
			if (e>hh) hh = e;
			
			s = RR[j];
			RR[j] = hh;
			SS[j] = e;
			
		}
	}
	
	SS[N]=RR[N];
	
	midh=HH[0]+RR[0];
	midj=0;
	type=1;
	for(j=0;j<=N;j++) {
		hh = HH[j] + RR[j];
		if(hh>=midh)
			if(hh>midh || (HH[j]!=DD[j] && RR[j]==SS[j])) {
				midh=hh;
				midj=j;
			}
	}
	
	for(j=N;j>=0;j--) {
		hh = DD[j] + SS[j] + g;
		if(hh>midh) {
			midh=hh;
			midj=j;
			type=2;
		}
	}
	
	/* Conquer recursively around midpoint  */
	
	
	if(type==1) {             /* Type 1 gaps  */
		diff(A,B,midi,midj,tb,g);
		diff(A+midi,B+midj,M-midi,N-midj,g,te);
	}
	else {
		diff(A,B,midi-1,midj,tb,0.0);
		del(2);
		diff(A+midi+1,B+midj,M-midi-1,N-midj,0.0,te);
	}
	
	return midh;       /* Return the score of the best alignment */
}

static void del(sint k)
{
	if(last_print<0)
		last_print = displ[print_ptr-1] -= k;
	else
		last_print = displ[print_ptr++] = -(k);
}