// *****************************************************************************
// Name: BranchAndTerminate.cpp
// By: Mark Voorhies
// On: 6/3/2003
// Time-stamp: <BranchAndTerminate.cpp 2003-06-19 11:24:35 Mark Voorhies>
// Impliments: BranchAndTerminate.h
// *****************************************************************************

#include "BranchAndTerminate.h"
#include "BoundingEnergy.h"
#include "InfDouble.h"

using namespace std;

int BoundingSingles(const DeeEnergy& energy, DeeTable& eliminated, 
		    const DeeSpace& space, const double& E_ref,
		    FILE *message_stream){
  int e_count = 0;
  vector<SeqPos> f(1);

  // Loop over positions

  for(unsigned int i_pos = 0; i_pos < space.NumPos(); ++i_pos){
    int num_valid = 0;

    //printf("Beginning bounding elimination for position %d\n", i_pos);

    // Loop over available resimers

    for(unsigned int r = 0; r < space.NumResimers(i_pos); ++r){
      if(eliminated.get(i_pos, r)){

	//printf("At resimer %d\n", r);	

	f[0] = SeqPos(i_pos, r);
	// Eliminate resimers for which E_bound > E_ref
	InfDouble e_temp(E_ref);

	if((e_temp = 
	    BoundingEnergy(energy, FixedSequence(f, &eliminated))) > E_ref){
	  eliminated.eliminate(i_pos, r);
	  ++e_count;
	  if(message_stream){
	    fprintf(message_stream,
		    "eliminated pos %s(%d) resimer %s(%d) (%s > %f)\n",
		    space.PosID(i_pos).c_str(), i_pos, 
		    space.ResimerID(i_pos, r).c_str(), r,
		    e_temp.String().c_str(), E_ref);
	  }
	}
	else{
	  //fprintf(message_stream, "Valid energy %f <= %f\n",
	  //	  e_temp, E_ref);
	  ++num_valid;
	}
      }
    }
    if(num_valid < 1){
      fprintf(stderr, "Error! No resimers remaining at %d\n", i_pos);
      abort();
    }
  }
  return e_count;
}

int BoundingDoubles(const DeeEnergy& energy, DeeTable& eliminated, 
		    const DeeSpace& space, const double& E_ref,
		    FILE *message_stream){
  int e_count = 0;
  vector<SeqPos> f(2);

  // Loop over pairs of positions

  for(unsigned int i_pos = 0; i_pos < space.NumPos(); ++i_pos){
    for(unsigned int j_pos = i_pos + 1; j_pos < space.NumPos(); ++j_pos){

      // Loop over available resimers

      for(unsigned int r = 0; r < space.NumResimers(i_pos); ++r){
	if(!eliminated.get(i_pos, r)){
	  continue;
	}
	f[0] = SeqPos(i_pos, r);
	for(unsigned int s = 0; s < space.NumResimers(j_pos); ++s){
	  if(!eliminated.get(i_pos, r, j_pos, s)){
	    continue;
	  }
	  //printf("i = %d r = %d j = %d s = %d\n", i_pos, r , j_pos, s);
	  //fflush(stdout);
	  f[1] = SeqPos(j_pos, s);

	  // Eliminate resimers for which E_bound > E_ref

	  InfDouble e_temp(E_ref);

	  if((e_temp = 
	      BoundingEnergy(energy, FixedSequence(f, &eliminated))) > E_ref){
	    eliminated.eliminate(i_pos, r, j_pos, s);
	    ++e_count;
	    if(message_stream){
	      fprintf(message_stream,
		      "%d,%d,%d,%d = DEP (%s > %f)\n",
		      i_pos, r, j_pos, s,
		      e_temp.String().c_str(), E_ref);
	    }
	  }
	}
      }
    }
  }
  return e_count;
}


int OptBoundingSingles(const DeeEnergy& energy, DeeTable& eliminated, 
		       const DeeSpace& space, const double& e_ref,
		       FILE *message_stream){
  int e_count = 0;
  
  // For each undetermined position, we will apply termination
  for(unsigned int k = 0; k < space.NumPos(); ++k){
    if(eliminated.AvailableResimers(k) <= 1){
      continue;
    }

    // Precalculate the component of the bounding energy that
    // is independent of the resimer at k

    vector<vector<InfDouble> > e_double(space.NumPos(),
					vector<InfDouble>(0));
    for(unsigned int i = 0; i < space.NumPos(); ++i){
      if(i == k){
	continue;
      }
      e_double[i].resize(space.NumResimers(i), 0);
      for(unsigned int r = 0; r < space.NumResimers(i); ++r){
	if(eliminated.get(i, r)){
	  // For each uneliminated i_r, i!=k

	  for(unsigned int j = 0; j < space.NumPos(); ++j){
	    if((j == k)||(j == i)){
	      continue;
	    }
	    // Find j_s that minimizes E_pair(i_r,j_s), excluding
	    // eliminated resimers and DEP's with i_r 
	    InfDouble min_s(1, true);
	    for(unsigned int s = 0; s < space.NumResimers(j); ++s){
	      if((eliminated.get(j,s))&&
		 (eliminated.get(i,r,j,s))){
		min_s = min(min_s, energy.pair_energy(i,r,j,s));

		// Note that we could have (eliminated.get(k,v,j,s) == false)
		// The worst case behavior that this would produce is a bounding
		// energy that is too low.  This is acceptable
	      }
	    }
	    e_double[i][r] += min_s;
	  }
	}
	else{
	  // Set eliminated resimers to infinite energy (redundant?)
	  e_double[i][r] = InfDouble(1,true);
	}
      }
    }

    // Calculate the bounding energy for each resimer at k, eliminating
    // resimers when the bounding energy exceeds the reference energy

    //printf("Beginning bounding elimination for position %d\n", k);
    for(unsigned int v = 0; v < space.NumResimers(k); ++v){
      if(!eliminated.get(k,v)){
	continue;
      }
      //printf("At resimer %d\n", v);
      InfDouble e_bound = 0.0;
      //printf("E(k,v) = %s\n", e_bound.String().c_str());
      for(unsigned int i = 0; i < space.NumPos(); ++i){
	if(i == k){
	  continue;
	}
	//printf("Adding position %d\n", i);
	InfDouble min_r(1, true); 
	for(unsigned int r = 0; r < space.NumResimers(i); ++r){
	  if((!eliminated.get(i,r))||
	     (!eliminated.get(i,r,k,v))){
	    //printf("Skipping %d,%d\n", i, r);
	    continue;
	  }
	  //printf("Checking resimer %d...%s ->", i, min_r.String().c_str());
	  min_r = min(min_r,
		      e_double[i][r]+double(2)*energy.pair_energy(i,r,k,v));
	  //printf("%s\n", min_r.String().c_str());
	}
	//printf("min_r(%d) = %s\n", i, min_r.String().c_str());
	e_bound += min_r;
      }
      if(e_bound > e_ref){
	if(message_stream){
	  fprintf(message_stream,
		    "eliminated pos %s(%d) resimer %s(%d) (%s > %f)\n",
		  space.PosID(k).c_str(), k, 
		  space.ResimerID(k, v).c_str(), v,
		  e_bound.String().c_str(), e_ref);
	eliminated.eliminate(k, v);
	++e_count;
	}
      }
      else{
	printf("keeping %s(%d) resimer %s(%d) (%s <= %f)\n",
	       space.PosID(k).c_str(), k, 
	       space.ResimerID(k, v).c_str(), v,
	       e_bound.String().c_str(), e_ref);
      }
    }
  }
  return e_count;
}

int OptBoundingDoubles(const DeeEnergy& energy, DeeTable& eliminated, 
		       const DeeSpace& space, const double& e_ref,
		       FILE *message_stream){
  int e_count = 0;

  // Lookup table for precalculation
  vector<vector<InfDouble> > e_double(space.NumPos(),
				      vector<InfDouble>(0));
  
  for(unsigned int i = 0; i < space.NumPos(); ++i){
    e_double[i].reserve(space.NumResimers(i));
  }

  // For each uneliminated pair, we will apply termination
  for(unsigned int k = 0; k < space.NumPos(); ++k){
    for(unsigned int m = k + 1; m < space.NumPos(); ++m){
      if((eliminated.AvailableResimers(k) <= 1)&&
	 (eliminated.AvailableResimers(m) <= 1)){
	// If there is only one choice at k and m, then we can't
	// do any eliminations.
	// (The stronger condition, that there is only one uneliminated
	//  pair common to both k and m, will produce the case we are
	//  testing for if logical elimination is applied).
	continue;
      }

      eliminated.LogicalPairs();

      // Precalculate the component of the bounding energy that
      // is independent of the resimers at k and m

      printf("Calculating lookup--table for %d,%d\n", k, m);
      for(unsigned int i = 0; i < space.NumPos(); ++i){
	if((i == k)||(i == m)){
	  continue;
	}
	e_double[i].resize(space.NumResimers(i), 0);
	for(unsigned int r = 0; r < space.NumResimers(i); ++r){
	  if(eliminated.get(i, r)){
	    // For each uneliminated i_r, i!=k or m

	    for(unsigned int j = 0; j < space.NumPos(); ++j){
	      if((j == k)||(j == m)||(j == i)){
		continue;
	      }
	      // Find j_s that minimizes E_pair(i_r,j_s), excluding
	      // eliminated resimers and DEP's with i_r 
	      InfDouble min_s(1, true);
	      for(unsigned int s = 0; s < space.NumResimers(j); ++s){
		if((eliminated.get(j,s))&&
		   (eliminated.get(i,r,j,s))){
		  min_s = min(min_s, energy.pair_energy(i,r,j,s));
		}
	      }
	      e_double[i][r] += min_s;
	    }
	    
	  }
	  else{
	    // Set eliminated resimers to infinite energy (redundant?)
	    e_double[i][r] = InfDouble(1,true);
	  }
	}
      }

      // Check for bad infinity flags
      for(unsigned int i = 0; i < space.NumPos(); ++i){
	if((i == k)||(i == m)){
	  continue;
	}
	for(unsigned int r = 0; r < space.NumResimers(i); ++r){
	  if(eliminated.get(i,r)&&
	     e_double[i][r].inf){
	    printf("Bad infinity flag! %d %d\n", i, r);
	  }
	}
      }

      // Calculate the bounding energy for each resimer pair at (k,m), 
      // eliminating resimers when the bounding energy exceeds the reference 
      // energy

      printf("eliminating");

      for(unsigned int v = 0; v < space.NumResimers(k); ++v){
	printf(".");
	if(!eliminated.get(k,v)){
	  continue;
	}
	for(unsigned int w = 0; w < space.NumResimers(m); ++w){
	  if((!eliminated.get(m,w))||
	     (!eliminated.get(k,v,m,w))){
	    continue;
	  }
	  InfDouble e_bound = energy.pair_energy(k,v,m,w)*double(2);
      
	  for(unsigned int i = 0; i < space.NumPos(); ++i){
	    if((i == k)||(i == m)){
	      continue;
	    }
	    InfDouble min_r(1, true); 
	    for(unsigned int r = 0; r < space.NumResimers(i); ++r){
	      if((!eliminated.get(i,r))||
		 (!eliminated.get(i,r,k,v))||
		 (!eliminated.get(i,r,m,w))){
		continue;
	      }
	      min_r = min(min_r,
			  e_double[i][r]+
			  (energy.pair_energy(i,r,k,v)+
			   energy.pair_energy(i,r,m,w))*double(2));
	    }
	    printf("min_r (i = %d k = %d v= %d m = %d w = %d): %s\n", 
		   i, k, v, m, w, min_r.String().c_str());
	    e_bound += min_r;
	  }
	  if(e_bound > e_ref){
	    if(message_stream){
	      fprintf(message_stream,
		      "eliminated pair (%d,%d,%d,%d) (%s > %f)\n",
		      m,w,k,v,
		      e_bound.String().c_str(), e_ref);
	    }
	    eliminated.eliminate(m, w, k, v);
	    ++e_count;
	  }
	}
      }
      printf("\n");
    }
  }
  return e_count;
}


/*

ExhaustiveSolution *BranchAndTerminate(const DeeEnergy& energy,
				       const DeeTable& eliminated,
				       const DeeSpace& space,
				       const int max_nodes,
				       FILE *message_stream){


  // ==================================================
  //  Step 1: Terminate rotamers at each position
  //          until no more rotamers can be terminated
  // ==================================================

  bool repeat = true;
  while(repeat){
    repeat = false;
    for(unsigned int i_pos = 0; i_pos < space.NumPos(); ++i_pos){
      if(Terminate(i_pos, energy, eliminated) > 0){
	repeat = true;
      }
    }
  }

  // ==================================================
  //  Step 2: Recursive Branch and Terminate search
  // ==================================================

  //  Initialize Branch & Bound variables

  int k = 0;  // Current depth
  bool init_flag = false; // True if E_bound and E_ref are initialized
  double E_bound = 0;
  double E_ref = 0;
  int count = 0, leaves = 0; // Nodes and leaves visited

  struct BandTchoice{
    unsigned int pos; // position
    unsigned int res; // resimer

    BandTchoice(const unsigned int pos_i = 0,
		const unsigned int res_i = 0) : pos(pos_i), res(res_i) {}
  };

  set<unsigned int> available_positions; // positions not in cur_path
  for(unsigned int i = 0; i < space.NumPos(); ++i){
    available_positions.insert(i);
  }
  // resimer choices at each position
  vector<vector<unsigned int> > S(space.NumPos()); 
  // Current path through the sequence positions
  vector<unsigned int> cur_path(space.NumPos());
  // Current set of resimer choices along the current path
  vector<unsigned int> a(space.NumPos());
  // Set of resimer choices for the best solution so far
  vector<unsigned int> best_solution(space.NumPos());
  // Path through the sequence positions for the best solution.
  vector<unsigned int> best_path(space.NumPos());

  // Choose the root position for the search
  cur_path[k] = ChoosePosition();
  available_positions.erase(cur_path[k]);
  
  for(unsigned int i = 0; i < space.NumResimers(cur_path[k]); ++i){
    if(eliminated.get(cur_pos, i)){
      S[k].push_back(BandTchoice(cur_pos, i));
    }
  }

  while(k > -1){

    // Loop over choices at current position

    while((k < k_max) &&     // S[k >= k_max].size() = 0
	  (S[k].size() > 0)
	  //  && E_bound < E_ref (impliment this!!!)
	  ){
      // Append a new rotamer to a

      a[k] = S[k].back();
      S[k].pop_back();
      ++count;

      if((message_stream)&&(count % 1000 == 0)){
      	fprintf(message_stream, "count = %d\n", count);
      }
      if((max_nodes > 0)&&(count > max_nodes)){
	fprintf(stderr, "max_nodes (%d) exceeded!\n", max_nodes);
	if(init_flag){
	  // Note: need to overload constructor to include path info
	  return new ExhaustiveSolution(best_energy, best_solution, &space,
					count, leaves, false);
	}
	else{
	  fprintf(stderr, "no solution found before exit!\n");
	  return 0;
	}
      }

      // update E_bound

      // check for solution
      if(k == k_max - 1){
	++leaves;
	// include path information
	cur_energy = energy.seq_energy(a);

	if((!init_flag)||(cur_energy < best_energy)){
	  best_energy = cur_energy;
	  // record solution
	  copy(a.begin(), a.end(), best_solution.begin());
	  copy(cur_path.begin(), cur_path.end(), best_path.begin());
	  init_flag = true;
	}

	E_ref = cur_energy; // make these the same variable?
      }
      ++k;
      
      // Generate S[k], applying termination

      if(k < k_max){
	for(set<unsigned int>::const_iterator i_pos = 
	      available_positions.begin();
	    i_pos != available_positions.end(); ++i_pos){
	  Terminate(*i_pos, energy, eliminated, 
		    a, cur_path.begin(), cur_path.begin()+k, 
		    available_positions);
	  }
	}
	cur_path[k] = ChoosePosition();
	available_positions.erase(cur_path[k]);

	// generate S[k]
	//printf("Generating S[%d]\n", k);
	for(unsigned int i = 0; i < space.NumResimers(cur_path[k]); ++i){
	  if(eliminated.get(cur_pos, i)){
	    // check this resimer for DEP's

	    bool e_flag = false;
	    for(int k2 = 0; k2 < k; ++k2){
	      if(!eliminated.get(k2, a[k2], k, i)){
		e_flag = true;
		break;
	      }
	    }
	    if(!e_flag){
	      S[k].push_back(i);
	    }
	  }
	}
      }
    }
    --k;
    // update E_bound
  }

  double cur_energy = 0.0, best_energy = 0.0;

  vector<int> a(space.NumPos()),    // Current solution
    best_solution(space.NumPos());  // Best solution so far
  vector<vector<int> > S(space.NumPos()); // Choices at current position


  // Initialize counter (number of nodes visited)

  // Initialize depth

  int k = 0, k_max = space.NumPos();

  // Begin backtrack search

  printf("Starting backtrack search\n");
      // advance
    
      // check for solution

      // descend

      ++k;

    }
    --k; // backtrack
  }

  if(!init_flag){
    fprintf(stderr, "All possible solutions eliminated!\n");
    return 0;
  }

  if(message_stream){
    fprintf(message_stream, "Visited %d nodes (%d leaves)\n", count, leaves);
  }
  // Done

  return new ExhaustiveSolution(best_energy, best_solution, &space,
				count, leaves, true);
}

*/