/* routines to rank side chains
 jan18,2000: Added routines for finding best distance between two fits
    that are to be connected by extendfit()
  rankaa is a routine I started to write to give a ranking to each
    a.a - NOt finished
 ANOTHER  simple procedure - simply store the value qual and nbad
    for each residue - go through sequence and sum up qual and 
     nbad for each possible starting point.
This approach has problem that small side chains will fit large side chain
    density.  Tried several procedures to correct for this
    (see rankside.c8sep99).  However, best seemed to be a simple
    conservative approach:
     For each residue
      Check if ring density (indicated by ibad[AA3] = 0;
        then if a.a. in seq. is not a ring, return 3=tvalue.
      require that tvalue/nsidechains < 0.3.
      This produces slight improvement.
   RESULT: IF 8 OR MORE RESIDUE - GET ALMOST PERFECT PREDICTION FOR
       endodm.
See rankside.notes for detail
*/
#include "drawwindow.h"  // calls global.h
#include "rankside.h"
#include "pdbtolev.h"
#include "peptide.h"
#include "assignseq.h"
#include "math.h"
#include "editbones.h"
#include "trace.h"

extern struct maiddatafile maidfile; // needed only for maidfile.denscale
extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern int fitnum[NFIT];// see trace.c++ for definition
extern struct extbones ext; // All the data for edited bones used to extend one end of a given fit
extern char*  nameres[30];


int readsequence(struct geometry& geo,char pdbsequence[])
/*
  The first line in file is the residue corresponding to first letter
     assigned to geo.seq[0]
 reads one letter sequence file ("pdbsequence") - assigns value 
   to geo.seq[i] (run time allocation)
    geo.seqnum = number of letters in sequence, 
    initially set to = -1: indicates sequence file
    not yet read.
  Returns 1 is sequcne successfully read, 0 if not
*/
{  
   FILE *fpseq;
   int count,resn,seq0;
   char reschar[6]; 

   if( (fpseq=fopen(pdbsequence,"r")) == NULL){
        return 0;
   }
   else{ // First just count letters in sequence
      fscanf(fpseq,"%d",&seq0); // seq0 = start number for pdb sequence
      count = seq0-1;
      if(count<0){
         printf("ERROR in readsequence - initial sequence number = %d\n",count);
         exit(0);
      }
      while(fscanf(fpseq,"%1s",reschar)!=EOF)/*this will skip over blanks*/
         count=count+1;
      geo.seqnum = count;
      geo.seq = new int[count+1]; // allocate seq
      fclose(fpseq);
      fpseq=fopen(pdbsequence,"r");// reopen and assign value
      fscanf(fpseq,"%d",&geo.seq[0]);
      count = geo.seq[0]-1;
      while(fscanf(fpseq,"%1s",reschar)!=EOF){/*this will skip over blanks*/
         count=count+1;
         addoneletter(reschar,&resn);
         geo.seq[count]=resn;
      }
      geo.maxseq = count;
      fclose(fpseq);
      return 1;
   } // else
}

void allocatefitqual(struct pepfit& fit)
/* Allocates fit.qualaa and fit.ibad that
    store qual and ibad for each of 20 a.a. for each fit residue
  NOTE: THIS ASSUMES fit is already allocated.
*/
{
   fmat2(fit.qualaa,fit.nalloc,VAL+5); // allocate new memory
   fmat2(fit.qualfirst,fit.nalloc,VAL+5); // allocate new memory
   fmat2(fit.quallast,fit.nalloc,VAL+5); // allocate new memory
   imat2(fit.ibad,fit.nalloc,VAL+5);
}

void freefitqual(struct pepfit& fit)
{
   freefmat2(fit.qualaa,fit.nalloc,VAL+5);// free memory array 3 dim array (mat)
   freefmat2(fit.qualfirst,fit.nalloc,VAL+5);// free memory array 3 dim array (mat)
   freefmat2(fit.quallast,fit.nalloc,VAL+5);// free memory array 3 dim array (mat)
   freeimat2(fit.ibad,fit.nalloc,VAL+5);// free memory array 3 dim array (mat)
}

void readresqual(char resqual[],int nfits[])
/* Reads residue quality file and assigns to value to fit
  Note: This routine must allocate qual and ibad since they are not
    allocated when fit is created.
  BE CAREFUL THAT NOT ALREADY ALLOCATED.
  Probably good idea to free them after use.
  nfit[0] = total # of fits in resqual
  nfit[i] = ifit for ith fit.
*/
{
   int i,resn,n0,truetype,rottype,ifit,resnum;
   float aveden;
   FILE *fp;

   if( (fp=fopen(resqual,"r")) == NULL){
         printf("ERROR: %s file cannot be found\n",resqual);
         exit(0);
   }
   fscanf(fp,"%d",&nfits[0]);
   if(nfits[0]>fitnum[0]){
      printf("ERROR: nfits>fitnum[0] in readresqual\n");
      exit(0);
   }
   if(nfits[0]>=MAXNF){
      printf("ERROR:  nfits > MAXNF in readresqual\n");
      exit(0);
   }
   for(i =1;i<= nfits[0];i++){
      fscanf(fp,"%d %d",&nfits[i],&n0); // start reading file resqual
      ifit = nfits[i];
      if(n0!=fit[ nfits[i] ].nres-1)
        printf("WARNING!!!! n0!=fit[ifit].nres in readresqual\n"); 
      allocatefitqual(fit[ifit]);
      for(resn=1;resn<=n0-1;resn++){
          fscanf(fp,"%d %d",&resnum,&truetype);
          if(resnum!=resn){
                printf("WARNING!!!! resnum( %d) !=resn (%d) in readresqual\n",resnum,resn); 
                printf("resn = %d, true type = %d %s:\n",resn,truetype,nameres[truetype]);
          }
          for(rottype = ALA; rottype<=VAL; rottype++){ // NOTE - new do not read AA#
#if 1 // read the detailed file  (qual = aveden since nbad not used in original )
             fscanf(fp,"%d %f %d %f %f %f",&rottype,&fit[ifit].qualaa[resn][rottype],
                &fit[ifit].ibad[resn][rottype],&aveden,&fit[ifit].qualfirst[resn][rottype],
                &fit[ifit].quallast[resn][rottype]);
#endif
#if 0 // read the short version
            fscanf(fp,"%d %f %d",&rottype,&fit[ifit].qualaa[resn][rottype],
                        &fit[ifit].ibad[resn][rottype]);
#endif
#if 0
          printf("   type = %d   %s:  qual = %5.3f nbad= %d aveden = %5.3f  first = %5.3f last = %5.3f\n",
            rottype,nameres[rottype],qual[resn][rottype][0],ibad[resn][rottype],
            qual[resn][rottype][1],qual[resn][rottype][2],qual[resn][rottype][3]);
#endif
#if 0
          printf("   type = %d   %s:  qual = %5.3f nbad= %d \n",
            rottype,nameres[rottype],fit[ifit].qualaa[resn][rottype],fit[ifit].ibad[resn][rottype]);
#endif
#if 0 // Do not need this because have compensate for bad ILE later on
          if(rottype==ILE){ // new - 13sep00 - because of problem with ILE, decrease ibad by 1
            if(fit[ifit].ibad[resn][rottype]>0)
               fit[ifit].ibad[resn][rottype]= fit[ifit].ibad[resn][rottype]-1;
          }
#endif
        }// rottype
      }// resn
    } // i
    fclose(fp);
}

int ringvalue(struct pepfit fit,int resn,int a1,int a2,int a3)
{
     int nbad,tbad;
     float iadd;
     float glnden,ringden,valden,aspden;

     valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
     ringden = (fit.qualaa[resn][PHE]+fit.qualaa[resn][TYR]+
                 fit.qualaa[resn][HIS] )/3.0;
     glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
     aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
     nbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];
     iadd= 0;
     if(nbad<=1){
       if(glnden<=1.2*ringden)
           iadd = 0;
       else
           iadd = 1;
     }
//     if(fit.qualaa[resn][ALA]> ringden)
//        iadd = iadd+1;
     if(valden>1.1*ringden) // change1 add 1.1
        iadd = iadd+1;
     if(aspden>1.2*ringden)
        iadd = iadd+1;
#if 0
     else if(aspden>ringden)
        iadd = iadd+1;
#endif
     tbad = fit.ibad[resn][a1] + fit.ibad[resn][a2] +fit.ibad[resn][a3];
     if(tbad==0)
         return 0+iadd;    
     if(tbad ==1)
         return 1+iadd;
     if(tbad==2)  
         return 2+iadd;
     if(tbad==3)  
         return 4+iadd;
     else if ((fit.ibad[resn][MET]==0)||(fit.ibad[resn][GLU]==0)||(fit.ibad[resn][GLN]==0) )
         return 3;
     else
         return 6;
}


int phevalue(struct pepfit fit,int resn,int a1,int a2,int a3)
// Idendical to ringvalue except use pheden for ringden
{
     int nbad;
     float iadd;
     float glnden,ringden,valden,aspden;

     valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
     ringden = fit.qualaa[resn][PHE];
     glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
     aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
     nbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];
     iadd= 0;
     if(nbad<=1){
       if(glnden<=1.2*ringden)
           iadd = 0;
       else
           iadd = 1;
     }
//     if(fit.qualaa[resn][ALA]> ringden)
//        iadd = iadd+1;
     if(valden>1.1*ringden) // change1 add 1.1
        iadd = iadd+1;
     if(aspden>1.2*ringden)
        iadd = iadd+1;
#if 0
     else if(aspden>ringden)
        iadd = iadd+1;
#endif
     if(fit.ibad[resn][a1]==0 ) 
         return 0+iadd;    
     else if( fit.ibad[resn][a2]==0 )
         return 1+iadd;
     else if( fit.ibad[resn][a3]==0 )
         return 1+iadd;
     else if ((fit.ibad[resn][MET]==0)||(fit.ibad[resn][GLU]==0)||(fit.ibad[resn][GLN]==0) )
         return 3;
     else
         return 6;
}

int tyrvalue(struct pepfit fit,int resn,int a1,int a2,int a3)
{
     int nbad;
     float iadd;
     float glnden,ringden,valden,aspden;

     valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
     ringden = fit.qualaa[resn][TYR];
     glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
     aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
     nbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];
     iadd= 0;
     if(nbad<=1){
       if(glnden<=1.2*ringden)
           iadd = 0;
       else
           iadd = 1;
     }
//     if(fit.qualaa[resn][ALA]> ringden)
//        iadd = iadd+1;
     if(valden>1.1*ringden) // change1 add 1.1
        iadd = iadd+1;
     if(aspden>1.2*ringden)
        iadd = iadd+1;
#if 0
     else if(aspden>ringden)
        iadd = iadd+1;
#endif
     if(fit.ibad[resn][a1]==0 ) 
         return 0+iadd;    
     else if( fit.ibad[resn][a2]==0 )
         return 1+iadd;
     else if( fit.ibad[resn][a3]==0 )
         return 1+iadd;
     else if ((fit.ibad[resn][MET]==0)||(fit.ibad[resn][GLU]==0)||(fit.ibad[resn][GLN]==0) )
         return 3;
     else
         return 6;
}

int hisvalue(struct pepfit fit,int resn,int a1,int a2,int a3)
{
     int nbad;
     float iadd;
     float glnden,ringden,valden;

     valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
     ringden = fit.qualaa[resn][HIS];
     glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
     nbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];
     iadd= 0;
     if(nbad<=1){
       if(glnden<=1.2*ringden)
           iadd = 0;
       else
           iadd = 1;
     }
     if(valden> 1.1*ringden) //change2
        iadd = iadd+1;
     if(fit.ibad[resn][a1]==0 ) 
         return 0+iadd;    
     else if( fit.ibad[resn][a2]==0 )
         return 1+iadd;
     else if( fit.ibad[resn][a3]==0 )
         return 1+iadd;
     else if ((fit.ibad[resn][MET]==0)||(fit.ibad[resn][GLU]==0)||(fit.ibad[resn][GLN]==0) )
         return 3;
     else
         return 6;
}

int glnvalue(struct pepfit fit,int resn)
{
   int iadd,nbad,nbadasp,iring;
   float glnden,aspden;

   iadd = 0;
   glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
   aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
   nbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];
   nbadasp = fit.ibad[resn][ASP]+fit.ibad[resn][ASN];
   iring = ringvalue(fit,resn,PHE,HIS,TYR);
   if(iring==0)
        iadd = iadd +2;
   else if(iring==1)
        iadd = iadd +1;
   if( glnden< 0.8*fit.qualaa[resn][ALA]) // change 3
        iadd = iadd+2;
   if(nbadasp==0){
      if(glnden<0.8*aspden)
        iadd = iadd+2;   
      else if(glnden< 0.9*aspden)
        iadd = iadd+1;   
   }
   if(nbad<=1){
     return nbad+iadd;
   }
   if(nbad==2){
//     return nbad+iadd+2; // OLD used up to 4oct00
       return nbad+iadd; // New used after to 4oct00
   }
#if 1
   if(2.0*fit.ibad[resn][SER]>=nbad)
      return 4;
#endif
   return 6;
}


int metvalue(struct pepfit fit,int resn)
/* similar to gln value but should have
   large asp den becuase that is where S is
*/
{
   int iadd,nbad,nbadasp,iring;
   float glnden,aspden,valden;

   iadd = 0;
   glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
   nbad = fit.ibad[resn][MET];
   nbadasp = fit.ibad[resn][ASP]+fit.ibad[resn][ASN];
   iring = ringvalue(fit,resn,PHE,HIS,TYR);
   if(iring==0)
        iadd = iadd +2;
   else if(iring==1)
        iadd = iadd +1;
   if( glnden<fit.qualaa[resn][ALA])
        iadd = iadd+2;
   if(nbadasp==0){
      if(aspden<valden)
        iadd = iadd+2;
   }
   if(nbad==0){
     return iadd;
   }
   if(nbad==1){
     return iadd+3;
   }
   return 6;
}

int aspvalue(struct pepfit fit,int resn)
{
   float glnden,aspden,valden;
   int glnbad,iadd;
   glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
   glnbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];

   if(phevalue(fit,resn,PHE,HIS,TYR)<=0)
      return 6;
   iadd = 0;
   if( (fit.ibad[resn][ASN]==0)||(fit.ibad[resn][ASP]==0) ||
              (fit.ibad[resn][LEU]==0)||(fit.ibad[resn][ILE]==0)  ){
      if(phevalue(fit,resn,PHE,HIS,TYR)==1)
         iadd = 1;
      if(aspden< 0.8*fit.qualaa[resn][ALA]) //change 4
         iadd = iadd+1;  
      if(aspden< 0.8*valden) // change5
         iadd = iadd+1;
      if(aspden< glnden){ //change6 - makes worse - changed back
        if(glnbad==0)
          iadd = iadd +2;
        else if(glnbad==1)
          iadd = iadd +1;
      }
      return iadd;
   }
   return 4;
}

int leuvalue(struct pepfit fit,int resn)
// Identical to aspvalue but us leuden
{
   float glnden,leuden,valden;
   int glnbad,iadd;
   glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   leuden = fit.qualaa[resn][LEU];
   glnbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];

   if(ringvalue(fit,resn,PHE,HIS,TYR)<=0)
      return 6;
   iadd = 0;
   if(fit.ibad[resn][LEU]==0){
      if(ringvalue(fit,resn,PHE,HIS,TYR)==1)
         iadd = 1;
      if(leuden< 0.8*fit.qualaa[resn][ALA]) //change 4
         iadd = iadd+1;  
      if(leuden< 0.8*valden) // change5
         iadd = iadd+1;
      if(leuden< glnden){ //change6 - makes worse - changed back
        if(glnbad==0)
          iadd = iadd +2;
        else if(glnbad==1)
          iadd = iadd +1;
      }
      return iadd;
   }
   return 4;
}

int ilevalue(struct pepfit fit,int resn)
{
   float glnden,aspden,valden;
   int glnbad,iadd;
   glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   aspden = fit.qualaa[resn][ASP];
   if(fit.qualaa[resn][ASN]>aspden)
      aspden = fit.qualaa[resn][ASN];
   if(fit.qualaa[resn][ILE]>aspden)
      aspden = fit.qualaa[resn][ILE];
   glnbad = fit.ibad[resn][GLN]+fit.ibad[resn][GLU];
   

   if(ringvalue(fit,resn,PHE,HIS,TYR)<=0)
      return 6;
   iadd = 0;
   
   if( (fit.ibad[resn][VAL]==0)||(fit.ibad[resn][THR]==0) ){
      if(ringvalue(fit,resn,PHE,HIS,TYR)==1)
         iadd = 2;
//      if(aspden<fit.qualaa[resn][ALA])
      if(valden< 0.8*fit.qualaa[resn][ALA]) // New 14sep00
         iadd = iadd+1; 
      if(aspden<glnden){
        if(glnbad==0)
          iadd = iadd +2;
        else if(glnbad==1)
          iadd = iadd +1;
      }
      return iadd;
   }
   return 4;
}

int valvalue(struct pepfit fit,int resn)
{
   int iadd;
   float valden;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
//   aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;

   if(ringvalue(fit,resn,PHE,TYR,HIS)<=1)
      return 6;
   if(glnvalue(fit,resn)<=1) 
      return 4;
   iadd = 0;
   if(valden < 0.75*fit.qualaa[resn][SER]) //change7 0.8 is worse than 0.7
     iadd = iadd+2;
   if( (fit.ibad[resn][VAL]==0)||(fit.ibad[resn][THR]==0) )
          return 0+iadd;
   if( (fit.ibad[resn][VAL]==1)||(fit.ibad[resn][THR]==1) )
          return 2+iadd;
   return 6;
}

int servalue(struct pepfit fit,int resn)
// New - 11dec00  qualify ringvalue by comparison with dnesith
{
   float serden,valden,ringden,aspden;
   ringden = (fit.qualaa[resn][PHE]+fit.qualaa[resn][HIS])/2.0;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   aspden = (fit.qualaa[resn][ASP]+fit.qualaa[resn][ASN])/2.0;
   serden = fit.qualaa[resn][SER];
      
   if(fit.ibad[resn][SER]==0){
      if(serden< ringden)
         return 6;
      if( (serden>1.8*valden)||(serden>1.8*aspden) )
         return 0;
      if( (serden>1.5*valden)||(serden>1.5*aspden) )
         return 2;
      if( (serden>1.0*valden)||(serden>1.0*aspden) )
         return 4;
      return 6;
  }
  return 2;  //ibad[ser]>0
}

int cysvalue(struct pepfit fit,int resn)
/* Similar to serval except removes aspvalue check and compares
    serden to alanineden  cystine den to alanine den - expect
    cystine den to be significantly greater.
*/
{
   float cysden,valden,aladen;
   aladen = fit.qualaa[resn][ALA];
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   cysden = fit.qualaa[resn][CYS];
#if 0
printf("ring = %d gln = %d asp = %d val= %d\n",
    ringvalue(fit,resn,PHE,HIS,TYR),glnvalue(fit,resn),aspvalue(fit,resn),
    valvalue(fit,resn) );
#endif
   if(ringvalue(fit,resn,PHE,HIS,TYR)<=1)
      return 6;
   if(ringvalue(fit,resn,PHE,HIS,TYR)<=2)
      return 4;
   if(glnvalue(fit,resn)==0) 
      return 6;
   if(glnvalue(fit,resn)==1) 
      return 4;
   if(glnvalue(fit,resn)==2) 
      return 3;
   if(fit.ibad[resn][CYS]==0 ){
      if(cysden> 1.5*aladen) 
        return 0;
      else
        return 2;
   }
   if( fit.ibad[resn][CYS]==1 ){
      if(cysden> 1.2*valden) 
        return 2;
      else
        return 4;
   }
   return 6;
}

int allbad(struct pepfit fit,int resn)
// New - 11dec00 - do not include Pro
// returns 1 if all residues at least one bad
{
   int i;
   for(i=ALA;i<=VAL;i++){
     if((i==GLY)||(i==PRO) )
       goto nexti;
     else{
       if(fit.ibad[resn][i]==0)
         return 0;
     }
     nexti:;
   }
   return 1;
}

int glyvalue(struct pepfit fit,int resn)
{
   float serden,valden;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   serden = (fit.qualaa[resn][SER]+fit.qualaa[resn][CYS])/2.0;
   if( fit.ibad[resn][ALA]>0) // New 15dec00
      return 0;
   if(fit.qualaa[resn][ALA]>2.0*valden)
      return 1; 
   if(ringvalue(fit,resn,PHE,HIS,TYR)<=1)
      return 6;
   if((glnvalue(fit,resn)<=1)&&(fit.ibad[resn][PHE]<=1) )
        return 6;
   if(aspvalue(fit,resn)<=1)
      return 4;
   if(valvalue(fit,resn)<=1)
      return 4;
   if( fit.ibad[resn][ALA]==0){
      if( fit.ibad[resn][SER]!=0)           
         return 0;
      else if(serden> 0.8*valden) // change10
         return 1;
      else
         return 2;
   }
   return 6;
}      
   
int alavalue(struct pepfit fit,int resn)
{
   float serden,valden;
   valden = (fit.qualaa[resn][VAL]+fit.qualaa[resn][THR])/2.0;
   serden = (fit.qualaa[resn][SER]+fit.qualaa[resn][CYS])/2.0;
   if( fit.ibad[resn][SER]>0){
      if(fit.ibad[resn][ALA]==0)
        return 0;
      else
        return 1;
   }
   if(ringvalue(fit,resn,PHE,HIS,TYR)<=1)
      return 6;
   if((glnvalue(fit,resn)<=1)&&(fit.ibad[resn][PHE]<=1) )
        return 6;
   if(aspvalue(fit,resn)<=1)
      return 4;
   if(valvalue(fit,resn)<=1)
      return 4;
#if 1
   if( fit.ibad[resn][ALA]==0){
      if(serden> 0.8*valden) // change10
         return 0;
      else if( fit.ibad[resn][SER]!=0)           
         return 0;
      else
         return 1;
   }
#endif
   if(allbad(fit,resn)==1)
//      return 0;
      return 1;  // new 11dec00 - return 1 - not perfect and use glyval
   return 6;
}      

int provalue(struct pepfit fit,int resn)
{
   if(fit.ibad[resn][PRO]==0){
      if(fit.qualaa[resn][PRO]> 0.7*fit.qualaa[resn][LEU])
         return 0;
      else if(fit.qualaa[resn][PRO]> 0.4*fit.qualaa[resn][LEU])
         return 1;
      else return 2;
   }
   return 6;
}

int lysvalue(struct pepfit fit,int resn)
/* modify on 9oct00 - return max value of 3 or 4
*/
{
  int iadd;
  float lysden = fit.qualaa[resn][LYS];
  float ringden = (fit.qualaa[resn][PHE]+fit.qualaa[resn][TYR]+
                 fit.qualaa[resn][HIS] )/3.0;
  float valden = (fit.qualaa[resn][VAL]+
                    fit.qualaa[resn][THR])/2.0;
  float aspden = (fit.qualaa[resn][ASP]+
                    fit.qualaa[resn][ASN])/2.0;
  float glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
  iadd = 0;
  if(lysden<ringden)
    iadd = iadd+2;
  else if(lysden< 1.1*ringden)
    iadd = iadd+1;
  if(fit.ibad[resn][LYS]>1)
    iadd = iadd +3;
  else if(lysden<0.8*valden){
    iadd = iadd+3;
  }
  else if(lysden<0.7*aspden){ // change14  0.8 to 0.7
    iadd = iadd+3;
  }
  else if(lysden<0.7*glnden){// change14  0.8 to 0.7
    iadd = iadd+3;  
  }
  else if(lysden< aspden){
    iadd = iadd+1;
  }
  else if(lysden< glnden){
    iadd = iadd+1;
  }
//printf(" lysvalue = %d\n",iadd);
  return iadd;
 }

int argvalue(struct pepfit fit,int resn)
{
  int iadd;
  float lysden = fit.qualaa[resn][LYS];
  float argden = fit.qualaa[resn][ARG];
  float valden = (fit.qualaa[resn][VAL]+
                    fit.qualaa[resn][THR])/2.0;
  float aspden = (fit.qualaa[resn][ASP]+
                    fit.qualaa[resn][ASN])/2.0;
  float glnden = (fit.qualaa[resn][GLN]+fit.qualaa[resn][GLU])/2.0;
  iadd = 0;
  if(argden<0.7*valden){ // change15
    iadd = iadd+4;
  }
  else if(argden<valden){
    iadd = iadd+2;
  }
  if(argden<0.7*aspden){
    iadd = iadd+4;
  }
  if(argden<aspden){
    iadd = iadd+2;
  }
  if(lysden<0.7*glnden){
    iadd = iadd+3;
  }
  else if(lysden<glnden){
    iadd = iadd+1;
  }
  if(fit.ibad[resn][ARG]>1)
    iadd = iadd +4;
  if(fit.ibad[resn][ARG]==1)
    iadd = iadd +2;
  return iadd;
 }

float findvalue(struct pepfit fit,int resn,int aatype)
/* assigns a numerical value for each residue, based on density:
   0:  indicates most likely to be correct
   1:  fair chance to be correct
   2,3 ..... 10 (no chance to be correct)
Input fit.qualaa[resn][ aatype ], and fit.ibad  and fit.qualfirst, fit.quallast
      where qualaa = aveden, qualfirst, quallast = aveden of first  and last part
   resn is the sidechain in the fit that is being tested
   aatype is the amino acid type in the sequence that might be assigned to this residue
First check for most distinctive parts in order: Rings, MET, etc.
*/
{
   float value;

   switch(aatype){
     case(PHE):
       value = phevalue(fit,resn,PHE,HIS,TYR);
       break;
     case(HIS):
       value = hisvalue(fit,resn,HIS,TYR,PHE);
       break;
     case(TYR):
       value = tyrvalue(fit,resn,TYR,PHE,HIS);
       break;
     case(TRP):
       if( (fit.qualaa[resn][TRP]>40)&&(fit.ibad[resn][TRP]==0) )
          value = 0;
       else if(ringvalue(fit,resn,PHE,HIS,TYR)<=1)
          value =2;
       else if ((fit.ibad[resn][MET]==0)||(fit.ibad[resn][GLU]==0)||(fit.ibad[resn][GLN]==0) )
         value = 3;
       else
         value = 6;
       break;
     case(MET): // essentially identical to GLN and GLU
       value = metvalue(fit,resn);
       break;
     case(GLN): //
       value = glnvalue(fit,resn);
       break;
     case(GLU): //
       value = glnvalue(fit,resn);
       break;
     case(ASP): //
       value = aspvalue(fit,resn);
       break;
     case(ASN): //
       value = aspvalue(fit,resn);
       break;
     case(LEU): //
       value = leuvalue(fit,resn);
       break;
     case(ILE): //
//       value = aspvalue(fit,resn);
       value = ilevalue(fit,resn);
       break;
     case(PRO):
       value = provalue(fit,resn);
       break;
     case(VAL):
       value = valvalue(fit,resn);
       break;
     case(THR):
       value = valvalue(fit,resn);
       break;
     case(SER):
       value = servalue(fit,resn);
       break;
     case(CYS):
       value = cysvalue(fit,resn);
       break;
     case(ALA):
       value = alavalue(fit,resn);
       break;
     case(GLY):
       value = glyvalue(fit,resn);
       break;
#if 0
     case(ARG):
       value = glnvalue(fit,resn)+2.0*fit.ibad[resn][ARG];
#endif
     case(ARG):
       value = lysvalue(fit,resn);
     case(LYS):
       value = lysvalue(fit,resn);
  }
  return value; 
}

void findqual(struct geometry& geo,int ifit,int seq0,float *tvalue)
{
   int iseq,resn,aatype,n0;
   float value;

    n0= fit[ifit].nres-1; // residues run from 0,1,2, .. n0, side chains from 1,n0-1
    *tvalue = 0;
    for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
        iseq= seq0+resn -1;
        aatype = geo.seq[iseq];
        value = findvalue(fit[ifit],resn,aatype);
        *tvalue =value + *tvalue; // the bigger tvalue, worse the fit
//printf("resn = %d  iseq = %d aatype = %d value = %5.3f\n",resn,iseq,aatype,value);
    }
    *tvalue = *tvalue*10.0/(n0-1);  //NEW - 18sep00, scale tvalue to value of 10 residues
}

void rankseq(struct geometry& geo,char resqual[])
/* 
THIS IS JUST USED FOR DEBUGGING findqual - called from fit.c++
routine for assigning sequence data to fits 
  Read data file resqual and assign values to fit
NOTE - specific for shad since starts test at seq0 = 12
*/
{
   int i,seq0,nfits[MAXNF],ifit,n0,count,nottop[MAXNF];
   int ii,ranki[MAXNF][MAXRANK+1],ranki2[MAXNF][MAXRANK+1],cor[30];
   float tvalue,rankf[MAXNF][MAXRANK+1],averank,avequal;
   int istart,iend,gapnum,tcount,seqbeg,seqend,seqend0;
   int irank0;
   float rankf0;

   readresqual(resqual,nfits); 
   if(geo.seqnum== -1){ // indicats pdb sequence not yet read - call readsequence
      if(readsequence(geo,maidfile.pdbsequence)!=1) //read into geo.seq indicates not successfully read
          printf("ERROR: %s cannot be found\n",maidfile.pdbsequence);
   }
   istart = 1;
   iend = nfits[0];
   seqbeg = geo.seq[0];
   seqend0 = -1; // if less than 0 than define  for each fit

#if 0  // Use this if just want to look at one fit = lookfit and/or one seq.
   int lookfit;
#if 0// Use this to look only at lookfit
   lookfit = 7;
   for(i=1;i<=nfits[0];i++){
     if(nfits[i]==lookfit){
       istart = iend =i;
       goto quitlookfit;
     }
   }
   quitlookfit:;
#endif
//   seqend0 = -1; // if less than 0 than define  for each fit
   seqend0 = 108; // then just look at this sequence
#endif

#if 0 // Use this for debugging findvalue - just prints the value for each reise
int truetype;
aatype = ARG;
//aatype = PHE;
printf("Testing Serine\n"); 
for(i=istart;i<=iend;i++){ // for each fit in X
 ifit = nfits[i];
 if( -fit[ifit].val1>0){ //NOW - this should always be the case - only write to resqual if assigned Only consider a fit if there is correct sequence value for it
   n0= fit[ifit].nres -1; //ASSUME THIS is correct
   for(resn=1;resn<=n0-1;resn++){// run through each residue with sidechains
     tcount=tcount+1;
     truetype= geo.seq[-fit[ifit].val1 + resn - 1]; // actual amino acid type
     if(truetype==aatype){
       printf("   ");
#if 0
printf("ring = %d gln = %d asp = %d val= %d  ser = %d\n",
    ringvalue(fit[ifit],resn,PHE,HIS,TYR),glnvalue(fit[ifit],resn),aspvalue(fit[ifit],resn),
    valvalue(fit[ifit],resn), servalue(fit[ifit],resn) );
#endif
     }
     tvalue = findvalue(fit[ifit],resn,aatype);
#if 1
     printf("fit %d res %d: value = %5.3f  Correct value = %s\n",
                    i,resn,tvalue,nameres[truetype]);
#endif
   }
 }
} // i
exit(0); 
#endif

tcount = 0;
printf("Running rankseq with istart = %d iend = %d geo.seq[0]=%d\n",istart,iend,geo.seq[0]);
for(i=istart;i<=iend;i++){ // for each fit in resqual
 ifit = nfits[i];
 if( -fit[ifit].val1>0){ // Only consider a fit if there is correct sequence value for it
   tcount=tcount+1;
   ranki[i][0]=0;
   n0= fit[ifit].nres -1; //ASSUME THIS is correct
   if(seqend0<0) // usual condition
      seqend = geo.seqnum - (n0-1);
   else{ // look at only sequence
      seqbeg=seqend = seqend0; 
   }
   for(seq0=seqbeg;seq0<=seqend;seq0++){ //shad starts with seq.# =12  pdb. seq. assigned to res 1
      tvalue = 0;
      findqual(geo,ifit,seq0,&tvalue);
      gapnum = 0;  // Not used here
      tvalue = -tvalue;
      rankqual(MAXRANK,ranki[i],ranki2[i],rankf[i],seq0,tvalue,gapnum);
   } // for seq0
 }// fit[i[].val>0
}
  count = 0; // counts number in top 20
  nottop[0] = 0; // counts number not in top, nottop[i] are fits not in top
  averank = 0;
  avequal = 0.0;
for(i=istart;i<=iend;i++){
  ifit = nfits[i];
  if( -fit[ifit].val1>0){ // Only consider a fit if there is correct sequence value for it
    cor[ifit]= -fit[ifit].val1;
    printf(" RANKING FOR FIT %d with nres = %d:  Correct iseq = %d\n",
            ifit,fit[ ifit ].nres-1,cor[ifit]);
    for(ii=1;ii<=ranki[i][0];ii++)
      printf("   seq. # = %d    qual = %5.3f\n",ranki[i][ii],rankf[i][ii]);
    irank0=0;  // counts ranking - testing if identical
    rankf0= -1000.0;
    for(ii=1;ii<=ranki[i][0];ii++){
      if(rankf[i][ii]!=rankf0){ // rankf if worse than prevous
         irank0 = irank0+1;
         rankf0=rankf[i][ii];
      }
      if(ranki[i][ii] == cor[ifit]){ // if correct seq
         averank = averank + irank0;
         avequal = avequal - rankf[i][ii];
         count = count+1;
         goto top20;
      }
    }
    nottop[0] = nottop[0]+1;
    nottop[ nottop[0] ] = i;
    top20:;
  }// fit[i[].val>0    
}
  averank = averank/count;
  avequal= avequal/count;
  printf(" Fraction in top 20 = %5.3f  ave rank = %5.3f ave quality = %5.3f\n",
             (float)count/tcount,averank,avequal);
  printf("Fits not in top:");
  for(i=1;i<=nottop[0];i++)
    printf("  %d",nottop[i]);
  printf("\n");
   for(i =1;i<= nfits[0];i++){
      ifit = nfits[i];
      freefitqual(fit[ifit]);
   }
}

void assignval1()
/* Routine that goes through each fit and assigns
  correct val1.
  Procedure - for each fit, find the pdb Ca that is closest
   to Ca1, - if within 1.5 Ang. - then sets fit[].val1 = - this value
*/
{
    int i,k,ip,is,ires,val1;
    float dist,mindist,fitpos[3],pdbpos[3];

    for(i=1;i<=fitnum[0];i++){
      mindist = 10000.0;
      for(k=0;k<=2;k++)
         fitpos[k] = fit[i].pos[1][1][k];
      for(ip=1;ip<=maidfile.npdb;ip++){ // each pdb file
        for(is=0;is<maidfile.pdb[ip].resdat[0][0][0];is++){ // all segments
          for(ires=1;ires<=maidfile.pdb[ip].resdat[is+1][0][0];ires++) { 
             for(k=0;k<=2;k++)
                pdbpos[k]=maidfile.pdb[ip].mainch[is][ires][1][k]; //mainch[0,1,2] = N,Ca,C
             dist = distv(fitpos,pdbpos);
             if(dist < mindist){
                mindist = dist;
                val1 = -atoi(maidfile.pdb[ip].resnumch[is][ires]);
             }
          }  
        }//is
      }//ip
      if(mindist < 1.5){
         fit[i].val1 = val1;
         printf(" Assigned fit %d  pdb val = %d\n",i,val1);
      }
      else{
         printf("Could not assign val1 to fit %d  mindist = %5.3f\n",i,mindist);
      }
    }//i
    savefit();
}

