/* EGAD: dee_utilities.cpp
   
Mark Voorhies, 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 functions for interfacing DEE with the EGAD datastructures.
   
*/


#include "dee_utilities.h"

extern float NON_INTERACT_PTR, FIXED_POSITION_PTR;
extern LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X;
extern LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X;

/* Initializes the lookuptable for dee */
void initialize_lookuptable_for_dee(PROTEIN *protein)
{
  int i, i_res, i_res_rot;
  int q;
  int j, j_res;
  int dee_index;
  short int resimer_index;
  GENE i_gene, j_gene;
  extern char *LOOKUP_TABLE_DIRECTORY;

  strcpy(LOOKUP_TABLE_DIRECTORY,protein->parameters.lookup_energy_table_directory);	

  i_gene = (GENE)malloc(sizeof(MENDEL));
  j_gene = (GENE)malloc(sizeof(MENDEL));
    
  i=1; dee_index=0;
  while(protein->var_pos[i].seq_position!=ENDFLAG)
    {
      protein->var_pos[i].number_of_resimers = 0;
      for(i_res=1;i_res<=protein->var_pos[i].number_of_choices; ++i_res)
	protein->var_pos[i].number_of_resimers += protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers;
      
      if(protein->var_pos[i].number_of_resimers > 1)
      {
	protein->var_pos[i].dee_index = dee_index;
	++dee_index;
      }

      if(protein->var_pos[i].resimer == NULL)
	protein->var_pos[i].resimer = (RESIMER *)calloc((protein->var_pos[i].number_of_resimers+2), sizeof(RESIMER));
      
      q=1; resimer_index=0;
      for(i_res=1;i_res<=protein->var_pos[i].number_of_choices;++i_res)
	for(i_res_rot=1;i_res_rot<=protein->var_pos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers;++i_res_rot)
	  {
	    protein->var_pos[i].resimer[q].i_res = i_res;
	    protein->var_pos[i].resimer[q].i_res_rot = i_res_rot;
	    ++q;
	  
            protein->ERESROT.resimer_index = resimer_index;
	    ++resimer_index;

		j=i+1;
		while(protein->var_pos[j].seq_position!=ENDFLAG)
		{
			for(j_res=1;j_res<=protein->var_pos[j].number_of_choices;++j_res)
			{
				 if(protein->ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X)
						if(protein->ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX == NULL)
							if(protein->var_pos[i].choice[i_res].in_use_flag==1 && protein->var_pos[j].choice[j_res].in_use_flag==1)
							{
								i_gene->choice_ptr = &protein->var_pos[i].choice[i_res];
      								i_gene->seq_position = protein->var_pos[i].seq_position;
      								i_gene->j_choice_index = i_res-1;
								i_gene->lookupRot = &i_gene->choice_ptr->lookup_res_ptr->lookupRot[1];
	  							i_gene->lookupRot_index = 1;
      								i_gene->varpos_ptr = &protein->var_pos[i];

      								j_gene->choice_ptr = &protein->var_pos[j].choice[j_res];
      								j_gene->seq_position = protein->var_pos[j].seq_position;
      								j_gene->j_choice_index = j_res-1;
								j_gene->lookupRot = &j_gene->choice_ptr->lookup_res_ptr->lookupRot[1];
	      							j_gene->lookupRot_index = 1;
								j_gene->varpos_ptr = &protein->var_pos[j];

      								i_gene->nextgene = NULL; j_gene->nextgene = NULL; 
	
								load_calc_save_lookupResX(i_gene, j_gene);
							}
			}
			++j;
		}

	  }
      ++i;
    }

	free_memory(i_gene); free_memory(j_gene);

}

/* returns the pair energy for resimer pair (i,i_resimer),(j,j_resimer); i,j = seq indexes.
		if i==j, and i_resimer == j_resimer, self energy for i,i_resimer returned
		else, return 0
*/ 
double get_energy(const PROTEIN *protein, const int i, const int i_resimer, 
		  const int j, const int j_resimer)
{
	int i_minus_j,j_minus_i;
        

 	if(j == i)
  	{
    		if(j_resimer == i_resimer)
		{
      			return get_self_energy(protein, i, i_resimer);
  		}
    		else
   		{
      			// By definition, different resimers at the same position can not see
      			// each other.
      			return 0.0;
    		}
  	}

	/* one or both residue choices are not in use */
  if(protein->var_pos[i].choice[RESIMER_i_res].in_use_flag==0 || protein->var_pos[j].choice[RESIMER_j_res].in_use_flag==0)
  {
	if(protein->var_pos[i].choice[RESIMER_i_res].in_use_flag==0 && protein->var_pos[j].choice[RESIMER_j_res].in_use_flag==0)
		return(2e5);	/* neither are in use */
	else
		return(1e5);
  }


  if(j > i)
    {
	j_minus_i = j-i;
      if(ERESROT_resimer.lookupX[j_minus_i].lookupResX != NON_INTERACT_LOOKUP_RES_X)
	{
	  if(ERESROT_resimer.lookupX[j_minus_i].lookupResX[RESIMER_j_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X)
	    {
	      if(ERESROT_resimer.lookupX[j_minus_i].lookupResX[RESIMER_j_res-1].lookupRotX[RESIMER_j_res_rot-1].energy_var_var != NON_INTERACT_PTR)
		{
		  return(ERESROT_resimer.lookupX[j_minus_i].lookupResX[RESIMER_j_res-1].lookupRotX[RESIMER_j_res_rot-1].energy_var_var);
		}
	    }
	} 
    }
  else // j < i
    {
	i_minus_j = i-j;
      if(ERESROT_resimer_j.lookupX[i_minus_j].lookupResX != NON_INTERACT_LOOKUP_RES_X)
	{
	  if(ERESROT_resimer_j.lookupX[i_minus_j].lookupResX[RESIMER_i_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X)
	    {
	      if(ERESROT_resimer_j.lookupX[i_minus_j].lookupResX[RESIMER_i_res-1].lookupRotX[RESIMER_i_res_rot-1].energy_var_var != NON_INTERACT_PTR)
		{
		  return(ERESROT_resimer_j.lookupX[i_minus_j].lookupResX[RESIMER_i_res-1].lookupRotX[RESIMER_i_res_rot-1].energy_var_var);
		}
	    }
	}
    }
    
  return(0);
}


/* returns internal and side-bkbn energy (self energy) for resimer i_resimer at seq index i */
double get_self_energy(const PROTEIN *protein, 
		       const int i, const int i_resimer)
{
  if(ERESROT_resimer.energy_var_fix!=FIXED_POSITION_PTR)
  {	
  	if(protein->var_pos[i].choice[RESIMER_i_res].in_use_flag==0)
		return(ERESROT_resimer.energy_var_fix+1e10);
	else	
		return(ERESROT_resimer.energy_var_fix);
  }
  else
	return(0.0);
}

/* for var_pos at some position, and an input_resimer at that position, find i_res and i_res_rot;
	used for converting from resimers back to residuetypes/rotamers 
*/
void get_i_res_i_res_rot_for_resimer(VARIABLE_POSITION *var_pos, int input_resimer, int *input_i_res, int *input_i_res_rot)
{
	int resimer, i_res, i_res_rot;
	
	
	i_res = 1; resimer= 0;
	while(i_res<=var_pos->number_of_choices)
	{
		i_res_rot = 1;
		while(i_res_rot<=var_pos->choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers)
		{
			++resimer;

			if(resimer > var_pos->number_of_resimers)
			{
				fprintf(stderr,"ERROR (1) Cannot find i_res, i_res_rot for resimer %d at seq_pos %d\n",	
						input_resimer,var_pos->seq_position);
				fprintf(stderr,"dumping core.....\n");
				abort();
			}
			
			if(resimer == input_resimer)
			{
				*input_i_res = i_res;
				*input_i_res_rot = i_res_rot;

				return;
			}

			++i_res_rot; 
		}

		++i_res;
	}

	fprintf(stderr,"ERROR (2) Cannot find i_res, i_res_rot for resimer %d at seq_pos %d\n",	input_resimer,var_pos->seq_position);
	fprintf(stderr,"dumping core.....\n");
	abort();

}


/* write logfile_line to logfile */
void dee_logfile(char logfile_line[], const char *logfile)
{
	FILE *log_file_ptr;
	extern int LOGFILE_FLAG;
    	time_t now;

	if(LOGFILE_FLAG != 0)
  	{
		if(logfile!=NULL)
		{
			now = time(NULL);
			log_file_ptr = fopen_file(logfile,"a");	   
  			fprintf(log_file_ptr,"dee\t%s  %s",logfile_line, ctime(&now));
			fclose(log_file_ptr);
		}
		else
		{	
			fprintf(stderr,"ERROR logfile undefined for LOGFILE_FLAG 1\n");
			fprintf(stderr,"dumping core.....\n");
			abort();
		}
  	}	
}