#include "../src-lib/dist.h"

/*  =======================                                  */
int CompCoor2To1DeselectMax
/*  =======================                                  */
(file_records * First, file_records * Second, int flag, bool SieveVerbosely)
{
  int ii, i1=0, i0=0;
  Statistic 
    s = {0} , *S = &s ;
  char 
    buf1 [BUFSIZ],
    buf2 [BUFSIZ];
  atom_record * WorstAtom ;

  EquivSecondToFirst(First,Second,flag);
  
  InitStats(S);
  S->name = "Distance" ;
  for (ii = 0;ii<Second->atomnum; ii++) {
    atom_record * a2 = Second->atoms[ii];
    if (!(a2->flags & flag)) {
      atom_record * a1 = a2->chaining ;
#     ifdef OFF
      /* The following must be true. The distances are set in EquivSecondToFirst. */
      double d = AtomDistance(a1,a2);
      if (d != a2->surface)
	STDERR("(d != a2->surface)");
#     endif

      if (a2->surface > S->max) {
	WorstAtom = a2 ;
      }
      IncrementStats(S,a2->surface);
    }
  }
  
  if (SieveVerbosely) {
    printf("CompCoor2To1DeselectMax():\n");
    printf("     1st         |    Dist    D/RMS  |      2nd\n");
    for (ii = 0;ii<Second->atomnum; ii++) {
      atom_record * a2 = Second->atoms[ii];
      if (!(a2->flags & flag)) {
	atom_record * a1 = a2->chaining ;
	double d = a2->surface ;
	sprintfAtomShortForm(buf2,a2);
	sprintfAtomShortForm(buf1,a1);
	printf("%s |  %6.2lf %1c %6.2lf  | %s\n",
	       buf1,d, (a2 == WorstAtom ? '*' : ' '), d/S->rms,buf2); 
      }
    }
    printf("Number-of-Atoms: %9.0f   RMS:           %9.4f\n",S->num,S->rms);
    printf("Max-Deviation:   %9.4f   Min-Deviation: %9.4f\n",S->max,S->min);
    printf("Avg-Deviation:   %9.4f   Num*RMS:       %9.4f\n",S->mean,S->num*S->rms);
  }
  
  SWITCH_ON(WorstAtom,FLAG_1_DESELECT|FLAG_2_DESELECT);
  SWITCH_ON(WorstAtom->chaining,FLAG_1_DESELECT|FLAG_2_DESELECT);
  return (int) S->num;
}


double SieveFit
(file_records *Moving, file_records *Fixed,
 double threshold, int MaxIters, bool SieveVerbosely)
{
  double rms, rms2;
  int i,NumAtoms;
  double ** RotMat = dmatrix(1,3,1,3) ;
  double *  TVect  = dvector(1,3);


  fprintf(stderr,"SieveFit(): threshold= %f  MaxIters= %d\n",threshold,MaxIters);
  
  for (i=0, rms=INFINITY;i<MaxIters && rms > threshold; i++) {
    rms = MyFit(RotMat,TVect,Moving,Fixed)  ;
    RotateCoordinates(Moving,RotMat);
    TranslateCoordinates(Moving,TVect);
    NumAtoms = CompCoor2To1DeselectMax(Fixed,Moving,FLAG_2_DESELECT, SieveVerbosely);
    fprintf(stderr," : Iteration= %4d  RMS= %10.6f  NumAtoms= %4d\n",i,rms,NumAtoms);
  }
}


main(int argc,char ** argv) 
{
  SAY_HI_FIRST;
  CITATION;
  do_initialization(0,0) ;
  {
    char 
      *MovingFileName = StringAfterOption(argc,argv,"Moving", "Moving.pdb"),
      *FixedFileName  = StringAfterOption(argc,argv,"Fixed", "Fixed.pdb");
    file_records
      *Moving = open_and_read_pdb(MovingFileName,0),
      *Fixed  = open_and_read_pdb(FixedFileName,0);
    double ** RotMat = dmatrix(1,3,1,3) ;
    double TVect[4];
    double rms ;
    
    FlagOffEverywhere(Fixed,FLAG_1_DESELECT|FLAG_2_DESELECT);
    FlagOffEverywhere(Moving,FLAG_1_DESELECT|FLAG_2_DESELECT);
    SieveFit(Moving,Fixed,
	     DoubleAfterOption(argc,argv,"Threshold",0.5),
	     (int)DoubleAfterOption(argc,argv,"MaxIters",200),
	     OPT("SieveVerbose"));
    if (OPT("ShowSelection")) {
      printf("SieveFit(): The final set of selected atoms is.... (* atom is included)\n");
      CompCoor2To1DeselectMax(Fixed,Moving,FLAG_2_DESELECT, true);
    } else if (!OPT("ShowSelection") || OPT("WriteMoved")) {
      write_pdb_file(stdout,Moving,0);
    }
  }			     
}