/* control.c++  Routines that control the entire process 
    of tracing, extending, and fixing etc.
*/
#include "drawwindow.h"
#include "pdbtolev.h"
#include "rankside.h"
#include "tordyn.h"
#include "trace.h" // includes a number of "extern .... " lines
#include "assignseq.h"
#include "peptide.h"
#include "extendfit.h"
#include "bone.h"


extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern int fitnum[NFIT];// see trace.c++ for definition
extern int verbose; //Defined in trace.c++  if = 1, print a lot of stuff
extern struct maiddatafile maidfile;
extern struct rotomer rot[30];// rotomer data for each amino acid
extern struct tordyn tordata;


void rerefine(int icontrol,struct griddata& den1,struct geometry& geo)
/*
  icontrol = 1:  rebuild side chains
  icontrol = 2:  refine fit using diamond force
  icontrol = 3:   rebuild and refine
  Refines the whole fit -
 First - runs build seq.
 Then refines using Diamond force
*/
{
   int irun,ifit,nmax,nbad,quitcond,val[2];
   int origmovecb,origendmove,origtorend,origusediamond;
   float origrfor,lastden;

   origmovecb = geo.movecb;
   geo.movecb = 0; // do not move Cb in refineonerot.
   tordata.useboneconstr = 0;  // if = 1 Use bone constraint on Ca (or N?) in dynamics
   origendmove = geo.endmove;
   geo.endmove = 6; // This seems to be better than default vallue (=4)
   origtorend = geo.torend;
   geo.torend = 6; // Number of residues whose density is averaged in avectermden
   origrfor = tordata.rfor;
   origusediamond = tordata.usediamond;
   nmax = 4;
//   tordata.usediamond = 1; // Use diamond force
   tordata.usediamond = 0; //Do Not Use diamond force
   geo.fnum = 1;  // starting residue in buildseq;
   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);
               goto quitrout;
          }
   }
   readfit(geo,maidfile.fitfile);
   if( (icontrol==1)||(icontrol==3) ){// rebuild side chain 
      for(ifit=1;ifit<=fitnum[0];ifit++){
        geo.quit = 0;
        if(fit[ifit].val1<0){
          printf(" ReBuilding residues for fit %d \n",ifit);
          geo.pickfitnum = ifit; // fit number that is used in buildseq
          geo.resnum = -fit[ifit].val1 ; // pdb sequence number
          val[0] =fit[ifit].val1; 
          val[1] =fit[ifit].val2; 
          buildseq(den1,geo,rot); // This returns the new fit
        }
      }
    }// icontrol
    if( (icontrol==2)||(icontrol==3) ){// refine - and rebuild every to times
      for(ifit=1;ifit<=fitnum[0];ifit++){
         geo.quit = 0;
         geo.fixend = 0; // run tordyn on all residues
         geo.oneside = 0;  // run on all sides
         val[0]= fit[ifit].val1;
         val[1]= fit[ifit].val2;     
         copyfittogeo(geo,fit[ifit]); //NOTE ASSUMES FIT IS NOT TOO LONG 
         if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
               tordata.adata[1] = 2;
//   srand48(7); // initialize randum number generator (drand48)
#if 0
         quitcond = 40; // run 30 rounds of refinement
//         tordata.rfor = 0.0;
         printf(" Refining fit %d using Diamond force and tordata.rfor/origrfor = %5.3f\n",
                  ifit,tordata.rfor/origrfor);
         tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);// "1" means use N term density as criteria for q;uiting
         avedenreg(den1,geo,geo.pos,&lastden,&avecb2,&nbad,geo.nend[0],
                  geo.nend[1]); // just applied to main chain atoms
         printf("lastden = %5.3f  nbad = %d\n",lastden,nbad);
#endif
#if 1 
         quitcond = 10; // run quitcond rounds of refinement
         for(irun = 0;irun<=nmax;irun++){ // 5 runs of 10 each
            tordata.rfor = (nmax-irun)*origrfor/nmax;
            if(tordata.usediamond==1)
                printf(" Refining fit %d using Diamond force and tordata.rfor/origrfor = %5.3f\n",
                  ifit,tordata.rfor/origrfor);
            else
                printf(" Refining fit %d NOT using Diamond force and tordata.rfor/origrfor = %5.3f\n",
                  ifit,tordata.rfor/origrfor);
            tordynbat(quitcond,den1,geo,tordata,&lastden,&nbad);// "1" means use N term density as criteria for q;uiting
            printf("  irun = %d lastden = %5.3f  nbad = %d\n",irun,lastden,nbad);
            if( (irun<nmax)&&(fit[ifit].val1<0) ){ // rebuild side chains if not last run
               printf(" ReBuilding residues for fit %d \n",ifit);
               copygeotofit(geo,fit[ifit],val); // Copy 
               geo.pickfitnum = ifit; // fit number that is used in buildseq
               geo.resnum = -fit[ifit].val1 ; // pdb sequence number
               buildseq(den1,geo,rot); // This returns the new fit
               copyfittogeo(geo,fit[ifit]); //NOTE ASSUMES FIT IS NOT TOO LONG 
               if(tordata.adata[1]>2) //This restarts interactive tor dyn each time called   
                    tordata.adata[1] = 2;
            }
            else if(irun==nmax)
               copygeotofit(geo,fit[ifit],val); // Copy 
         } // irun
#endif
         val[0]= fit[ifit].val1;
         val[1]= fit[ifit].val2;     
         copygeotofit(geo,fit[ifit],val); // Copy 
         savefit();
      } // refine
   } // icontrol
   savefit();
   tordata.usediamond=origusediamond;
   geo.torend = origtorend;
   geo.endmove=origendmove;
   geo.movecb = origmovecb;
   tordata.rfor = origrfor;
   quitrout:;
}

void retrace(struct griddata& den1,struct geometry& geo,struct bonedata& bone1)
/* Try tracing use less restrictive conditions
*/
{
   int origsavlen,origsetval;
   float origminavecb;
 
   origsavlen = geo.savlen;
   origsetval = geo.setval;
   origminavecb = geo.minavecb;
//   geo.minavecb = 0.5*geo.minavecb;
   printf(" Reducing geo.savlen by 1 from %d to %d\n",geo.savlen,geo.savlen-1);
   geo.savlen = geo.savlen -1;
   geo.setval = -1;
   geo.newfit = 2;  //if = 1, start from last bone pt, if = 2, use startpt = 1 -  Use old fit file for setting points.
//geo.newfit = 0; // USED ONLY FOR DEBUGGING
   printf(" Using geo.newfit = %d\n",geo.newfit);
   if( bone1.iexist == 0){
         printf("ERROR:  No bone file listed in %s\n",maidfile.name);
         exit(0);
   }
   readbonefile(bone1);
   tracepep(den1,geo,bone1);
   geo.savlen = origsavlen;
   geo.newfit = 0;
   geo.setval = origsetval;
   geo.minavecb=origminavecb;
}
