/* EGAD: HQM_rotamers.cpp
   
Navin Pokala and Tracy Handel
Dept. of Molecular and Cell Biology
University of California, Berkeley
Copyright (C) 2003 Regents of the University of California
GNU Public License
Aug 12 2003

Absolutely no warranties are made or are implied with the use of this program or its parts.

This file contains the function HQM_rotamers
*/



#include "HQM_rotamers.h"

int has_this_position_been_used(int k, int *already_used_positions)
{
	int i;
		
	i=1;
	while(already_used_positions[i]!=ENDFLAG)
	{
		if(k==already_used_positions[i])
			return(1);
		++i;
	}	
	return(0);
}

/* exhaustive downhill line minimization for the sidechains defined by CHROMOSOME chr; used as a quench at the 
    end of GA or MC. chr must be soluble.
*/
/*
     Select a position at random
	Place every residue-rotamer at that position,  evaluate energy, and keep the best
	Repeat until no more downhill moves occur.
      */
      
void HQM_rotamers(PARAMETERS *parameters, CHROMOSOME *chr, VARIABLE_POSITION *var_pos)
{
   int neighborOptimizeFlag = 0;
   int GenX = 0;
   GENE i_gene;
   int i, k, previous_seqpos, num_moving_positions;
   int best_res, best_res_rot;
   static int *indicies_of_floating_positions=NULL, *already_used_positions=NULL;
   static int max_residues=0;
   CHROMOSOME dummychr, dummychr2;
   double best_energy, current_energy, fabs_dE;
	extern int SOLUBILITY_CUTOFF_FLAG;
	double fixed_sasa, fixed_hphob_sasa, fixed_charge, fixed_E_transfer;
	double sasa_other, sasa_hphob_other, charge_other, E_transfer_other;
   int num_isoenergy_safety, safety_convergence;

   
	if(is_this_chr_soluble(chr)==0)
	    {
		failure_report("ERROR HQM_rotamers.cpp: HQM_rotamers assumes that the chr input is soluble","abort");

		/*
	    	fprintf(stderr,"ERROR HQM_rotamers.cpp: HQM_rotamers assumes that the chr input is soluble\n");
		abort();
		*/
	    }


   best_res = ENDFLAG; best_res_rot = ENDFLAG;
   
 
   if(max_residues< MAX_RESIDUES)
   {
	max_residues=MAX_RESIDUES;
	if(indicies_of_floating_positions != NULL)
	{
		free_memory(indicies_of_floating_positions);
		free_memory(already_used_positions);
	}
	indicies_of_floating_positions = (int *)calloc(MAX_RESIDUES,sizeof(int));
	already_used_positions = (int *)calloc(MAX_RESIDUES,sizeof(int));
   }

	for(i=1;i<=parameters->numMovingPositions+1;++i)
		already_used_positions[i] = ENDFLAG;

	i=1; k=0; num_moving_positions=0;
   	while(var_pos[i].seq_position!=ENDFLAG)
   	{
   		 if(var_pos[i].fixed_flag == 0)  /* not fixed, so add to the list */
		 if(var_pos[i].dead_ended_flag == 0)  /* not dead-ended, so add to the list */
		 {
			++k;
			++num_moving_positions;
			indicies_of_floating_positions[k] = i;
		 }
		 ++i;
   	}
    
	if(num_moving_positions == 0)
	{
		randomize_sidechains(chr,0);
		CHROMOSOME_to_lookupEnergy(chr,&parameters->fixedatoms_energy);       
		return;	
	}

	inoculate_sidechains(&dummychr, var_pos,0);  
  inoculate_sidechains(&dummychr2, var_pos,0); 
   dummychr.bkbngenes = NULL;
   dummychr.first_bkbngene = NULL;
   dummychr2.bkbngenes = NULL;
   dummychr2.first_bkbngene = NULL;

	parameters->hqm_convergence = num_moving_positions;
   	safety_convergence = 1000*num_moving_positions;
     


     i_gene = (MENDEL *)malloc(sizeof(MENDEL));
     
	copyCHROMOSOME(dummychr, (*chr));  /* dummychr is the working chromosome */
	copyCHROMOSOME(dummychr2, (*chr));  /* dummychr2 is used for evaluating energies  */



	/* calculate fixed_postion contributions to sasa and charge */
 	if(SOLUBILITY_CUTOFF_FLAG == 1)
		fixed_position_sasa_charge(&dummychr, &fixed_sasa, &fixed_hphob_sasa, &fixed_charge, &fixed_E_transfer);

     GenX=0; previous_seqpos=0; num_isoenergy_safety=0;
     while(GenX<parameters->hqm_convergence)
     { 


	/* pick a random position */
	k=1;
	if(num_moving_positions>1)
	{
		k = randint(num_moving_positions);
		while(has_this_position_been_used(k,already_used_positions)==1 || k==0)
			k = randint(num_moving_positions);
		i=1;
		while(already_used_positions[i]!=ENDFLAG)
			++i;
		already_used_positions[i]=k;

	}

	i = indicies_of_floating_positions[k];
	previous_seqpos = k;

	 dummychr.genes = dummychr.firstgene; 
	 while(dummychr.genes->seq_position!=var_pos[i].seq_position)
	 {
	    dummychr.genes = dummychr.genes->nextgene;
	 }   
	 copyGENE(i_gene,dummychr.genes); 
	 neighborOptimizeFlag = 1;   

	/* current energy at this position */
	current_energy = lookup_energy_sasa_charge_other(i, i_gene, dummychr.firstgene, &sasa_other, &sasa_hphob_other, &charge_other, &E_transfer_other) + iROT->energy_var_fix;

	/*
	if(DEBUG == 1)
		printf("before:\t%d\t%lf\t%d\t%d\n",i,current_energy,i_gene->j_choice_index+1,i_gene->lookupRot_index);
	 */

	/* place all rotamers at this position and find the best one */
	best_energy = find_best_res_res_rot(i, var_pos, dummychr.firstgene, fixed_sasa, fixed_hphob_sasa, fixed_charge, fixed_E_transfer, &best_res, &best_res_rot);
 
/*
	if(DEBUG == 1)
	{
		printf("after:\t%d\t%lf\t%d\t%d\t%lf\n",i,best_energy,best_res,best_res_rot,best_energy-current_energy);
		if(fabs(best_energy-current_energy)>EPS)
		{
			DEBUG=2;

			printf("\t\tbefore:\n");
			lookup_energy_sasa_charge_other(i, i_gene, dummychr.firstgene, &sasa_other, &sasa_hphob_other, &charge_other, &E_transfer_other);

			CHROMOSOME_to_lookupEnergy(&dummychr, &parameters->fixedatoms_energy);
			
			printf("total_energy = %lf\n",dummychr.energy);

			i_gene->choice_ptr = &var_pos[i].choice[best_res];
			i_gene->j_choice_index = best_res-1;
			i_gene->lookupRot_index = best_res_rot;
			i_gene->chi = i_gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[i_gene->lookupRot_index].chi;
			i_gene->lookupRot = &(i_gene->choice_ptr->lookup_res_ptr->lookupRot[i_gene->lookupRot_index]);

			while(dummychr.genes->seq_position!=var_pos[i].seq_position)
	    			dummychr.genes = dummychr.genes->nextgene;

			copyGENE(dummychr.genes, i_gene); 
			dummychr.genes = dummychr.firstgene;
			
			printf("\t\tafter:\n");
			
			lookup_energy_sasa_charge_other(i, i_gene, dummychr.firstgene, &sasa_other, &sasa_hphob_other, &charge_other, &E_transfer_other);
			CHROMOSOME_to_lookupEnergy(&dummychr, &parameters->fixedatoms_energy);

			printf("total_energy = %lf\n",dummychr.energy);


			while(dummychr.genes->seq_position!=var_pos[i].seq_position)
	    			dummychr.genes = dummychr.genes->nextgene;


			DEBUG=1;
		}
	}
*/


	fabs_dE = fabs(best_energy - current_energy);

	// printf("hqm %d %lf\n",i,fabs_dE);

		if(fabs_dE > EPS)	/* energy has  dropped */
		{	
	    		i_gene->choice_ptr = &var_pos[i].choice[best_res];
			i_gene->j_choice_index = best_res-1;
			i_gene->lookupRot_index = best_res_rot;
			i_gene->chi = i_gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[i_gene->lookupRot_index].chi;
			i_gene->lookupRot = &(i_gene->choice_ptr->lookup_res_ptr->lookupRot[i_gene->lookupRot_index]);

			dummychr2.genes = dummychr2.firstgene; 
	 		while(dummychr2.genes->seq_position!=i_gene->seq_position)
	 		{
	    			dummychr2.genes = dummychr2.genes->nextgene;
	 		}   
	 		copyGENE(dummychr2.genes, i_gene); 
			
			CHROMOSOME_to_lookupEnergy(&dummychr2,&parameters->fixedatoms_energy);

			/* has it actually dropped? */
			if(dummychr2.energy < dummychr.energy)
			{
			
				// printf("hqm drop %d %lf\n",i,dummychr2.energy - dummychr.energy);


				dummychr.energy = dummychr2.energy;
				copyGENE(dummychr.genes, i_gene); 
				dummychr.genes = dummychr.firstgene;
		
				neighborOptimizeFlag=0; GenX=0;


				/* reset the already_used_positions array since the energy has dropped */
				i=1;
				while(already_used_positions[i]!=ENDFLAG)
				{
					already_used_positions[i] = ENDFLAG;
					++i;
				}
			}
			else
			{
				dummychr2.genes = dummychr2.firstgene; 
	 			while(dummychr2.genes->seq_position!=i_gene->seq_position)
	    				dummychr2.genes = dummychr2.genes->nextgene;

				copyGENE(dummychr2.genes, dummychr.genes);
				dummychr2.energy = dummychr.energy;
			} 

		}	

		if(fabs_dE > EPS2)
				num_isoenergy_safety=0;
		else
		{	
			++num_isoenergy_safety;
			if(num_isoenergy_safety>safety_convergence)	/* bail out */
				GenX =  parameters->hqm_convergence + 20;
		}

	 if(neighborOptimizeFlag == 1)
		++GenX;

     }
     
	copyCHROMOSOME((*chr), dummychr);  
	CHROMOSOME_to_lookupEnergy(chr,&parameters->fixedatoms_energy);       

     free_memory(i_gene);
     
     free_CHROMOSOME(&dummychr);
	free_CHROMOSOME(&dummychr2);
}