/* EGAD: general_GA.h
   
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 a general-purpose genetic algorithm with continuous moves

*/

#include "general_GA.h"

struct proto_plasmidGENE
{
        double transform_element;
	double upper_bound;
	double lower_bound;
        struct proto_plasmidGENE *next;
};

typedef struct proto_plasmidGENE plasmidMENDEL;        
typedef plasmidMENDEL *plasmidGENE;                    

typedef struct
{
        plasmidGENE genes;
        plasmidGENE firstgene;
        double energy;
        double mating_freq;
	int n;
} PLASMID;

#define copyPLASMID(b,a)        \
{       \
        (b).genes = (b).firstgene;   \
        (a).genes = (a).firstgene;   \
        while((b).genes->next != NULL)     \
        {       \
                (b).genes->transform_element = (a).genes->transform_element;  \
                (b).genes = (b).genes->next;  \
                (a).genes = (a).genes->next;  \
        }       \
        (b).genes = (b).firstgene;   \
        (a).genes = (a).firstgene;   \
	(b).energy = (a).energy;	\
}

#define mutate_plasmidGENE(gene)	\
{						\
        (gene)->transform_element = rand_double_window((gene)->lower_bound,(gene)->upper_bound);		\
}

void mutagenize_PLASMID(PLASMID *chr, double mut_freq)
{
        chr->genes = chr->firstgene;
        
        while(chr->genes->next!=NULL)
        {
        	if(dice()<=mut_freq)
			mutate_plasmidGENE(chr->genes);

                chr->genes = chr->genes->next;
        }

        chr->genes = chr->firstgene;
}

void randomize_PLASMID(PLASMID *chr)
{
        chr->genes = chr->firstgene;
        
        while(chr->genes->next!=NULL)
        {

        	mutate_plasmidGENE(chr->genes);

                chr->genes = chr->genes->next;
        }

        chr->genes = chr->firstgene;
}

void initialize_PLASMID_array(PLASMID *chr, int n, double *initial_soln, double *lower_bound, double *upper_bound, int popsize)
{
	int i, j;

	
		for(i=1;i<=popsize;++i)
		{
			chr[i].n = n;

			chr[i].genes = (plasmidMENDEL *)malloc(sizeof(plasmidMENDEL));
			chr[i].firstgene = chr[i].genes;

			chr[i].genes->next = NULL;

		
			for(j=1;j<=n;++j)
			{
				chr[i].genes->transform_element = initial_soln[j];
				chr[i].genes->lower_bound = lower_bound[j];
				chr[i].genes->upper_bound = upper_bound[j];


				mutate_plasmidGENE(chr[i].genes);		

				chr[i].genes->next = (plasmidMENDEL *)malloc(sizeof(plasmidMENDEL));
				chr[i].genes = chr[i].genes->next;
				chr[i].genes->next = NULL;

			}
			chr[i].genes = chr[i].firstgene;
		}
	
}


inline void score_PLASMID(PLASMID *plasmid, double (* func)(double [], POWELL *), POWELL *powell)
{
	static double *soln=NULL;
	static int n_max=0;
	int i, logfile_flag;
	time_t start_time;
	extern int LOGFILE_FLAG;

	if(plasmid->n > n_max)
	{	
		if(soln!=NULL)
			free_memory(soln);

		n_max = plasmid->n;	

		soln = (double *)calloc(n_max+2,sizeof(double));
	}

	i=1;
	plasmid->genes = plasmid->firstgene;
	while(plasmid->genes->next!=NULL)
	{
		soln[i] = plasmid->genes->transform_element;
		++i; plasmid->genes = plasmid->genes->next;
	}
	plasmid->genes = plasmid->firstgene;
	
	plasmid->energy = (*func)(soln,powell);

	
		start_time =time(NULL); 
		logfile_flag = LOGFILE_FLAG;
		LOGFILE_FLAG = 0;
		powell_minimization(soln, plasmid->n, plasmid->n, &plasmid->energy, func, powell, start_time);
		LOGFILE_FLAG = logfile_flag;
	 

	i=1;
	plasmid->genes = plasmid->firstgene;
	while(plasmid->genes->next!=NULL)
	{
		plasmid->genes->transform_element = soln[i];
		++i; plasmid->genes = plasmid->genes->next;
	}
	plasmid->genes = plasmid->firstgene;

}

void mating_freq(PLASMID *plasmid, int popsize)
{
        int i, half_popsize;
        double under;
        
	half_popsize = popsize/2;
	
        
        
	under=0;
        for(i=1;i<=half_popsize;++i)
        {
                plasmid[i].mating_freq = 1.0/(double)i;
                under += plasmid[i].mating_freq; 
        }
        
        for(i=1;i<=half_popsize;++i)
                plasmid[i].mating_freq = plasmid[i].mating_freq/under;
}

void sort_PLASMID(PLASMID *plasmid, int *first, int *last)
{
  int left, right;
  PLASMID pivot, dummy;
  

  left = *first;
  right = *last;
  
  pivot = plasmid[(left+right)/2];
  
  while(left <= right)
    {
      while(plasmid[right].energy > pivot.energy)
        --right;
      while(plasmid[left].energy < pivot.energy)
        ++left;
      if(left <= right)
        {
          dummy = plasmid[left];
          plasmid[left] = plasmid[right];
          plasmid[right] = dummy;
          ++left;
          --right;
        }
    }
  if(*first < right)
    sort_PLASMID(plasmid, first,&right);
  if(*last > left)
    sort_PLASMID(plasmid, &left,last);
}


void mate(PLASMID *plasmid, int popsize, double mut_freq)
{
        int i;
        static double *mate_freq=NULL;
        int A,B;  /* indexes for parents */
        int a,b;  /* indexes for children */
        int newpopsize, half_popsize;
        plasmidGENE dummygene;

	half_popsize = popsize/2;

        if(mate_freq == NULL)
        {
                mate_freq = (double *)calloc(half_popsize+4,sizeof(double));
        }

        mate_freq[0] = 0; mate_freq[1] = 0; mate_freq[2] = 0; 
        for(i=1;i<=half_popsize;++i)  
                mate_freq[i] = plasmid[i].mating_freq + mate_freq[i-1];

        newpopsize = half_popsize;
        while(newpopsize <= popsize-2)
        {
                A=1; B=1;
                while(A==B) /* while(plasmid[A].score == plasmid[B].score) */
                {
                        A = roll_loaded_dice(mate_freq, half_popsize);     /* selecting parent A */
                        B = roll_loaded_dice(mate_freq, half_popsize);     /* selecting parent B */
                }
                ++newpopsize; a=newpopsize; 
                ++newpopsize; b=newpopsize;

                copyPLASMID(plasmid[a],plasmid[A]);
                copyPLASMID(plasmid[b],plasmid[B]);

                plasmid[a].genes = plasmid[a].firstgene;
                plasmid[b].genes = plasmid[b].firstgene;
                i=1;
                while(plasmid[a].genes->next!=NULL)
                {
                        if(dice()<=0.5)
                        {
                                dummygene = plasmid[a].genes->next;
                                plasmid[a].genes->next = plasmid[b].genes->next;
                                plasmid[b].genes->next = dummygene;
                        }

                        if(dice()<=mut_freq)
                                mutate_plasmidGENE(plasmid[a].genes);

                        if(dice()<=mut_freq)
                                mutate_plasmidGENE(plasmid[b].genes);
                
                        plasmid[a].genes = plasmid[a].genes->next;
                        plasmid[b].genes = plasmid[b].genes->next;

                        ++i;
                }
                plasmid[a].genes = plasmid[a].firstgene;
                plasmid[b].genes = plasmid[b].firstgene;
   
        }

}

/* 
        if chr[i] is isoenergetic with any chr[1]->chr[i-1], returns 1; else 0;
*/
int is_plasmid_i_isoenergetic_with_others(int i, PLASMID *chr, int popsize)
{
        int n;

        for(n=1;n<=popsize;++n)
                if(i!=n)
                        if(fabs(chr[i].energy - chr[n].energy)<=EPS)
                                return(1);

        return(0);
}

void get_rid_of_isoenergetic_plasmids(PLASMID *chr, double mut_freq, int pop_size, double (*func)(double [], POWELL *), POWELL *powell )
{
        int i;
        int flag;
        int ctr, local_ctr;

        flag=0; ctr=0;
        while(flag==0 && ctr<=pop_size)
        {
                flag=1;
                for(i=1;i<=pop_size;++i)
                {
                        local_ctr=0;
                        while(is_plasmid_i_isoenergetic_with_others(i, chr, pop_size)==1 && local_ctr<=pop_size)
                        {
                                flag=0;
                                mutagenize_PLASMID(&chr[i], local_ctr*mut_freq);
                                score_PLASMID(&chr[i], func, powell);
                                ++local_ctr;
                        }       
                }
                ++ctr;
        }
}

void general_purpose_GA(double *soln, double *lower_bound, double *upper_bound, int n, int max_iterations, double *score, double (*func)(double [], POWELL *), POWELL *powell, time_t start_time)
{
	extern double MAX_OPTIMIZATION_TIME;
	extern int LOGFILE_FLAG, NUM_FUNCTION_CALLS, GET_PID;
 	extern char *CURRENT_WORKING_DIRECTORY;
	static char *filename=NULL, *escape_hatch_filename=NULL;
	static int n_max=0;
	FILE *file_ptr;
	static PLASMID *plasmid=NULL;
	time_t now;
	int popsize,j;
	int half_popsize,first,i,iter,cycle;
	double E_previous;


	popsize = powell->protein->parameters.pop_size;

	half_popsize = popsize/2;
	first=1;

	if(escape_hatch_filename==NULL)
	{
		escape_hatch_filename = (char *)calloc(MAXLINE,sizeof(char));
		sprintf(escape_hatch_filename,"%s/escape.gp_GA.%d",CURRENT_WORKING_DIRECTORY,GET_PID);
	}

	if(n > n_max)
	{
		if(plasmid!=NULL)
			free_memory(plasmid);
		
		plasmid = (PLASMID *)calloc(popsize+2,sizeof(PLASMID));
		initialize_PLASMID_array(plasmid, n, soln, lower_bound, upper_bound, popsize);

		n_max=n;
	}


	

	now = time(NULL);
	if(LOGFILE_FLAG != 0)
		{
			if(filename==NULL)
				filename = (char *)calloc(MAXLINE,sizeof(char));
			sprintf(filename, "%s.log", powell->protein->parameters.output_prefix);
	    		file_ptr = fopen_file(filename, "a");
	    		fprintf(file_ptr, "gp_ga start n = %d\t\t", n);
			fprintf(file_ptr, "%f\t%s",  (*score), ctime(&now) );
			fprintf(file_ptr,"To exit ga minimization and move to the next step:");
			fprintf(file_ptr,"\ttouch %s\n",escape_hatch_filename);
	    		fclose(file_ptr);
		}

	for(i=1;i<=popsize;++i)
	{
		plasmid[i].genes = plasmid[i].firstgene;
		j=1;
       	 	while(plasmid[i].genes->next!=NULL)
        	{
			plasmid[i].genes->transform_element = soln[j];
			plasmid[i].genes->lower_bound = lower_bound[j];
			plasmid[i].genes->upper_bound = upper_bound[j];

                	plasmid[i].genes = plasmid[i].genes->next;
			++j;
        	}
        	plasmid[i].genes = plasmid[i].firstgene;

		if(i>1)
			randomize_PLASMID(&plasmid[i]);

		score_PLASMID(&plasmid[i], func, powell);

		  
		if(LOGFILE_FLAG != 0)
		{
			sprintf(filename, "%s.log", powell->protein->parameters.output_prefix);
	    		file_ptr = fopen_file(filename, "a");
			now = time(NULL);
	    		fprintf(file_ptr, "initialize %d/%d\t%lf\t%s",i,popsize,plasmid[i].energy,ctime(&now));
			fclose(file_ptr);
		}
	}

	get_rid_of_isoenergetic_plasmids(plasmid, powell->protein->parameters.mutation_freq, popsize, func, powell);

	sort_PLASMID(plasmid, &first, &popsize);
	

	iter=0; E_previous=1e10; cycle=0;
	while(iter<=max_iterations)
	{
		++cycle;

		mating_freq(plasmid, half_popsize);
		mate(plasmid, popsize, powell->protein->parameters.mutation_freq);


		for(i=half_popsize;i<=popsize;++i)
		{
			score_PLASMID(&plasmid[i], func, powell);
			
		 	 /*
			if(LOGFILE_FLAG != 0)
			{
				sprintf(filename, "%s.log", powell->protein->parameters.output_prefix);
	    			file_ptr = fopen_file(filename, "a");
				now = time(NULL);
	    			fprintf(file_ptr, "%d\t%lf\t%s",i,plasmid[i].energy,ctime(&now));
				fclose(file_ptr);
			}
		 		*/	 
			

		}

		get_rid_of_isoenergetic_plasmids(plasmid, powell->protein->parameters.mutation_freq, popsize, func, powell);
		sort_PLASMID(plasmid, &first, &popsize);
		

		if(LOGFILE_FLAG != 0)
		{
			sprintf(filename, "%s.log", powell->protein->parameters.output_prefix);
	    		file_ptr = fopen_file(filename, "a");
			now = time(NULL);
	    		fprintf(file_ptr, "ga iteration %d\t%d\t%d\t%f\t%s", cycle,iter, NUM_FUNCTION_CALLS, plasmid[1].energy,
											ctime(&now) );

	    		fclose(file_ptr);

		
		}

		if(NUM_FUNCTION_CALLS >= powell->protein->parameters.max_function_calls || 
				difftime(now,start_time) >= MAX_OPTIMIZATION_TIME ||
				does_this_file_exist(escape_hatch_filename)==1)
		{
			iter = max_iterations + 2;

			if(does_this_file_exist(escape_hatch_filename)==1)
			{
				NUM_FUNCTION_CALLS = powell->protein->parameters.max_function_calls+10;
				rm_file(escape_hatch_filename);
			}

			i=1;
			plasmid[1].genes = plasmid[1].firstgene;
			while(plasmid[1].genes->next!=NULL)
			{
				soln[i] = plasmid[1].genes->transform_element;
				++i; plasmid[1].genes = plasmid[1].genes->next;
			}

			*score = plasmid[1].energy;

			if(LOGFILE_FLAG != 0)
			{
			sprintf(filename, "%s.log", powell->protein->parameters.output_prefix);
	    		file_ptr = fopen_file(filename, "a");
			now = time(NULL);
	    		fprintf(file_ptr, "ga done final score = %lf\t%s", plasmid[1].energy, ctime(&now) );

	    		fclose(file_ptr);
			}


			return;
		}

		if(fabs(E_previous - plasmid[1].energy)<EPS)
			++iter;
		if(plasmid[1].energy < E_previous)
		{
			iter = 0;
			E_previous = plasmid[1].energy;
		}	

		if(fabs(plasmid[1].energy - plasmid[popsize].energy)<EPS)
			iter = max_iterations + 10;

	}

	i=1;
			plasmid[1].genes = plasmid[1].firstgene;
			while(plasmid[1].genes->next!=NULL)
			{
				soln[i] = plasmid[1].genes->transform_element;
				++i; plasmid[1].genes = plasmid[1].genes->next;
			}


	*score = plasmid[1].energy;

	if(LOGFILE_FLAG != 0)
			{
			sprintf(filename, "%s.log", powell->protein->parameters.output_prefix);
	    		file_ptr = fopen_file(filename, "a");
			now = time(NULL);
	    		fprintf(file_ptr, "ga done final score = %lf\t%s", plasmid[1].energy, ctime(&now) );

	    		fclose(file_ptr);
			}

	return;
}