/*Secondary Structure Prediction Program  SIMPA96 */
/* Author Jonathan Levin E-Mail jon@glycine.jouy.inra.fr*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#define mxl 200
#define totres 150000
#define mxres 1500
#define mxseq 500
#define seg 15
#define cut 13
#define conf "Hb."
#define h 4
#define b 2
#define nchar 60
char dbs[] ="dbks.seq";
char dbo[] ="dbks.obs";
int readseq (char *, char [],char [][70], int [], int []); 
void find (char [],char [],int [],int ,int [][3],
char [][70],char *, int *[],int ,int []);
void regul(float (*)[3],int, char *);
float big(float *);
int ibig(float *);
float max(float,float);
void fmax(char *, int *, float *);
void consensus (char [][mxseq],int ,int ,float (*)[3],int *,
char *,char *);
void HH (char [][mxseq],int,int,float *);
void pseudo(char [][mxseq],int,int,int *, char *);
void sipred (char *,char *,int,char *,int *, int, FILE *);
int mat[676];
float dc[3]={.86,1.10,.81};
int main (void)
{
    FILE *fi1,*fi2,*fi3;
    int dird=0,i,j,k,kk,dnseq,dnres[mxseq],dcres[mxseq],tnseq,tnres[mxseq],
    tcres[mxseq],tusf,ali,l,ii,ncl,*ptseq[mxres],chk[mxres][3],mc,
    nrts[mxseq],scor[mxres],scs[mxres];
    float (*prop)[3],tprop[mxres][3],uprop[mxres][3],tmp[3];
    char dseq[totres],dobs[totres],dnames[mxseq][70],usf[70], 
    tseq[totres],tnames[mxseq][70],yes[3],ucf[70],line[mxl],
    clus[mxres][mxseq],prs[mxres],rpr[mxres*mxseq],mpr[mxres][mxseq],
    pr[mxres],lpr[mxres],file3[70],uof[70];
    char file1[70]="./";
    char file2[70]="/usr/local/lib/";
/* data entry */
    if((fi1=fopen("./mat.data","r"))== NULL)
    {
	if((fi1=fopen("/usr/local/lib/mat.data","r"))== NULL)
	{
	    printf(" can't open mat.data file from . or /usr/local/lib\n");
	    exit(1);
	}
	dird=1;
    }
    for(i=0;i<26;i++)
    {
	for(k=0;k<26;k++)fscanf(fi1,"%2d",&mat[i*26+k]);
	while(getc(fi1)!='\n');
    }
    fclose(fi1);
    if(dird)
    {
	strcpy(file3,file2);
        strcat(file3,dbs);
	dnseq=readseq(file3,dseq,dnames,dnres,dcres);
	strcpy(file3,file2);
	strcat(file3,dbo);
	dnseq=readseq(file3,dobs,dnames,dnres,dcres);
    }
    else
    {
	strcpy(file3,file1);
	strcat(file3,dbs);
	dnseq=readseq(file3,dseq,dnames,dnres,dcres);
	strcpy(file3,file1);
	strcat(file3,dbo);
	dnseq=readseq(file3,dobs,dnames,dnres,dcres);
    }
/* enter user files */
    printf(" enter sequence file name (PIR,fasta format)\n");
    gets(usf);
    tnseq=readseq(usf,tseq,tnames,tnres,tcres);
    if(tnseq<1)
    {
	printf(" no sequences found in file %s\n",usf);
	exit(1);
    }
    else
    for(i=0;i<tnseq;i++)
    if(tnres[i]<20)
    {
	printf(" sequence %s contains less than 20 AA\n",tnames[i]);
	exit(1);
    }
    printf(" %d sequences found in file %s\n",tnseq,usf);
    tusf=1;
    while(tusf)
    {
	printf(" is this correct? (y/n)\n");
	gets(yes);
	if(yes[0]=='y')tusf=0;
	else if(yes[0]=='n')exit(1);
    }
    yes[0]='\0';
    printf(" is there a clustal alignment file? (y/n)\n");
    gets(yes);
    ali=0;
    if(yes[0]=='y')
    {
	printf(" enter clustal file name \n");
	gets(ucf);
	ali=1;
	if((fi2=fopen(ucf,"r"))== NULL)
	{
	    printf(" can't open clustal file: %s\n",ucf);
	    exit(1);
	}
/* reads clustal file */
	j=0;k=0;
	for(i=0;i<199;i++)line[i]=' ';
	while((line[17]<'A'||line[17]>'Z')&&line[17]!='-')
	if(fgets(line,mxl,fi2)==NULL)
	{
	    printf("Error reading clustal file %s\n",ucf);
	    printf("Are you sure this file is ");
	    printf("in CLUSTAL W(1.60) format ?\n");
	    exit(1);
	}
	else k++;
	while(line[17]>='A'&&line[17]<='Z'||line[17]=='-')
	{
	    line[17]=' ';
	    if(fgets(line,mxl,fi2)==NULL)
	    {
		printf("Error reading clustal file %s\n",ucf);
		printf("Are you sure this file is ");
		printf("in CLUSTAL W(1.60) format ?\n");
		exit(1);
	    }
	    j++;
	}
	if(j!=tnseq)
	{
	    printf("ERROR no of seqs in seq file ");
	    printf("%d no of seq clustal file %d\n",tnseq,j);
	    exit(1);
	}
	rewind(fi2);
	for(i=0;i<k;i++)
	if(fgets(line,mxl,fi2)==NULL)
	{
	    printf("Error reading clustal file %s\n",ucf);
	    printf("Are you sure this file is ");
	    printf("in CLUSTAL W(1.60) format ?\n");
	    exit(1);
	}
	l=0;ii=0;
	while(l==0)
	{
	    for(i=0;i<tnseq;i++)
	    {
		j=ii*60;
		k=16;
		while(line[k]!='\n'){clus[j+k-16][i]=line[k];k++;}
		if(fgets(line,mxl,fi2)==NULL)
		{
		     printf("Error reading clustal file %s\n",ucf);
		     printf("Are you sure this file is ");
		     printf("in CLUSTAL W(1.60) format ?\n");
		     exit(1);
		}
	    }
	    if(fgets(line,mxl,fi2)==NULL)l=1;
	    line[16]=' ';
	    while((line[16]<'A'||line[16]>'Z')&&line[16]!='-')
	    if(fgets(line,mxl,fi2)==NULL){l=1;line[16]='-';}
	    ii++;
	    ncl=j+k-16;
	}
    fclose(fi2);
/* check to see if same number of residues in seq file and clus file */
    for(i=0;i<tnseq;i++)
    {
	k=tcres[i];
	for(j=0;j<ncl;j++)
	if(clus[j][i]!='-')
	if(clus[j][i]==tseq[k])k++;
	else
	{
	    printf("In file %s sequence no %3d res no. %4d '%c'\nis ",
	    usf,i+1,k-tcres[i]+1,tseq[k]);
	    printf("different from\nfile");
	    printf(" %s sequence no %3d res no. %4d '%c'\n",ucf,i+1,j+1,
	    clus[j][i]);
	    exit(1);
	}
    }
    
    }
    printf("enter prediction ouput file name\n");
    gets(uof);
    if((fi3=fopen(uof,"w"))== NULL)
    {
	printf("can't open output file: %s\n",uof);
	exit(1);
    }
    else printf(" opening prediction ouput file: %s\n", uof);
    fprintf(fi3,"SIMPA96 secondary structure prediction program\n\n");
    fprintf(fi3,"H = helix b = strand . = other\n\n");
    fprintf(fi3,"Scale indicates prediction strength: ");
    fprintf(fi3,"0-weak, 2-medium, 5-strong 9-very strong\n\n");
/* do prediction */
     prop=calloc(tcres[tnseq]*3,sizeof(float));
     for(ii=0;ii<tnseq;ii++)
     {
	 for(i=0;i<tnres[ii];i++)ptseq[i]=&mat[(tseq[tcres[ii]+i]-65)*26];
	 find (dseq,dobs,dnres,dnseq,chk,dnames,tnames[ii],ptseq,
	 tnres[ii],dcres);
	 for(i=0;i<tnres[ii];i++)
	 for(k=0;k<3;k++)prop[tcres[ii]+i][k]=chk[i][k]*dc[k];
	 printf("%s\n",tnames[ii]);
     }
/* regularize structure or do multiple alignment prediction */
     if(ali)
     {
	kk=0,mc=0;
	for(i=0;i<tnseq;i++)
	{
	    for(k=0;k<tnres[i];k++)for(j=0;j<3;j++)
	    tprop[k][j]=prop[mc+k][j];
	    regul(tprop,tnres[i],prs);
	    for(k=0;k<tnres[i];k++)rpr[mc+k]=prs[k];
	    mc+=tnres[i];
	}
	for(i=0;i<ncl;i++)for(k=0;k<3;k++)uprop[i][k]=0;
	nrts[0]=0;k=0;
	for(i=0;i<tnseq;i++)
	{
	    k+=tnres[i];
	    nrts[i+1]=k;
	}
	for(i=0;i<ncl;i++)
	{
	    for(j=0;j<tnseq;j++)
	    if(clus[i][j]!='-')
	    {
		for(k=0;k<3;k++)uprop[i][k]+=prop[nrts[j]][k];
		mpr[i][j]=rpr[nrts[j]];
		nrts[j]++;
	    }
	    else mpr[i][j]='-';
	    for(k=0;k<3;k++)tmp[k]=uprop[i][k];
	    fmax(&pr[i],&scor[i],tmp);
	}
	consensus(clus,tnseq,ncl,prop,tnres,pr,lpr);
	for(i=0;i<ncl;i++)
	if(clus[i][0]!='-')
	{
	    prs[kk]=pr[i];
	    scs[kk]=scor[i];
	    kk++;
	}
	for(i=0;i<tnseq;i++)
	for(j=0;j<25;j++)
	if(tnames[i][j]=='\n'||tnames[i][j]=='\0')tnames[i][j]=' ';
	fprintf(fi3,"multiple sequence prediction\n\n");
	kk=ncl/nchar;
	for(k=0;k<kk;k++)
	{
	    for(i=0;i<19;i++)putc(' ',fi3);
	    for(i=k*nchar+10;i<=k*nchar+nchar;i=i+10)fprintf(fi3,"%10d",i);
	    putc('\n',fi3);
	    for(j=0;j<16;j++)putc(tnames[0][j],fi3);fprintf(fi3,"   ");
	    for(j=k*nchar;j<k*nchar+nchar;j++)putc(clus[j][0],fi3);
	    putc('\n',fi3);
	    for(i=0;i<tnseq;i++)
	    {
		for(j=0;j<16;j++)putc(tnames[i][j],fi3);
		fprintf(fi3,"   ");
		for(j=k*nchar;j<k*nchar+nchar;j++)putc(mpr[j][i],fi3);
		putc('\n',fi3);
	    }
	    fprintf(fi3,"consensus pred     ");
	    for(j=k*nchar;j<k*nchar+nchar;j++)putc(pr[j],fi3);
	    putc('\n',fi3);
	    fprintf(fi3,"scale              ");
	    for(j=k*nchar;j<k*nchar+nchar;j++)fprintf(fi3,"%1d",scor[j]);
	    putc('\n',fi3);putc('\n',fi3);putc('\n',fi3);
	}
	if(kk*nchar!=ncl)
	{
	    for(i=0;i<19;i++)putc(' ',fi3);
	    for(i=kk*nchar+10;i<=ncl;i=i+10)fprintf(fi3,"%10d",i);
	    putc('\n',fi3);
	    for(j=0;j<16;j++)putc(tnames[0][j],fi3);fprintf(fi3,"   ");
	    for(j=kk*nchar;j<ncl;j++)putc(clus[j][0],fi3);putc('\n',fi3);
	    for(i=0;i<tnseq;i++)
	    {
		for(j=0;j<16;j++)putc(tnames[i][j],fi3);fprintf(fi3,"   ");
		for(j=kk*nchar;j<ncl;j++)putc(mpr[j][i],fi3);
		putc('\n',fi3);
	    }
	    fprintf(fi3,"consensus pred     ");
	    for(j=kk*nchar;j<ncl;j++)putc(pr[j],fi3);putc('\n',fi3);
	    fprintf(fi3,"scale              ");
	    for(j=kk*nchar;j<ncl;j++)fprintf(fi3,"%1d",scor[j]);
	    putc('\n',fi3);putc('\n',fi3);putc('\n',fi3);
	    fprintf(fi3,"consensus prediction\n");
	}
	sipred(tnames[0],tseq,tcres[0],lpr,scs,tnres[0],fi3);
     }
     else
     {
	 ii=0;mc=0;
	 fprintf(fi3,"secondary structure prediction without homologues\n\n");
	 for(i=0;i<tnseq;i++)
	 {
	     for(k=0;k<tnres[i];k++)
	     {
		 for(j=0;j<3;j++)tmp[j]=prop[ii][j];
		 ii++;
		 fmax(&prs[k],&scs[k],tmp);
	     }
	     for(k=0;k<tnres[i];k++)
	     for(j=0;j<3;j++)tprop[k][j]=prop[mc+k][j];
	     mc+=tnres[i];
	     regul(tprop,tnres[i],prs);
	     sipred(tnames[i],tseq,tcres[i],prs,scs,tnres[i],fi3);
	 }
     }
     fclose(fi3);
}
/* function sipred prints output seq,pred,scale */
void sipred (char *name,char *seq,int cres,char *pr,int *scale,int nres,
FILE *fi3)
{
    int j,k,kk;
    fprintf(fi3,"\n%s\n\n",name);
    kk=nres/nchar;
    for(k=0;k<kk;k++)
    {
	for(j=0;j<12;j++)putc(' ',fi3);
	for(j=k*nchar+10;j<=k*nchar+nchar;j=j+10)fprintf(fi3,"%10d",j);
	putc('\n',fi3);
	fprintf(fi3,"sequence    ");
	for(j=k*nchar;j<k*nchar+nchar;j++)putc(seq[j+cres],fi3);
	putc('\n',fi3);
	fprintf(fi3,"prediction  ");
	for(j=k*nchar;j<k*nchar+nchar;j++)putc(pr[j],fi3);
	putc('\n',fi3);
	fprintf(fi3,"scale       ");
	for(j=k*nchar;j<k*nchar+nchar;j++)fprintf(fi3,"%1d",scale[j]);
	putc('\n',fi3);putc('\n',fi3);putc('\n',fi3);
    }
    if(kk*nchar!=nres)
    {
	for(j=0;j<12;j++)putc(' ',fi3);
	for(j=kk*nchar+10;j<=nres;j=j+10)fprintf(fi3,"%10d",j);
	putc('\n',fi3);
	fprintf(fi3,"sequence    ");
	for(j=kk*nchar;j<nres;j++)putc(seq[j+cres],fi3);
	putc('\n',fi3);
	fprintf(fi3,"prediction  ");
	for(j=kk*nchar;j<nres;j++)putc(pr[j],fi3);putc('\n',fi3);
	fprintf(fi3,"scale       ");
	for(j=kk*nchar;j<nres;j++)fprintf(fi3,"%1d",scale[j]);
	putc('\n',fi3);putc('\n',fi3);
    }
}
/* function consensus tries to make a consensus prediction */
void consensus (char clus[][mxseq],int nseq,int ncl,float prop[][3],
int *nres,char *pr,char *lpr)
{
    char mpr[mxres],psu[mxres];
    float tprop[mxres][3],mprop[mxres][3],tmp[3],x,wei[mxseq];
    int i,j,k,l,nrts[mxres],maxp[mxres][3],norm[mxres],ici[26],t;
    HH(clus,nseq,ncl,wei);
    pseudo(clus,nseq,ncl,nres,psu);
    k=0;nrts[0]=0;
    for(i=0;i<nseq;i++)
    {
	k+=nres[i];
	nrts[i+1]=k;
    }
    for(i=0;i<ncl;i++)for(j=0;j<3;j++)tprop[i][j]=0;
    for(i=0;i<ncl;i++)for(j=0;j<3;j++)maxp[i][j]=0;
    for(i=0;i<ncl;i++)norm[i]=0;
    for(i=0;i<ncl;i++)
    {
	for(j=0;j<26;j++)ici[j]=0;t=0;
	for(j=0;j<nseq;j++)
	if(clus[i][j]!='-')ici[clus[i][j]-65]++;
	for(j=0;j<26;j++)if(ici[j]!=0)t++;
	for(j=0;j<nseq;j++)
	if(clus[i][j]!='-')
	{
	    x=0;
	    for(k=0;k<3;k++)tmp[k]=prop[nrts[j]][k];
	    k=ibig(tmp);
	    maxp[i][k]++;
	    for(k=0;k<3;k++)x+=tmp[k];
	    for(k=0;k<3;k++)tmp[k]=tmp[k]/x;
	    tmp[1]=tmp[1]*1.04;
	    for(k=0;k<3;k++)tprop[i][k]+=tmp[k]*wei[j];
	    nrts[j]++;
	    norm[i]++;
	}
    }
    for(i=0;i<ncl;i++)for(j=0;j<3;j++)
    tprop[i][j]=tprop[i][j]/(float)norm[i];
    j=0;
    for(i=0;i<ncl;i++)
    if(clus[i][0]!='-')
    {
	for(k=0;k<3;k++)mprop[j][k]=tprop[i][k];
	j++;
    }
    regul(tprop, ncl, pr);
    regul(mprop, nres[0], lpr);    
}
void HH (char clus[][mxseq],int nseq,int ncl,float *wei)
{
    int i,j,ici[26],t;
    float u;
    for(i=0;i<ncl;i++)wei[i]=0;
    for(i=0;i<ncl;i++)
    {
	for(j=0;j<26;j++)ici[j]=0;t=0;
	for(j=0;j<nseq;j++)
	if(clus[i][j]!='-')ici[clus[i][j]-65]++;
	for(j=0;j<26;j++)if(ici[j]!=0)t++;
	for(j=0;j<nseq;j++)
	if(clus[i][j]!='-')wei[j]+=1/(ici[clus[i][j]-65]*(float)t);
    }
    u=0;
    for(i=0;i<nseq;i++)u+=wei[i];
    for(i=0;i<nseq;i++)wei[i]=wei[i]/u;
}
void pseudo(char clus[][mxseq],int nseq,int ncl,int *nres, char *psu)
{
    int i, j, nco[mxseq];
    for(i=0;i<nseq;i++)nco[i]=0;
    for(i=0;i<ncl;i++)psu[i]='X';
    for(i=0;i<ncl;i++)
    {
	for(j=0;j<nseq;j++)
	{
	    if(clus[i][j]=='-')
	    {
		if(nco[j]!=0&&nco[j]<nres[j])psu[i]='-';
	    }
	    else nco[j]++;
	}
    }
}
/* calculates scale values and chosen conformation */
void fmax(char *c, int *sc, float tmp[])
{
int i,j;
float t=0,u=0;
for (i=0;i<3;i++)
if(tmp[i]>=t)
  {
    t=tmp[i];
    j=i;
  }
  if(j==0)*c='H';
  else if(j==1)*c='b';
  else if(j==2)*c='.';
  tmp[j]=0;
for (i=0;i<3;i++)
if(tmp[i]>=u)
  {
    u=tmp[i];
    j=i;
  }
 *sc=(1-u/t)*10;
 if(*sc>9)*sc=9;
}
/* regularisation using dynamic programming */
void regul(float (*prop)[3],int nres, char pr[])
{
  float dyn[mxres][3],H,B;
  int i,j,k,l;
  for(i=0;i<nres;i++)for(j=0;j<3;j++)dyn[i][j]=0;
/* start filling dyn first and last residues will be coil*/
  for(i=1;i<h+1;i++)dyn[i][0]=0;
  for(i=1;i<b+1;i++)dyn[i][1]=0;
  for(i=1;i<h+1;i++)dyn[h][0]+=prop[i][0];
  for(i=1;i<b+1;i++)dyn[b][1]+=prop[i][1];
  dyn[1][2]=prop[1][2];
  if(h>=b)
   {
     for(i=2;i<b+1;i++)dyn[i][2]=dyn[i-1][2]+prop[i][2];
     for(i=b+1;i<h+1;i++)
       {
         B=0;
         for(j=i-b+1;j<=i;j++)B+=prop[j][1];
         dyn[i][1]=max(big(&dyn[i-b][0])+B,dyn[i-1][1]+prop[i][1]);
         dyn[i][2]=big(&dyn[i-1][0])+prop[i][2];
       }
   }
  else
   {    
     for(i=2;i<h+1;i++)dyn[i][2]=dyn[i-1][2]+prop[i][2];
     for(i=h+1;i<b+1;i++)
       {
         H=0;
         for(j=i-h+1;j<=i;j++)H+=prop[j][0];
         dyn[i][0]=max(big(&dyn[i-h][0])+H,dyn[i-1][0]+prop[i][0]);
         dyn[i][2]=big(&dyn[i-1][0])+prop[i][2];
       }
   }
/* fill up dyn from max(h,b) to end */
  l=(h>=b)?h:b;
  for(i=l+1;i<nres-1;i++)
   {
     H=0;
     for(j=i-h+1;j<=i;j++)H+=prop[j][0];
     dyn[i][0]=max(big(&dyn[i-h][0])+H,dyn[i-1][0]+prop[i][0]);
     B=0;
     for(j=i-b+1;j<=i;j++)B+=prop[j][1];
     dyn[i][1]=max(big(&dyn[i-b][0])+B,dyn[i-1][1]+prop[i][1]);
     dyn[i][2]=big(&dyn[i-1][0])+prop[i][2];
   }
/* go back through dyn to choose best path */
  i=nres-2;
  k=ibig(&dyn[i][0]);
  pr[i]=conf[k];
  while(i>0)
   {
     if(i<0)
       {
         printf(" error at end of dyn **** i<0 ");
         exit(1);
       }
     i--;  
     j=ibig(&dyn[i][0]);
     if(j==k)pr[i]=conf[j];
     else if(k==0)
       {
         for(l=i;l>i-h+1;l--)pr[l]=conf[k];
         i=i-h+1;
         k=ibig(&dyn[i][0]);
         pr[i]=conf[k];
       }
     else if(k==1)
       {
         for(l=i;l>i-b+1;l--)pr[l]=conf[k];
         i=i-b+1;
         k=ibig(&dyn[i][0]);
         pr[i]=conf[k];
       }
     else if(k==2)
       {
         k=ibig(&dyn[i][0]);
         pr[i]=conf[k];
        }
    }
pr[0]='.';
pr[nres-1]='.';
}         
/* function big returns largest value in row of dyn */ 
float big(float dyn[])
{
  int i;
  float a=0;
  for(i=0;i<3;i++)if(a<dyn[i])a=dyn[i];
  return a;
}
/* function ibig returns index of biggest value */
int ibig(float dyn[])
{
  int i,j=0;
  float a=0;
  for(i=0;i<3;i++)
  if(a<dyn[i]){a=dyn[i];j=i;}
  return j;
}
/* function max returns biggest value */
float max(float x,float z)
{
return(x<z)?z:x;
}
/* fills up propensity matrix */
void find (char dseq[],char dobs[],int dnres[],int dnseq,int chk[][3],
char dnames[][70],char *tname, int *tseq[],int tres,int dcres[])
{
  int ii,i,j,ik,k,iseq[mxres],iobs[mxres],tsiz,dsiz,ijk,it;
  for(i=0;i<mxres;i++)for(j=0;j<3;j++)chk[i][j]=0;
  for(ii=0;ii<dnseq;ii++)
  {
    for(i=0;i<dnres[ii];i++)
      {
         iseq[i]=dseq[dcres[ii]+i]-65;
         if(dobs[dcres[ii]+i]=='H')iobs[i]=0;
         else if(dobs[dcres[ii]+i]=='E')iobs[i]=1;
         else if(dobs[dcres[ii]+i]=='C')iobs[i]=2;
      }
    if(tres<=dnres[ii])
      {
         tsiz=tres-seg-1;
	 dsiz=dnres[ii]-seg-1;
	 for(ik=0;ik<=dsiz;ik++)
	   {
	       tsiz++;
	       if(tsiz>=dsiz)tsiz=dsiz;
	       ijk=0;it=0;
	       for(k=0;k<seg;k++)ijk+=*(tseq[(it+k)]+iseq[ik+k]);
	       if(ijk>=cut)
	       for(k=0;k<seg;k++)chk[it+k][iobs[ik+k]]+=ijk;
	       for(j=ik;j<tsiz;j++)
	       {
		   ijk=ijk-*(tseq[it]+iseq[j])
		   +*(tseq[it+seg]+iseq[j+seg]);
		   it++;
		   if(ijk>=cut)
		   for(k=0;k<seg;k++)chk[it+k][iobs[j+1+k]]+=ijk;
	       }
	   }
	 tsiz=tres-seg;
	 for(ik=1;ik<tsiz;ik++)
	 {
	     ijk=0;it=0;
	     for(k=0;k<seg;k++)ijk+=*(tseq[(ik+k)]+iseq[it+k]);
	     if(ijk>=cut)
	     for(k=0;k<seg;k++)chk[ik+k][iobs[it+k]]+=ijk;
	     for(j=ik;j<tsiz;j++)
	     {
		 ijk=ijk-*(tseq[j]+iseq[it])
		 +*(tseq[j+seg]+iseq[it+seg]);
		 it++;
		 if(ijk>=cut)
		 for(k=0;k<seg;k++)chk[j+1+k][iobs[it+k]]+=ijk;
	     }
	 }
       }
     else
       {
         tsiz=tres-seg-1;
	 dsiz=dnres[ii]-seg-1;
	 for(ik=0;ik<=tsiz;ik++)
	 {
	     dsiz++;
	     if(dsiz>=tsiz)dsiz=tsiz;
	     ijk=0;it=0;
	     for(k=0;k<seg;k++)ijk+=*(tseq[(ik+k)]+iseq[it+k]);
	     if(ijk>=cut)
	     for(k=0;k<seg;k++)chk[ik+k][iobs[it+k]]+=ijk;
	     for(j=ik;j<dsiz;j++)
	     {
		 ijk=ijk-*(tseq[j]+iseq[it])
		 +*(tseq[j+seg]+iseq[it+seg]);
		 it++;
	         if(ijk>=cut)
		 for(k=0;k<seg;k++)chk[j+1+k][iobs[it+k]]+=ijk;
	     }
	 }
	 dsiz=dnres[ii]-seg;
	 for(ik=1;ik<dsiz;ik++)
	 {
	     ijk=0;it=0;
	     for(k=0;k<seg;k++)ijk+=*(tseq[(it+k)]+iseq[ik+k]);
	     if(ijk>=cut)
	     for(k=0;k<seg;k++)chk[it+k][iobs[ik+k]]+=ijk;
	     for(j=ik;j<dsiz;j++)
	     {
		 ijk=ijk-*(tseq[it]+iseq[j])
		 +*(tseq[it+seg]+iseq[j+seg]);
		 it++;
		 if(ijk>=cut)
		 for(k=0;k<seg;k++)chk[it+k][iobs[j+1+k]]+=ijk;
	     }
	 }
       }
  }
  return;
}
/* reads seq (or obs) files in pir or my format */ 
int readseq (char *file, char seq[],char names[][70],int nres[],int cres[])  
{
  FILE *fi1;
  char line[mxl];
  int i,j,k,nseq=0,cr=0;
  cres[0]=0; 
  if((fi1=fopen(file,"r"))== NULL)
     {
       printf(" cant open file %s\n",file);
       exit(1);
     }
  while (line[0]!='>'&&line[0]!='!')
  if(fgets(line,mxl-1,fi1)==NULL)
  {
      printf("file %s is not in FastA format\n",file);
      exit(1);
  }
  if(line[0]=='!')
    {
      rewind(fi1);
      while (fgets(line,mxl-1,fi1)!=NULL)
        {
          if(line[0]=='!')
          {
          i=0;j=0,k=0;
          while(i<68&&line[i+1]!='\n'){names[nseq][i]=line[i+1];i++;}
          names[nseq][i]='\n';names[nseq][i+1]='\0';
          while(k==0)
            {
              fgets(line,mxl-1,fi1);
              i=0;
              while(line[i]!='\n'&&i<199)
              if(line[i]>='A'&&line[i]<='Z')
                {
                  seq[cr+j++]=line[i];
                  i++;
                }
              else if (line[i]=='@')
                {
                   nres[nseq]=j;
		   cr+=j;
		   cres[nseq+1]=cr;
                   nseq++;
                   k=1;
                   i=200;
                }
              else i++;
            }
           }
        }
      fclose(fi1);
      return nseq;
    }
  else
    {
      rewind(fi1);
      while (fgets(line,mxl-1,fi1)!=NULL)
        {
          if(line[0]=='>')
          {
          i=0;j=0,k=0;
          while(i<68&&line[i+1]!='\n'){names[nseq][i]=line[i+1];i++;}
          names[nseq][i]='\n';names[nseq][i+1]='\0';
          while(k==0)
            {
              if(fgets(line,mxl-1,fi1)!=NULL)
                {
                  if (line[0]=='>')
                    {
                       nres[nseq]=j;
		       cr+=j;
		       cres[nseq+1]=cr;
                       nseq++;
                       i=0;j=0,k=0;
                       while(i<68&&line[i+1]!='\n')
                      {names[nseq][i]=line[i+1];i++;}
                       names[nseq][i]='\n';names[nseq][i+1]='\0';
                    }
                  else
                   {
		      i=0;
                      while(line[i]!='\n'&&i<199)
                      if(line[i]>='A'&&line[i]<='Z')
                        {
                          seq[cr+j++]=line[i];
                          i++;
                        }
                      else i++;
                   }
                }
              else
                {
                   nres[nseq]=j;
		   cr+=j;
		   cres[nseq+1]=cr;
                   nseq++;
                   k=1;
                   fclose(fi1);
                   return nseq;
                }
            }
          }
        }
    }
}