/* EGAD: multistate_objective_function.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
Dec 12 2003

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

This file contains the multistate_objective_function ; called by multistate.cpp: CHROMOSOME_array_to_fitness.
Returns the energy or whatever score is being minimized during the optimization process.

These function may be replaced with custom functions, C++ objects etc as needed. However, the value must 
be returned to the multistate optimization via 
	double multistate_objective_function(CHROMOSOME *chr, ENERGY *energy, int num_competitor_structures) 

energy is an array of ENERGY; energy[0] is the ENERGY of the target structure, others are for structured
	competitors. These energies have been baseline subtracted (ala at all variable seq positions)

chr is the chromosome of interest. 
chr->energy_list is the ENERGY (w/ respect to the ala baseline) for the target state

num_competitor_structures is the number of structured non-target states being considered explicitly

MULTISTATE_OBJ_FUNCT_INPUTFILENAME is read by input_stuff.cpp: input_stuff, but is not actually opened there.
	The user may create functions that open and read this file if needed....see the example below


*/

#include "multistate_objective_function.h"

/*

The default here returns the energy difference between the target state and the boltzmann-average energy of the non-target 
	states, including the unfolded, aggregate, and specificity reference states.

*/

double target_non_target_ensemble_energy_gap(ENERGY *energy, int num_competitor_structures)
{
	double score;
	int i;
	extern int MULTISTATE_INDEX_STABILITY, MULTISTATE_INDEX_SPECIFICITY, MULTISTATE_INDEX_SOLUBILITY;
	extern double RT, OVERALL_ENERGY_SCALE;

	score = 0;
	
	// non-structured reference states; unf, specificity, solubility

	score += exp(-(OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_STABILITY].E_rss - energy[0].E_structure))/RT);
	score += exp(-(OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_SPECIFICITY].E_specificity - energy[0].E_structure))/RT);
	score += exp(-(OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_SOLUBILITY].E_solubility - energy[0].E_structure))/RT);
			
	// other structured states

	for(i=1;i<=num_competitor_structures;++i)
		score += exp(-(OVERALL_ENERGY_SCALE*(energy[i].E_structure - energy[0].E_structure))/RT);
			
	// G_non_target_ensemble = -RTlnQ , where Q = partition function (score to this point)
	//	dG_gap = G_target - Gnon_target_ensemble = RTlnQ (G_target is defined as 0 when calc'd Q)

	 /*
	printf("unf\t%lf\t%lf\t%lf\n",OVERALL_ENERGY_SCALE*energy[MULTISTATE_INDEX_STABILITY].E_rss, 
							  OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_STABILITY].E_rss - energy[0].E_structure),  
		exp(-(OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_STABILITY].E_rss - energy[0].E_structure))/RT));

	printf("specif\t%lf\t%lf\t%lf\n",OVERALL_ENERGY_SCALE*energy[MULTISTATE_INDEX_SPECIFICITY].E_specificity, 
							OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_SPECIFICITY].E_specificity - energy[0].E_structure),  
				exp(-(OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_SPECIFICITY].E_specificity - energy[0].E_structure))/RT));

	printf("solub\t%lf\t%lf\t%lf\n",OVERALL_ENERGY_SCALE*energy[MULTISTATE_INDEX_SOLUBILITY].E_solubility, 
							OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_SOLUBILITY].E_solubility - energy[0].E_structure),  
					exp(-(OVERALL_ENERGY_SCALE*(energy[MULTISTATE_INDEX_SOLUBILITY].E_solubility - energy[0].E_structure))/RT));

	printf("competitor\t%lf\t%lf\t%lf\n",OVERALL_ENERGY_SCALE*energy[1].E_structure,
									OVERALL_ENERGY_SCALE*(energy[1].E_structure - energy[0].E_structure),
									exp(-(OVERALL_ENERGY_SCALE*(energy[1].E_structure - energy[0].E_structure))/RT));

	printf("target\t%lf\t%lf\t%lf\n",OVERALL_ENERGY_SCALE*energy[0].E_structure,
								OVERALL_ENERGY_SCALE*(energy[0].E_structure-energy[0].E_structure),
								exp(-(OVERALL_ENERGY_SCALE*(energy[0].E_structure - energy[0].E_structure))/RT));

	printf("Q = %lf\n",score);
	 */
	 

	score = RT*log(score);

	 /*
	printf("RTlnQ = %lf\n",score);
	 */

	return(score);
}

/*

// 
	We are trying to design a binding protein w/ a target Kd
	We want to optimize the total energy of the target structure (including solubility, specificity), 
	but penalize with a harmonic function if the gap between the target and competitor (unbound) structures is 
	less than or greater than the target Kd

double binding_protein_objective_function(ENERGY *energy, int num_competitor_structures)
{
	double score, target_Kd;
	double dG;
	char *line=NULL;
	FILE *file_ptr;
	static double harmonic_constant, target_energy;
	static int first_time_flag = 0;
	extern double OVERALL_ENERGY_SCALE, RT;
	extern char *MULTISTATE_OBJ_FUNCT_INPUTFILENAME; // name is read by input_stuff

	if(first_time_flag == 0)
	{
		first_time_flag = 1;	

		// get values from some file - this function reads these values or whatever

		line = (char *)calloc(MAXLINE,sizeof(char));
	
		file_ptr = fopen_file(MULTISTATE_OBJ_FUNCT_INPUTFILENAME,"r");
		
		// not the most user-friendly format, but this is just an example
		
		fgets(line,MAXLINE,file_ptr);

		sscanf(line, "%lf %lf",&harmonic_constant, &target_Kd); 
	
		fclose(file_ptr);

		free_memory(line);

		target_energy = RT*log(target_Kd); 		
	}

	score = target_non_target_ensemble_energy_gap(energy, num_competitor_structures);

	dG = OVERALL_ENERGY_SCALE*(energy[0].E_structure - energy[1].E_structure);

	if(dG < 0)	// target structure has lower energy than competitor
		score += harmonic_constant*pow((dG - target_energy),2); // penalize if too tight or not tight enough
	else	// target structure has higher energy than competitor
		score += harmonic_constant*pow((dG + target_energy),2);


	return(score);
}

*/

#include "sequence_CHROMOSOME_stuff.h"

// the actual function called by multistate.cpp: CHROMOSOME_array_to_fitness
double multistate_objective_function(CHROMOSOME *chr, ENERGY *energy, int num_competitor_structures)
{
	// return(your_favorite_function_here());
        // for example, for the binding protein function sketched above,
        //      return(binding_protein_objective_function(energy, num_competitor_structures));

	static int first_time_flag=0;
	static ENERGY starting_energy;
	static double starting_pseudo_dG;
	double penalty;
	extern double OVERALL_ENERGY_SCALE, ENERGY_ERROR;
	extern double VDW_CUTOFF;
	extern int MULTISTATE_INDEX_STABILITY;

	energy[0].dG_gap = 	target_non_target_ensemble_energy_gap(energy, num_competitor_structures);
	chr->energy_list.dG_gap = energy[0].dG_gap;

	energy[0].pseudo_dG = energy[MULTISTATE_INDEX_STABILITY].pseudo_dG;
	chr->energy_list.pseudo_dG = energy[MULTISTATE_INDEX_STABILITY].pseudo_dG;
	
	// initialize some variables 
	if(first_time_flag == 0)
	{
		starting_energy = energy[0];
		starting_pseudo_dG = energy[MULTISTATE_INDEX_STABILITY].pseudo_dG;
		first_time_flag = 1;
	}
	
	penalty = 0;

	if(OVERALL_ENERGY_SCALE*energy[0].E_total > OVERALL_ENERGY_SCALE*starting_energy.E_total + 2*ENERGY_ERROR)
		penalty += VDW_CUTOFF*OVERALL_ENERGY_SCALE*(energy[0].E_total - starting_energy.E_total);

	if(energy[MULTISTATE_INDEX_STABILITY].pseudo_dG > starting_pseudo_dG + 2*ENERGY_ERROR)
		penalty += VDW_CUTOFF*(energy[MULTISTATE_INDEX_STABILITY].pseudo_dG - starting_pseudo_dG);

	
	// printf("gap  = %lf\t penalty = %lf\t total = %lf\n",energy[0].dG_gap , penalty,energy[0].dG_gap + penalty);

	return(energy[0].dG_gap + penalty);
}