/*
schedule_analysis

Author: Darcy Cook


*/

#define PHI_ERROR_MEASURE (0.0001)
#define CONVERGE_ERROR (1E-4)
#define PI_ERROR (1E-4)
#define EMPTY_SCHEDULE_SLOT (-1)
#define ANALYZED_RATIO_MAX	(0.999)
#define INIT_VAL (0.65)
#define MAX_ITERATION (400)
#define MAX_MU (5)
#define MIN_MU (0.5)
#define MAX_NUM_MOVES (30)
#define MIN_NUM_MOVES (1)
#define K_const (2)
#define RESTART_THRESHOLD (5)
#define WALK_FILE_NAME "walk_CSA3.txt"
#define RESULT_FILE_NAME "result_CSA3.txt"
#define N_const (0.9)

#include <stdio.h>
#include <string.h>
#include <sys/time.h>
#include <stdlib.h>
#include <mpi.h>
#include <math.h>
#include "matrix.h"

/*******************
create_service_matrix

Description: This function creates the service process matrices. The memory for the created matricies is assumed to have been
			 previously allocated.
			 
Parameters: S_ptr - this is a pointer to the service matrix
			S0_ptr - this is a pointer to the end of the process vector
			Beta_ptr - this is a pointer to the beginning of the process vector
			service_prob - this is the probability of moving to the next phase of the service matrix
*******************/
void create_service_matrix(Matrix *S_ptr, Matrix *S0_ptr, Matrix *Beta_ptr, double service_prob)
{
	int i;
	
	//first set all the elements in the three matricies to 0
	all_zeros(S_ptr);
	all_zeros(S0_ptr);
	all_zeros(Beta_ptr);
	
	//set the first element of Beta to 1
	get_matrix_element(Beta_ptr, 0, 0) = 1;
	
	//set the last element in the S0 vector to the service probability
	get_matrix_element(S0_ptr, S0_ptr->rows-1, 0) = service_prob;
	
	for (i = 0; i < S_ptr->rows; i++)
	{
		//set the diagonal to the probability of staying in the same state (1 - service probability)
		get_matrix_element(S_ptr, i, i) = 1 - service_prob;
		if (i != (S_ptr->rows-1))
		{
			//set the diagonal row above the main diagonal to the probability of progressing to the next state
			get_matrix_element(S_ptr, i, i+1) = service_prob;
		}
	}
}

/*******************
create_vacation_process

Description: This function creates the vacation process matrices for a specified processor. It is assumed that the memory
			 for all of the matrices used has been allocated previous to this function.
			 
Parameters: Vacation_ptr - this is a pointer to the matrix containing the entire vacation processes
			S_ptr - this is a pointer to the service matrix
			S0_ptr - this is a pointer to the end of the process vector
			Beta_ptr - this is a pointer to the beginning of the process vector
			S0_mult_Beta_ptr - this is a pointer to a matrix which is S0*Beta
			Alpha_ptr - this is a pointer to a vector containing the probablity of memory arrivals for each processor
			Phi_ptr - this is a pointer to a vector containing the probability of a memory arrival within the vacation time for each processor
			num_procs - this is the number of total processors in the system
			curr_proc - this is the index of the processor that the vacation process is currently being built for
			service_length - this is the length of the service process
			vacation_order - this is the length of the vacation process
			
Returns: 1 if the vacation process was built successfully, and 0 otherwise
*******************/
int create_vacation_process(Matrix* Vacation_ptr, Matrix *S_ptr, Matrix *S0_ptr, 
							 Matrix *Beta_ptr, Matrix *S0_mult_Beta_ptr, Matrix *Alpha_ptr, Matrix *Phi_ptr, 
							 int num_procs, int curr_proc, int service_length, int vacation_order, int print_out)
{
	int curr_row_index, curr_col_index, i, j;
	double temp;
		
	//first set the Vacation process to all zeros
	all_zeros(Vacation_ptr);
	
	/* first create the starting vector of the vacation process	
	Here is the equivalent Matlab code for what is being implemented (4 processor case is shown)
	
	v_st = [ (1-pa(mod(i,N)+1)) pa(mod(i,N)+1)*B 0 zeros(1,K) 0 zeros(1,K) ];
	*/
	get_matrix_element(Vacation_ptr, 0, 1) = 1 - get_matrix_element(Phi_ptr, 0, ((curr_proc+1) % num_procs));	
	scalar_matrix_block_mult(get_matrix_element(Phi_ptr, 0, ((curr_proc+1) % num_procs)), Beta_ptr, Beta_ptr->rows, Beta_ptr->cols, 0, 0, Vacation_ptr, 0, 2);	
	
	if (print_out == 1)
	{
		printf("v_st = \n");
		print_matrix_block(Vacation_ptr, 0, 1, 1, vacation_order);
	}
	
	/* next the vector representing the end of vacation (absorbtion) is created
	here is the Matlab equivalent for the 4 processor case
	
	V0 = [ 0;         
           zeros(K,1);
           0;         
           zeros(K,1);
           1-a(mod(i+2,N)+1);  
           zeros(K,1)];
	*/
	//this is the index that will be set
	curr_row_index = 1 +(num_procs-2)*(service_length+1);
	get_matrix_element(Vacation_ptr, curr_row_index, 0) = 1 - get_matrix_element(Alpha_ptr, 0, ((curr_proc+num_procs-1) % num_procs));
	
	if (print_out == 1)
	{
		printf("V0 = \n");
		print_matrix_block(Vacation_ptr, 1, 0, vacation_order, 1);
	}
	
	/* next the vacation service portion of the matrix is created
	  here is the equivalent Matlab code
	/*V = [ 0                     a(mod(i,N)+1)*B         (1-pa(mod(i+1,N)+1))*(1-a(mod(i,N)+1))   pa(mod(i+1,N)+1)*(1-a(mod(i,N)+1))*B 0                                   zeros(1,K);       
           (1-a(mod(i,N)+1))*S0  a(mod(i,N)+1)*(S0*B)+S  zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K);       
           0                     zeros(1,K)              0                                   a(mod(i+1,N)+1)*B                       (1-pa(mod(i+2,N)+1))*(1-a(mod(i+1,N)+1)) pa(mod(i+2,N)+1)*(1-a(mod(i+1,N)+1))*B;      
           zeros(K,1)            zeros(K,K)              (1-a(mod(i+1,N)+1))*S0              a(mod(i+1,N)+1)*(S0*B)+S                zeros(K,1)                          zeros(K,K);       
           0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   a(mod(i+2,N)+1)*B;    
           zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              (1-a(mod(i+2,N)+1))*S0              a(mod(i+2,N)+1)*(S0*B)+S    ];

     %V is the matrix that represents the vacation before it is done
        V = [ 0                     a(mod(i,N)+1)*B         (1-pa(mod(i+1,N)+1))*(1-a(mod(i,N)+1))   pa(mod(i+1,N)+1)*(1-a(mod(i,N)+1))*B ;       
              (1-a(mod(i,N)+1))*S0  a(mod(i,N)+1)*(S0*B)+S  zeros(K,1)                          zeros(K,K)                             ;       
              0                     zeros(1,K)              0                                   a(mod(i+1,N)+1)*B                     ;      
              zeros(K,1)            zeros(K,K)              (1-a(mod(i+1,N)+1))*S0              a(mod(i+1,N)+1)*(S0*B)+S              ];      
                     
               */
    
     curr_row_index = 1;
     curr_col_index = 1;
     for (i = 0; i < (num_procs-1); i++)  
     {
	     /*each row consists of a matricies [ A B ] starting at the curr_row_index and the curr_col_index, except the last row which has just the A matrix */
	     
	     /*first create the matrix  A =[ 0    				alpha*beta
	                                    (1-alpha)*S0		alpha*(S0*beta)+S ] */
	 	//this is alpha*beta
	    scalar_matrix_block_mult(get_matrix_element(Alpha_ptr, 0, ((curr_proc+i+1) % num_procs)), Beta_ptr, Beta_ptr->rows, Beta_ptr->cols, 0, 0, Vacation_ptr, curr_row_index, curr_col_index+1);
	 	//this is (1-alpha)*S0
	    scalar_matrix_block_mult((1-get_matrix_element(Alpha_ptr, 0, ((curr_proc+i+1) % num_procs))), S0_ptr, S0_ptr->rows, S0_ptr->cols, 0, 0, Vacation_ptr, curr_row_index+1, curr_col_index);					   
	 	//this is alpha*(S0+beta)
	    scalar_matrix_block_mult(get_matrix_element(Alpha_ptr, 0, ((curr_proc+i+1) % num_procs)), S0_mult_Beta_ptr, S0_mult_Beta_ptr->rows, S0_mult_Beta_ptr->cols, 0, 0, Vacation_ptr, curr_row_index+1, curr_col_index+1);
	 	//this is adding S to the result of alpha*(S0+beta), which is alpha*(S0+beta)+S
	    matrix_block_addition(S_ptr->rows, S_ptr->cols, 
	 						Vacation_ptr, curr_row_index+1, curr_col_index+1, 
	 						S_ptr, 0, 0,
						 	Vacation_ptr, curr_row_index+1, curr_col_index+1);
		
		//the last row does not contain the B matrix
		if (i < (num_procs-2))
		{
			/*now create the matrix  B =[ (1-phi)*(1-alpha)		phi*(1-alpha)*beta
	        	                          0's				 	0's					 ] */
			//this is (1-phi)*(1-alpha)
	    	get_matrix_element(Vacation_ptr, curr_row_index, curr_col_index+1+service_length) = 
	    				(1 - get_matrix_element(Phi_ptr, 0, ((curr_proc+i+2) % num_procs))) * (1 - get_matrix_element(Alpha_ptr, 0, ((curr_proc+i+1) % num_procs)));
			//this is phi*(1-alpha)
			temp = get_matrix_element(Phi_ptr, 0, ((curr_proc+i+2) % num_procs)) * (1 - get_matrix_element(Alpha_ptr, 0, ((curr_proc+i+1) % num_procs)));	   
	    	//now multiply temp by beta to create phi*(1-alpha)*beta
	    	scalar_matrix_block_mult(temp, Beta_ptr, Beta_ptr->rows, Beta_ptr->cols, 0, 0, Vacation_ptr, curr_row_index, curr_col_index+2+service_length);
		}
		
	    //now increase the row and column indicies for the next row
	    curr_row_index += 1 + service_length;
	    curr_col_index += 1 + service_length;	    
     }
     
    if (print_out == 1)
	{
		printf("V = \n");
		print_matrix_block(Vacation_ptr, 1, 1, vacation_order, vacation_order);
	}
                  
    return 1;          
}

/*******************
calc_phi

Description: This function calculates the value of phi for a given alpha and vacation process.
			 
Parameters: Vacation_ptr - this is a pointer to the matrix containing the entire vacation process
			curr_alpha - this is the value of alpha that the phi calculation is based on.
			full_vacation_order - this is the length of the vacation process
			
Returns: the value of phi calculated
*******************/
double calc_phi(Matrix* Vacation_ptr, double curr_alpha, int full_vacation_order)
{
	/* this is the cumlative of the probability distribution of the vacation ending */
	double sum = 0;
	double alpha_cdf; //this is the cumlative distribution of alpha at time n
	double total;
	double phi;
	int n, i, j, m;	
	
	/* this is matrix that contains the probabilities of first arriving at state i from 
	   state j of the vacation process */
	Matrix F, F_new;
	
	create_matrix(&F, Vacation_ptr->rows, Vacation_ptr->cols);
	create_matrix(&F_new, Vacation_ptr->rows, Vacation_ptr->cols);
	
	//set F = Vacation
	copy_matrix(Vacation_ptr, &F);	
	//this is the probability of going from state 0,0 to 0,0 in one time quanta
	sum = get_matrix_element(&F, 0, 0);
	alpha_cdf = 0;
	phi = sum*curr_alpha;
	n = 2;
		
	while (sum < (1-PHI_ERROR_MEASURE))
	{		
		for (i = 0; i < full_vacation_order; i++)
		{	for (j = 0; j < full_vacation_order; j++)
			{
				total = 0;
				for (m = 0; m < full_vacation_order; m++)
				{
					total += get_matrix_element(Vacation_ptr, i, m)*get_matrix_element(&F, m, j);
				}
				get_matrix_element(&F_new, i, j) = total - get_matrix_element(&F, j, j)*get_matrix_element(Vacation_ptr, i, j);				
			}
		}
		//add the probability of finishing vacation in n steps to the sum
		sum += get_matrix_element(&F_new, 0, 0);
		//calculate the probability of the first arrival in step n
		alpha_cdf += pow((1-curr_alpha), (double)(n-1));
		//calculate the new value of phi
		phi += get_matrix_element(&F_new, 0, 0)*curr_alpha*alpha_cdf;
		//increment the number of steps for the next iteration		
		n++;
		// F = F_new
		copy_matrix(&F_new, &F);	
	}
	
	free_matrix(&F);
	free_matrix(&F_new);
	
	return phi;
}

/*******************
convergence_algorithm

Description: This function determines the vacation processes by an algorithm that is iterated until the values
			 of Phi converge.
			 
Parameters: Vacation_ptr - this is a pointer to the matrix containing the entire vacation processes
			S_ptr - this is a pointer to the service matrix
			S0_ptr - this is a pointer to the end of the process vector
			Beta_ptr - this is a pointer to the beginning of the process vector
			S0_mult_Beta_ptr - this is a pointer to a matrix which is S0*Beta
			Alpha_ptr - this is a pointer to a vector containing the probablity of memory arrivals for each processor
			Phi_ptr - this is a pointer to a vector containing the probability of a memory arrival within the vacation time for each processor
			num_procs - this is the number of total processors in the system
			service_length - this is the length of the service process
			vacation_order - this is the length of the vacation process
			
Returns: the number of iterations through the loop
*******************/
int convergence_algorithm(Matrix* Vacation_ptr, Matrix *S_ptr, Matrix *S0_ptr, 
							 Matrix *Beta_ptr, Matrix *S0_mult_Beta_ptr, Matrix *Alpha_ptr, Matrix *Phi_ptr, 
							 int num_procs, int service_length, int vacation_order)
{
	int num_loops = 0; //this is the number of iterations needed for convergence
	int i;
	double total_error;
	Matrix Old_phi;
		
	create_matrix(&Old_phi, 1, Phi_ptr->cols);
	
	total_error = 1000; //set the error to a large value to ensure that the while loop is entered the first time
	
	while (total_error > CONVERGE_ERROR)
	{
		//copy the phi vector to the old phi vector
		copy_matrix(Phi_ptr, &Old_phi);
		total_error = 0;
		
		for (i = 0; i < num_procs; i++)
		{
			create_vacation_process(Vacation_ptr, S_ptr, S0_ptr, Beta_ptr, S0_mult_Beta_ptr, Alpha_ptr, Phi_ptr, 
							 		num_procs, i, service_length, vacation_order, 0);
			//based on the vacation process, determine the new value of phi
			get_matrix_element(Phi_ptr, 0, i) = calc_phi(Vacation_ptr, get_matrix_element(Alpha_ptr, 0, i), vacation_order+1);		
			//now determine the difference between the new value of phi and the previous value of phi for the current processor
			if (get_matrix_element(Phi_ptr, 0, i) < get_matrix_element(&Old_phi, 0, i))
			{
				total_error += get_matrix_element(&Old_phi, 0, i) - get_matrix_element(Phi_ptr, 0, i);
			}
			else
			{
				total_error += get_matrix_element(Phi_ptr, 0, i) - get_matrix_element(&Old_phi, 0, i);
			}						
		}		
		num_loops++;
	}
	
	free_matrix(&Old_phi);
	
	return num_loops;
}
	
/*******************
build_probability_transition_matrix

Description: This function puts together the full probability transition matrix for a particular
			 processor.
			 
Parameters: P_ptr - this is a pointer to thge probability transition matrix that is to be built
			Vacation_ptr - this is a pointer to the matrix containing the entire vacation process for the current processor
			S_ptr - this is a pointer to the service matrix
			S0_ptr - this is a pointer to the end of the process vector
			Beta_ptr - this is a pointer to the beginning of the process vector
			S0_mult_Beta_ptr - this is a pointer to a matrix which is S0*Beta
			curr_alpha - this is the alpha value for the current processor
			Alpha_ptr - this is a pointer to a vector containing the probablity of memory arrivals for each processor
			Phi_ptr - this is a pointer to a vector containing the probability of a memory arrival within the vacation time for each processor
			service_length - this is the length of the service process
			vacation_order - this is the length of the vacation process
			
Returns: the order of the P matrix (which is a square matrix)
*******************/
int build_probability_transition_matrix(Matrix* P_ptr, Matrix* Vacation_ptr, Matrix *S_ptr, Matrix *S0_ptr, 
							 Matrix *Beta_ptr, Matrix *S0_mult_Beta_ptr, double curr_alpha, 
							 int service_length, int vacation_order)
{
	Matrix V0_mult_Beta;
	
	create_matrix(&V0_mult_Beta, vacation_order, service_length);
	
	//first set all elements to 0
	all_zeros(P_ptr);
	
	/*
	This is the Matlab version of the matrix that will be built
	P = [ 0               alpha*Beta          (1-alpha)*v_st          zeros(1,orderV);
          (1-alpha)*S0    alpha*(S0*Beta)+S   zeros(K,orderV)         zeros(K,orderV);
          (1-alpha)*V0    alpha*(V0*Beta)     (1-alpha)*V             alpha*V;
          zeros(orderV,1) V0*Beta            zeros(orderV,orderV)    V ];
    
	*/
	
	//this implements alpha*Beta 
	scalar_matrix_block_mult(curr_alpha, Beta_ptr, Beta_ptr->rows, Beta_ptr->cols, 0, 0, 
							P_ptr, 0, 1);
	
	//this implements (1-alpha)*v_st
	scalar_matrix_block_mult((1-curr_alpha), Vacation_ptr, 1, vacation_order, 0, 1, 
							P_ptr, 0, 1 + service_length);
	
	//this implements (1-alpha)*S0
	scalar_matrix_block_mult((1-curr_alpha), S0_ptr, S0_ptr->rows, S0_ptr->cols, 0, 0, 
							P_ptr, 1, 0);
							
	//this implements alpha*(S0*Beta)
	scalar_matrix_block_mult(curr_alpha, S0_mult_Beta_ptr, S0_mult_Beta_ptr->rows, S0_mult_Beta_ptr->cols, 0, 0, 
							P_ptr, 1, 1);
	
	//this adds S to alpha*(S0*Beta) to create alpha*(S0*Beta)+S
	matrix_block_addition(service_length, service_length, P_ptr, 1, 1, S_ptr, 0, 0, P_ptr, 1, 1);
	
	//this implements (1-alpha)*V0
	scalar_matrix_block_mult((1-curr_alpha), Vacation_ptr, vacation_order, 1, 1, 0, P_ptr, 1+service_length, 0);
	
	//this creates V0*Beta
	matrix_block_mult(Vacation_ptr, 1, 0, vacation_order, 1, Beta_ptr, 0, 0, Beta_ptr->rows, Beta_ptr->cols, 
					  &V0_mult_Beta, 0, 0);

	//this implements alpha*(V0*Beta)
	scalar_matrix_block_mult(curr_alpha, &V0_mult_Beta, vacation_order, service_length, 0, 0, P_ptr, 1+service_length, 1);
	
	//this implements (1-alpha)*V
	scalar_matrix_block_mult((1-curr_alpha), Vacation_ptr, vacation_order, vacation_order, 1, 1, P_ptr, 1+service_length, 1+service_length);
	
	//this implements alpha*V
	scalar_matrix_block_mult(curr_alpha, Vacation_ptr, vacation_order, vacation_order, 1, 1, P_ptr, 1+service_length, 1+service_length+vacation_order);
	
	//this places V0*Beta in the right spot, since P_ptr will have all zeros in this location, V0*Beta can be addded to all zeros to copy it
	matrix_block_addition(vacation_order, service_length, &V0_mult_Beta, 0, 0, P_ptr, 1+service_length+vacation_order, 1, P_ptr, 1+service_length+vacation_order, 1);
		
	//this places V in the correct spot
	matrix_block_addition(vacation_order, vacation_order, Vacation_ptr, 1, 1, P_ptr, 1+service_length+vacation_order, 1+service_length+vacation_order, 
						  P_ptr, 1+service_length+vacation_order, 1+service_length+vacation_order);
	
	free_matrix(&V0_mult_Beta);
	
	return (1 + service_length + 2*vacation_order);
}						 

/*******************
find_stationary_prob_vector

Description: This function determines the stationary probability vector of the matrix P, and stores 
			 it in the matrix Pi. Pi will then contain the probability of the being in each particular
			 state of P.
			 
Parameters: Pi_ptr - this is a vector containing the stationary probability
			P_ptr - this is a pointer to thge probability transition matrix that is to be built
			p_order - this is the ordewr of the square P matrix
			
Returns: nothing
*******************/
void find_stationary_prob_vector(Matrix* Pi_ptr, Matrix* P_ptr, int p_order)
{
	Matrix Pi_prime;
	Matrix* Pi_prime_ptr;
	
	double total_error;
	
	int i;
	
	create_matrix(&Pi_prime, 1, Pi_ptr->cols);
	Pi_prime_ptr = &Pi_prime;
	
	//first set all elements to 0
	all_zeros(Pi_ptr);
	//then set the first element to 1
	get_matrix_element(Pi_ptr, 0, 0) = 1;
	
	//set the error to a large value to ensure the loop enters the first time
	total_error = 1000;
	
	//loop until the stationary probability matrix has converged
	while (total_error > PI_ERROR)
	{
		//Pi_prime = Pi * P
		matrix_block_mult(Pi_ptr, 0, 0, 1, p_order, P_ptr, 0, 0, p_order, p_order,
					   	  Pi_prime_ptr, 0, 0);
					   	  
		total_error = 0;
					   	  
		for (i=0; i < p_order; i++)
		{
			//determine the difference between the Pi and Pi_prime to see if they have converged
			if (get_matrix_element(Pi_ptr, 0, i) < get_matrix_element(Pi_prime_ptr, 0, i))
			{
				total_error += get_matrix_element(Pi_prime_ptr, 0, i) - get_matrix_element(Pi_ptr, 0, i);
			}
			else
			{
				total_error += get_matrix_element(Pi_ptr, 0, i) - get_matrix_element(Pi_prime_ptr, 0, i);
			}
		}

		//now copy the Pi_prime matrix to the Pi matrix; Pi = Pi_prime
		copy_matrix(Pi_prime_ptr, Pi_ptr);
	}
	
	free_matrix(Pi_prime_ptr);
}

/*******************
determine_prob_of_waiting_for_others

Description: This function determines the probability of waiting for memory access because other
			 processors were currently accessing global memory from the stationary probability vector.
			 
Parameters: Pi_ptr - this is a vector containing the stationary probability
			service_length - this is the length of the service process
			vacation_order - this is the length of the vacation process
			prob_trans_matrix_order - this is the order of the probability transition matrix
			
Returns: the probability of this processor waiting for memory access due to memory contention
*******************/
double determine_prob_of_waiting_for_others(Matrix* Pi_ptr, int service_length, int vacation_order, int prob_trans_matrix_order)
{
	double prob_of_waiting_for_others = 0;
	int i;
	
	for (i = 1 + service_length + vacation_order; i < prob_trans_matrix_order; i++)
	{
		prob_of_waiting_for_others += get_matrix_element(Pi_ptr, 0, i);
	}
	
	return (prob_of_waiting_for_others);
}

/*******************
determine_start_times

Description: This function determines the start times for all of the tasks based on their
			 scheduling and task execution times.
			 
Parameters: Task_start_times_ptr - this is the matrix that contains the task start times
			Task_schedule_ptr - this is the matrix containing the task scheduling
			Task_durations_ptr - this is the matrix containing the task execution times 
			
Returns: nothing
*******************/
void determine_start_times(Matrix *Task_start_times_ptr, Matrix *Task_schedule_ptr, Matrix *Task_durations_ptr)
{
	int i, j;
	double curr_time;
	
	for (i = 0; i < Task_schedule_ptr->rows; i++)
	{
		j = 0;
		curr_time = 0;
		while ((j < Task_schedule_ptr->cols) && (get_matrix_element(Task_schedule_ptr, i, j) != EMPTY_SCHEDULE_SLOT))
		{
			get_matrix_element(Task_start_times_ptr, 0, (int)get_matrix_element(Task_schedule_ptr, i, j)) = curr_time;
			curr_time += get_matrix_element(Task_durations_ptr, 0, (int)get_matrix_element(Task_schedule_ptr, i, j));
			j++;
		}
	}
}

/*******************
assign_alphas

Description: This function determines which tasks are currently being analyzed based on the 
			 curr_time, and assigns the associated Alphas from the current tasks to the 
			 Alpha matrix.
			 
Parameters: Alpha_ptr - This is the Matrix that contains the arrival probabilities 
			Task_schedule_ptr - this is the matrix containing the task scheduling
			Task_alphas_ptr - this is the vector containing the alphas to be analyzed
			Curr_Tasks - this is a vector containing the task number of the tasks being analyzed
			Prev_analyzed_ratio_ptr - this is the vector containing the ratios that each of the
						tasks have been analyzed.
			
Returns: the number of tasks to be analyzed in the current partition.
*******************/
int assign_alphas(Matrix* Alpha_ptr, Matrix *Task_schedule_ptr, Matrix *Task_alphas_ptr, 
					Matrix *Curr_tasks, Matrix *Prev_analyzed_ratio_ptr)
					
{
	int i, j;
	int current_task;
	int num_curr_tasks;
	
	num_curr_tasks = 0;
		
	for (i = 0; i < Task_schedule_ptr->rows; i++)
	{
		j = 0;
		
		current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		
		//while the index is less than the size of the Task_schedule_ptr AND
		//the task is not the empty schedule slot AND
		//the current task has been completely analyzed
		while ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT) 
			   && (get_matrix_element(Prev_analyzed_ratio_ptr, 0, current_task) >= ANALYZED_RATIO_MAX))
		{			
			j++;
			current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		}
	
		//check to see that there are still tasks that are not finished executing 
		if ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT))
		{
			get_matrix_element(Alpha_ptr, 0, num_curr_tasks) = get_matrix_element(Task_alphas_ptr, 0, current_task);
			get_matrix_element(Curr_tasks, 0, num_curr_tasks) = current_task; 
			num_curr_tasks++;
		}	
	}
	
	return num_curr_tasks;
}

/*******************
update_task_parameters

Description: This function updates the parameters for the tasks that were just analyzed to account
			 for the analysis of the partition.
			 
Parameters: Task_start_times_ptr - this is the matrix that contains the task start times
			Task_schedule_ptr - this is the matrix containing the task scheduling
			Task_durations_ptr - this is the matrix containing the task execution times 
			curr_time - this is the current time of the anaylsis
			Curr_Tasks - this is a vector containing the task number of the tasks being analyzed
			num_curr_tasks - this is the number of current tasks that were analyzed 
							 in the current partition.
			Prob_of_wait_while_others_served_ptr - this is a vector holding the probability of each
							task waiting for memory access.
			Prev_analyzed_ratio_ptr - this is the vector contain the ratio of analysis for each task
			
Returns: the next start time of the partition to be analyzed
*******************/
double update_task_parameters(Matrix *Task_start_times_ptr, Matrix *Task_schedule_ptr, 
							  Matrix *Task_durations_ptr, double curr_time, Matrix *Curr_tasks,
							  int num_curr_tasks, Matrix *Prob_of_wait_while_others_served_ptr,
							  Matrix *Prev_analyzed_ratio_ptr, Matrix *Original_task_durations_ptr)
{
	double partition_time;
	double temp;
	double finish_time;
	double ratio_of_task_analyzed;
	double task_time_remaining;
	
	int i;
	int curr_task;
	
	//initialize the partition time to some extremely large number so that the first pass will be less
	partition_time = 1E12; 	
	
	for (i = 0; i < num_curr_tasks; i++)
	{
		//get the current task number
		curr_task = (int)get_matrix_element(Curr_tasks, 0, i); 
		temp = get_matrix_element(Original_task_durations_ptr, 0, curr_task)*(1 - get_matrix_element(Prev_analyzed_ratio_ptr, 0, curr_task)) 
				/ (1 - get_matrix_element(Prob_of_wait_while_others_served_ptr, 0, i));
		
		//now check if the current task is the first to be finished		
		if (temp < partition_time)
		{
			partition_time = temp;
		} 		
	}
	
	//determine the finish time for the current partition
	finish_time = curr_time + partition_time;
	
	//printf("partition time = %f, curr time = %f, finish_time = %f\n", partition_time, curr_time, finish_time);
	
	for (i = 0; i < num_curr_tasks; i++)
	{
		//get the current task number
		curr_task = (int)get_matrix_element(Curr_tasks, 0, i);
		//calculate the ratio of the task analyzed in the current partition
		ratio_of_task_analyzed = partition_time*(1 - get_matrix_element(Prob_of_wait_while_others_served_ptr, 0, i)) 
								 / get_matrix_element(Original_task_durations_ptr, 0, curr_task);
								 
		//update the Prev_analyzed_ratio_ptr vector to account for the current partition
		get_matrix_element(Prev_analyzed_ratio_ptr, 0, curr_task) = get_matrix_element(Prev_analyzed_ratio_ptr, 0, curr_task) + ratio_of_task_analyzed;
		
		//determine the time remaining that need to be analyzed	
		task_time_remaining = get_matrix_element(Original_task_durations_ptr, 0, curr_task)*(1 - get_matrix_element(Prev_analyzed_ratio_ptr, 0, curr_task));
		
		//update the task duration
		get_matrix_element(Task_durations_ptr, 0, curr_task) = finish_time - get_matrix_element(Task_start_times_ptr, 0, curr_task) + task_time_remaining;		
	}
		
	//now determine the new start times for each of the tasks
	determine_start_times(Task_start_times_ptr, Task_schedule_ptr, Task_durations_ptr);
	
	return (finish_time);
}

void print_task_times_and_analyzed_ratios(Matrix *Task_start_times_ptr, Matrix *Task_durations_ptr, 
										  Matrix *Prev_analyzed_ratio_ptr)
{
	int i = 0;
	
	printf("Task     Start      Duration      Analyzed Ratio\n"); 
	
	for (i = 0; i < Task_start_times_ptr->cols; i++)
	{
		printf("%d   %f    %f    %f\n", i, get_matrix_element(Task_start_times_ptr, 0, i), get_matrix_element(Task_durations_ptr, 0, i), get_matrix_element(Prev_analyzed_ratio_ptr, 0, i));
	}
}

/*******************
calc_total_processing_period

Description: This function calculates the total processing period length.
			 
Parameters: Task_start_times_ptr - this is the matrix that contains the task start times
			Task_schedule_ptr - this is the matrix containing the task scheduling
			Task_durations_ptr - this is the matrix containing the task execution times 
			Proc_total_task_times_ptr - this is a vector containing the total execution time
									of each processor in the processing period. The largest
									of these values is the total processing period.
			
Returns: the total processing period
*******************/
double calc_total_processing_period(Matrix *Task_schedule_ptr, Matrix *Task_start_times_ptr, Matrix *Task_durations_ptr, 
									Matrix *Proc_total_task_times_ptr)
{
	int i, j, current_task;
	double period;
	
	//first clear the vector
	all_zeros(Proc_total_task_times_ptr);
	
	period = 0;
		
	for (i = 0; i < Task_schedule_ptr->rows; i++)
	{
		j = 0;
		
		current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		
		//while the index is less than the size of the Task_schedule_ptr AND
		//the task is not the empty schedule slot
		while ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT))			   
		{				
			get_matrix_element(Proc_total_task_times_ptr, 0, i) = get_matrix_element(Task_start_times_ptr, 0, current_task) + get_matrix_element(Task_durations_ptr, 0, current_task);
			j++;
			current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		}
		
		if (get_matrix_element(Proc_total_task_times_ptr, 0, i) > period)
		{
			period = get_matrix_element(Proc_total_task_times_ptr, 0, i);
		}
	}	
	
	return period;
}

/*******************
perturb_task_schedule

Description: This function makes a small change to the task schedule. This works by choosing
			 one of the tasks at random, and then moving the task in one of 8 possible 
			 directions: up-left, up, up-right, left, right, down-left, down, or down-right.
			 In some cases all options may not be valid. Movements within the same row,
			 (i.e. same processor) require swapping with another task, where movements to
			 different rows do not require swapping. An example of a task movement is
			 shown.  Given the initial schedule of:
			 			 			
			 Proc0:  Task0 Task1 Task2
			 Proc1:  Task3 Task4 Task5
			 Proc2:  Task6 Task7 Task8
			 
			 If task 1 was the task chosen to be moved, then a left movement would look like:
			 
			 Proc0:  Task1 Task0 Task2
			 Proc1:  Task3 Task4 Task5
			 Proc2:  Task6 Task7 Task8
			 
			 A down-right movement would look like:
			 
			 Proc0:  Task0 Task2
			 Proc1:  Task3 Task4 Task1 Task5 
			 Proc2:  Task6 Task7 Task8
			 
			 An up-left movement would look like:
			 
			 Proc0:  Task0 Task2
			 Proc1:  Task3 Task4 Task5
			 Proc2:  Task1 Task6 Task7 Task8
			 
			 Wrap-around movements are possible for up and down, but not for left and right.
			 A task can be moved to a blank slot, if the slot before it is filled, but not
			 if the slot in front is empty.  For example given the original example schedule
			 if task2 is moved down-right, the result would be:
			 
			 Proc0:  Task0 Task1 
			 Proc1:  Task3 Task4 Task5 Task2
			 Proc2:  Task6 Task7 Task8
			 
			 Which is valid. But if task2 was chosen again at the next iteration, a down-right 
			 movement would not be allowed, because it would result in a gap between Task8 and
			 Task2, which is not allowed.
			 
Parameters: Task_schedule_ptr - this is the matrix containing the task scheduling
			num_tasks - this is the total number of tasks in the schedule
			
Returns: nothing
*******************/
void perturb_task_schedule(Matrix *Task_schedule_ptr, int num_tasks)
{
	int i, j;
	int curr_task_row, curr_task_col, new_task_row, new_task_col;
	int perturbed_task, current_task, swapped_task;
	int direction;
	double rand_num;
	int valid_move, task_found;
	double prob_sum;
	
	//generate a random number
	rand_num = ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) ) * num_tasks;
	perturbed_task = (int)rand_num;
	
	while (perturbed_task >= num_tasks)
	{
		perturbed_task--; 
	}
	
	//find the row and column of the chosen task
	task_found = 0;
	i = 0;
	while (task_found == 0) //this assumes that the task will definitely exist
	{
		j = 0;
		
		current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
				
		//while the index is less than the size of the Task_schedule_ptr AND
		//the task is not the empty schedule slot AND
		//the task has not been found
		while ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT) && (task_found == 0))			   
		{			
			if (current_task == perturbed_task)
			{	
				task_found = 1;
				curr_task_row = i;
				curr_task_col = j;
			}			
			j++;
			current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		}
		i++;
	}
	
	valid_move = 0; //this is a flag to specify if the move chosen is valid
	
	//repeat until a valid move is chosen
	while (valid_move == 0)
	{
		rand_num = ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) ) * 4;
		
		//determine the direction to be moved
		direction = (int)rand_num;
		if (direction >= 4)
		{
			direction = 3;
		}
		
		valid_move = 1;
		
		switch (direction)
		{			
			case 0: //up				
				//up is always valid, but need to find the valid coordinates of the new slot
			
				new_task_row = curr_task_row - 1;
				new_task_col = curr_task_col;
				
				//wrap around
				if (new_task_row < 0)
				{
					new_task_row = Task_schedule_ptr->rows-1;
				} 
				
				//move the column to the next available spot
				while (((new_task_col > 0) && ((int)get_matrix_element(Task_schedule_ptr, new_task_row, new_task_col-1) == EMPTY_SCHEDULE_SLOT)))
				{
					new_task_col--;
				}				
				
				printf("Made an up move with task %d\n", perturbed_task); 				
				break;			
			case 1: //left move
				new_task_row = curr_task_row;
				new_task_col = curr_task_col - 1;
							
				//check the invalid cases			
				if (new_task_col < 0)
				{
					valid_move = 0;
				}				
				
				if (valid_move == 1)
				{
					printf("Made a left move with task %d\n", perturbed_task); 
				}
				break;
			case 2: //right move
				new_task_row = curr_task_row;
				new_task_col = curr_task_col + 1;
							
				//check the invalid cases			
				if ((new_task_col >= Task_schedule_ptr->cols) ||
					((int)get_matrix_element(Task_schedule_ptr, new_task_row, new_task_col) == EMPTY_SCHEDULE_SLOT))
				{
					valid_move = 0;
				}
				
				if (valid_move == 1)
				{
					printf("Made a right move with task %d\n", perturbed_task); 
				}
				break;
			case 3: //down move
				//a down move is always valid, but need to find the valid coordinates of the new slot
			
				new_task_row = curr_task_row + 1;
				new_task_col = curr_task_col;
				
				//wrap around
				if (new_task_row >= Task_schedule_ptr->rows)
				{
					new_task_row = 0;
				} 
				
				//move the column to the next available spot
				while (((new_task_col > 0) && ((int)get_matrix_element(Task_schedule_ptr, new_task_row, new_task_col-1) == EMPTY_SCHEDULE_SLOT)))
				{
					new_task_col--;
				}
				printf("Made a down move with task %d\n", perturbed_task); 				
				break;			
			default:
				printf("This should never happen\n");						
		}
	}
	
	//in the case of a right or left movement, the task will be switched with another task
	if (new_task_row == curr_task_row)
	{
		swapped_task = (int)get_matrix_element(Task_schedule_ptr, new_task_row, new_task_col);
		get_matrix_element(Task_schedule_ptr, new_task_row, new_task_col) = (double)perturbed_task;
		get_matrix_element(Task_schedule_ptr, curr_task_row, curr_task_col) = (double)swapped_task;	
	}
	else //this means that a movement to another row (processor) will be done
	{
		//first shift all of the tasks to the right of the new task in the new row to the right to make room for the moved task
		for (i = Task_schedule_ptr->cols-1; i > new_task_col; i--)
		{
			//shift the task to the right
			get_matrix_element(Task_schedule_ptr, new_task_row, i) = get_matrix_element(Task_schedule_ptr, new_task_row, i-1);
		}
		
		//now assign the moved task to its new spot
		get_matrix_element(Task_schedule_ptr, new_task_row, new_task_col) = (double)perturbed_task;
		
		//now shift all the tasks in the old row to the left
		for (i = curr_task_col; i < (Task_schedule_ptr->cols - 1); i++)
		{
			get_matrix_element(Task_schedule_ptr, curr_task_row, i) = get_matrix_element(Task_schedule_ptr, curr_task_row, i+1);
		} 
				
		get_matrix_element(Task_schedule_ptr, curr_task_row, (Task_schedule_ptr->cols - 1)) = EMPTY_SCHEDULE_SLOT;		
	}	
}

/*******************
analyze_current_processing_period

Description: This function calculates the total processing period length.
			 
Parameters: num_procs - this is the number of processors
			num_tasks - this is the number of tasks
			service_length - this is the length of the service process
			Orig_taks_durations_ptr - this is a matrix containing the original task durations
			Task_schedule_ptr - this is the matrix containing the task scheduling
			Task_start_times_ptr - this is the matrix that contains the task start times
			Task_alphas_ptr - this is a matrix constining the alphas for each of the tasks
			S_ptr,
			S0_ptr,
			Beta_ptr,
			S0_mult_Beta_ptr,
			Proc_total_task_times_ptr - this is a vector containing the total execution time
									of each processor in the processing period. The largest
									of these values is the total processing period.
			
Returns: the total processing period
*******************/
double analyze_current_processing_period(int num_procs, int num_tasks, int service_length, Matrix* Orig_task_durations_ptr, Matrix* Task_schedule_ptr,
										 Matrix* Task_alphas_ptr, Matrix* S_ptr, Matrix* S0_ptr, Matrix* Beta_ptr, Matrix* S0_mult_Beta_ptr,
										 Matrix *Proc_total_task_times_ptr)
{
	Matrix Alpha, Prev_analyzed_ratio, Curr_tasks, Prob_of_wait_while_others_served,
		   Vacation, Phi, P, Pi, Task_durations, Task_start_times;
	int num_curr_tasks, vacation_order, prob_trans_matrix_order, i, j, current_task;
	double curr_time, processing_period;
	
	create_matrix(&Prob_of_wait_while_others_served, 1, num_procs);	
	create_matrix(&Curr_tasks, 1, num_procs);
	create_matrix(&Prev_analyzed_ratio, 1, num_tasks);	
	create_matrix(&Task_durations, 1, num_tasks);
	create_matrix(&Task_start_times, 1, num_tasks);
	
	//this matrix contains the ratio of how much of a task was analyzed for each task 
	all_zeros(&Prev_analyzed_ratio);
	
	vacation_order = (service_length + 1) * (num_procs - 1);
	prob_trans_matrix_order = 1 + service_length + 2*vacation_order;
		
	create_matrix(&Alpha, 1, num_procs);
	create_matrix(&Phi, 1, num_procs);
				
	//initialize the start time
	curr_time = 0.0;
		
	//make a copy of the original task durations into the task durations
	copy_matrix(Orig_task_durations_ptr, &Task_durations);	
	
	//determines the initial start times for all of the tasks
	determine_start_times(&Task_start_times, Task_schedule_ptr, &Task_durations);	
	
	//determine the number of processors that have tasks that need to be analyzed in the current partition
	num_curr_tasks = assign_alphas(&Alpha, Task_schedule_ptr, Task_alphas_ptr, &Curr_tasks, &Prev_analyzed_ratio);
	
	//while there is more than one task to analyze in the current partition
	// if there is only 1 (or less) this means that there will be no more memory contention, and so no more analysis is required
	while (num_curr_tasks > 1)
	{
		#ifdef PRINT_TASK_INFO					  			   
		printf("Num current tasks = %d\n", num_curr_tasks);
		printf("Current time = %f\n", curr_time);
		printf("Current Tasks = ");		
		print_matrix_block(&Curr_tasks, 0, 0, 1, num_curr_tasks);
		#endif
		
		vacation_order = (service_length + 1) * (num_curr_tasks - 1);
		prob_trans_matrix_order = 1 + service_length + 2*vacation_order;
		
		//initialize the phi vector to the alpha vector
		copy_matrix(&Alpha, &Phi);
			
		create_matrix(&Vacation, vacation_order+1, vacation_order+1);
		create_matrix(&P, prob_trans_matrix_order, prob_trans_matrix_order);
		create_matrix(&Pi, 1, prob_trans_matrix_order);
			
		//implement the convergence algorithm to determine the value of Phi
		convergence_algorithm(&Vacation, S_ptr, S0_ptr, Beta_ptr, S0_mult_Beta_ptr, &Alpha, &Phi, num_curr_tasks, service_length, vacation_order);
	
		//now that the Phi vector has been determined, the stationary probability vector for each processor can be calculated
		for (i = 0; i < num_curr_tasks; i++)
		{
			//create the vacation matrix for processor i
			create_vacation_process(&Vacation, S_ptr, S0_ptr, Beta_ptr, S0_mult_Beta_ptr, &Alpha, &Phi, num_curr_tasks, i, service_length, vacation_order, 0);
						
			//put together the probability transistion matrix
			build_probability_transition_matrix(&P, &Vacation, S_ptr, S0_ptr, Beta_ptr, S0_mult_Beta_ptr, get_matrix_element(&Alpha, 0, i), service_length, vacation_order);
		
			//find the stationary probability vector of P
			find_stationary_prob_vector(&Pi, &P, prob_trans_matrix_order);		
		
			get_matrix_element(&Prob_of_wait_while_others_served, 0, i) = determine_prob_of_waiting_for_others(&Pi, service_length, vacation_order, prob_trans_matrix_order);		
		}
		
		free_matrix(&Vacation);	
		free_matrix(&P);
		free_matrix(&Pi);		
		
		//update the task parameters and set the current time to the beginning of the next partition
		curr_time = update_task_parameters(&Task_start_times, Task_schedule_ptr, &Task_durations, curr_time, &Curr_tasks,
							  			   num_curr_tasks, &Prob_of_wait_while_others_served, &Prev_analyzed_ratio, Orig_task_durations_ptr);
		
		#ifdef PRINT_TASK_INFO					  			   
		print_task_times_and_analyzed_ratios(&Task_start_times, &Task_durations, &Prev_analyzed_ratio);
		#endif
		
		//determine the number of processors that have tasks that need to be analyzed in the current partition
		num_curr_tasks = assign_alphas(&Alpha, Task_schedule_ptr, Task_alphas_ptr, &Curr_tasks, &Prev_analyzed_ratio);
	}
		
	processing_period = calc_total_processing_period(Task_schedule_ptr, &Task_start_times, &Task_durations, Proc_total_task_times_ptr);
	
	printf("Task schedule\n\n");
	for (i = 0; i < num_procs; i++)
	{
		j = 0;
		for (i = 0; i < Task_schedule_ptr->rows; i++)
		{
			j = 0;
		
			current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		
			printf("Proc %d: ", i);
			//while the index is less than the size of the Task_schedule_ptr AND
			//the task is not the empty schedule slot AND
			//the start time of the current task plus the duration of the current task is less than the current time
			while ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT))				   
			{			
				printf("%d ", current_task);
				j++;
				current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);				
			}
			printf("   Proc time = %f\n", get_matrix_element(Proc_total_task_times_ptr, 0, i));
		}
	}	
#ifdef PRINT_TASK_INFO
	
	printf("\nTask durations\n");
	for (i = 0; i < num_procs; i++)
	{
		j = 0;
		for (i = 0; i < Task_schedule_ptr->rows; i++)
		{
			j = 0;
		
			current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		
			printf("Proc %d: ", i);
			//while the index is less than the size of the Task_schedule_ptr AND
			//the task is not the empty schedule slot AND
			//the start time of the current task plus the duration of the current task is less than the current time
			while ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT))				   
			{			
				printf("%f ",get_matrix_element(&Task_durations, 0, current_task));
				j++;
				current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);				
			}
			printf("\n");
		}
	}
	
	printf("\nTask start times\n");
	for (i = 0; i < num_procs; i++)
	{
		j = 0;
		for (i = 0; i < Task_schedule_ptr->rows; i++)
		{
			j = 0;
		
			current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);
		
			printf("Proc %d: ", i);
			//while the index is less than the size of the Task_schedule_ptr AND
			//the task is not the empty schedule slot AND
			//the start time of the current task plus the duration of the current task is less than the current time
			while ((j < Task_schedule_ptr->cols) && (current_task != EMPTY_SCHEDULE_SLOT))				   
			{			
				printf("%f ",get_matrix_element(&Task_start_times, 0, current_task));
				j++;
				current_task = (int)get_matrix_element(Task_schedule_ptr, i, j);				
			}
			printf("\n");
		}
	}	
#endif
	
	free_matrix(&Prev_analyzed_ratio);
	free_matrix(&Alpha);
	free_matrix(&Phi);	
	free_matrix(&Prob_of_wait_while_others_served);
	free_matrix(&Curr_tasks);	
	free_matrix(&Task_durations);
	free_matrix(&Task_start_times);
	
	return (processing_period);
}

/*******************
generate_chaos

Description: This function creates chaotic values over the whole bifurcation map range.
			 
Parameters: Chaos_ptr - this is a matrix containing the chaos values
			Avg_chaos_val_ptr - this is a vector containing the average values of chaos

Returns: nothing
*******************/
void generate_chaos(Matrix* Chaos_ptr, Matrix* Avg_chaos_val_ptr)
{
	int i, j;
	double mu;	
	double offset, gain;
		
	offset = INIT_VAL;
	gain = INIT_VAL;
	
	for (i = 0; i < MAX_ITERATION; i++)
	{		
		if (i == 0)
		{
			get_matrix_element(Chaos_ptr, i, 0) = INIT_VAL;			
		}
		else
		{
			get_matrix_element(Chaos_ptr, i, 0) = get_matrix_element(Chaos_ptr, i-1, MAX_ITERATION-1);
		}
		
		get_matrix_element(Avg_chaos_val_ptr, 0, i) = get_matrix_element(Chaos_ptr, i, 0);
		
		mu = (((double)MAX_ITERATION - (double)i - (double)1) / ((double)MAX_ITERATION-(double)1))*((double)MAX_MU - (double)MIN_MU) + (double)MIN_MU;
		
		for (j = 1; j < MAX_ITERATION; j++)
		{
			get_matrix_element(Chaos_ptr, i, j) = N_const * get_matrix_element(Chaos_ptr, i, j-1) - 2 * tanh(mu*get_matrix_element(Chaos_ptr, i, j-1))
											*exp(-3*get_matrix_element(Chaos_ptr, i, j-1)*get_matrix_element(Chaos_ptr, i, j-1)); 			
			get_matrix_element(Avg_chaos_val_ptr, 0, i) = get_matrix_element(Avg_chaos_val_ptr, 0, i) + get_matrix_element(Chaos_ptr, i, j);
			
			if (offset > get_matrix_element(Chaos_ptr, i, j))
			{	
				//the offset will be the minimum value found
				offset = get_matrix_element(Chaos_ptr, i, j);
			}
		
			if (gain < get_matrix_element(Chaos_ptr, i, j))
			{
				//the gain will be the maximum value found for now, and then calculated at the end
				gain = get_matrix_element(Chaos_ptr, i, j);
			}
		}
		get_matrix_element(Avg_chaos_val_ptr, 0, i) /= MAX_ITERATION; 
	}
	
	gain = gain - offset;
	
	for (i = 0; i < MAX_ITERATION; i++)
	{
		for (j = 1; j < MAX_ITERATION; j++)
		{
			get_matrix_element(Chaos_ptr, i, j) = (get_matrix_element(Chaos_ptr, i, j) - offset) / gain;
		}
		get_matrix_element(Avg_chaos_val_ptr, 0, i) = (get_matrix_element(Avg_chaos_val_ptr, 0, i) - offset) / gain;		
	}		
}
							 
int main(int argc, char **argv) {
	int num_procs = 3; //this is the number of processors in the system
	int num_tasks = 16; //this is the number of tasks to be scheduled
	Matrix single_proc_exec_time, S, S0, Beta, S0_mult_Beta, Orig_task_durations, 
	Task_alphas, Task_schedule, New_task_schedule, Optimal_task_schedule,
	Proc_total_task_times, Chaos, Avg_chaos_val;			
	double service_prob = 0.95;
	int service_length = 18;
	int i, j, current_task, num_perturbs;
	double last_processing_period, curr_processing_period, optimal_processing_period;
	double prob_threshold, rand_num;
	int step, new_path, reject_count;
	FILE *walk_file_ptr;
	FILE *result_file_ptr;
	
	walk_file_ptr = fopen(WALK_FILE_NAME,"w");
	if (walk_file_ptr == NULL)
	{
		printf("Could not open the walk file\n");
		exit(0);			
	}
	
	result_file_ptr = fopen(RESULT_FILE_NAME,"w");
	if (result_file_ptr == NULL)
	{
		printf("Could not open the result file\n");
		exit(0);			
	}
	
	create_matrix(&S, service_length, service_length);
	create_matrix(&S0, service_length, 1);
	create_matrix(&Beta, 1, service_length);
	create_matrix(&S0_mult_Beta, service_length, service_length);	
	create_service_matrix(&S, &S0, &Beta, service_prob);
		
	create_matrix(&Orig_task_durations, 1, num_tasks);
	create_matrix(&Task_alphas, 1, num_tasks);
	create_matrix(&Task_schedule, num_procs, num_tasks);
	create_matrix(&New_task_schedule, num_procs, num_tasks);
	create_matrix(&Optimal_task_schedule, num_procs, num_tasks);	
	create_matrix(&Proc_total_task_times, 1, num_procs);
	create_matrix(&Chaos, MAX_ITERATION, MAX_ITERATION);
	create_matrix(&Avg_chaos_val, 1, MAX_ITERATION);	
		
	//the following is a list of the task times when executed on a single processor
	get_matrix_element(&Orig_task_durations, 0, 0) = 500;
	get_matrix_element(&Orig_task_durations, 0, 1) = 800;
	get_matrix_element(&Orig_task_durations, 0, 2) = 1500;
	get_matrix_element(&Orig_task_durations, 0, 3) = 700;
	get_matrix_element(&Orig_task_durations, 0, 4) = 250;
	get_matrix_element(&Orig_task_durations, 0, 5) = 350;
	get_matrix_element(&Orig_task_durations, 0, 6) = 3000;
	get_matrix_element(&Orig_task_durations, 0, 7) = 450;
	get_matrix_element(&Orig_task_durations, 0, 8) = 600;
	get_matrix_element(&Orig_task_durations, 0, 9) = 1800;
	get_matrix_element(&Orig_task_durations, 0, 10) = 100;
	get_matrix_element(&Orig_task_durations, 0, 11) = 1100;
	get_matrix_element(&Orig_task_durations, 0, 12) = 900;
	get_matrix_element(&Orig_task_durations, 0, 13) = 950;
	get_matrix_element(&Orig_task_durations, 0, 14) = 850;
	get_matrix_element(&Orig_task_durations, 0, 15) = 2100;
		
	//the following is a list of the task memory request arrivals probability
	get_matrix_element(&Task_alphas, 0, 0) = 0.25;
	get_matrix_element(&Task_alphas, 0, 1) = 0.05;
	get_matrix_element(&Task_alphas, 0, 2) = 0.01;
	get_matrix_element(&Task_alphas, 0, 3) = 0.025;
	get_matrix_element(&Task_alphas, 0, 4) = 0.35;
	get_matrix_element(&Task_alphas, 0, 5) = 0.05;
	get_matrix_element(&Task_alphas, 0, 6) = 0.04;
	get_matrix_element(&Task_alphas, 0, 7) = 0.15;
	get_matrix_element(&Task_alphas, 0, 8) = 0.075;
	get_matrix_element(&Task_alphas, 0, 9) = 0.025;
	get_matrix_element(&Task_alphas, 0, 10) = 0.1;
	get_matrix_element(&Task_alphas, 0, 11) = 0.08;
	get_matrix_element(&Task_alphas, 0, 12) = 0.095;
	get_matrix_element(&Task_alphas, 0, 13) = 0.16;
	get_matrix_element(&Task_alphas, 0, 14) = 0.09;
	get_matrix_element(&Task_alphas, 0, 15) = 0.06;	
		
	//first set all schedule slots to the empty schedule slot flag
	set_all_elements(&Task_schedule, EMPTY_SCHEDULE_SLOT);
	//now assign the tasks to their schedule slots	
	 							/*	 Processor,  Task	*/
	get_matrix_element(&Task_schedule, 		0, 		0   ) = 0;
	get_matrix_element(&Task_schedule, 		0, 		1   ) = 1;
	get_matrix_element(&Task_schedule, 		0, 		2   ) = 2;
	get_matrix_element(&Task_schedule, 		0, 		3   ) = 3;
	get_matrix_element(&Task_schedule, 		0, 		4   ) = 4;
	get_matrix_element(&Task_schedule, 		1, 		0   ) = 5;
	get_matrix_element(&Task_schedule, 		1, 		1   ) = 6;
	get_matrix_element(&Task_schedule, 		1, 		2   ) = 7;
	get_matrix_element(&Task_schedule, 		1, 		3   ) = 8;
	get_matrix_element(&Task_schedule, 		1, 		4   ) = 9;
	get_matrix_element(&Task_schedule, 		2, 		0   ) = 10;
	get_matrix_element(&Task_schedule, 		2, 		1   ) = 11;
	get_matrix_element(&Task_schedule, 		2, 		2   ) = 12;
	get_matrix_element(&Task_schedule, 		2, 		3   ) = 13;
	get_matrix_element(&Task_schedule, 		2, 		4   ) = 14;
	get_matrix_element(&Task_schedule, 		2, 		5   ) = 15;
	
	
	//copy the starting schedule into the optimal to begin with
	copy_matrix(&Task_schedule, &Optimal_task_schedule);
				
	// the matrix S0_mult_Beta is the result of S0 multiplied by Beta, since this calculation is needed several times, it is done once and
	// the result is saved in this matrix to reduce the processing overhead
	matrix_mult(&S0, &Beta, &S0_mult_Beta);	
	
	//first calculate the processing period of the original schedule
	optimal_processing_period = analyze_current_processing_period(num_procs, num_tasks, service_length, &Orig_task_durations, &Task_schedule,
										 &Task_alphas, &S, &S0, &Beta, &S0_mult_Beta, &Proc_total_task_times);		
										 
	last_processing_period = optimal_processing_period;
	fprintf(result_file_ptr, "%f\n", last_processing_period);
	
	printf("Original processing period is: %f\n", optimal_processing_period);
	
	//generate the chaos values
	generate_chaos(&Chaos, &Avg_chaos_val);
	
	step = 0; //initialize the time
	new_path = 0; //this is a flag specifying if a new path will be chosen	
	reject_count = 0;
	
	while (step < MAX_ITERATION)
	{
		// a new path is taken if the path last time was rejected
		// in this case a number of steps according to the chaotic number will be chosen
		if (new_path == 1)
		{
			rand_num = ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) ) * MAX_ITERATION;
			if ((int)rand_num >= MAX_ITERATION)
			{
				rand_num -= 1;
			}
		
			printf("\n \n rand = %f, chaos = %f, avg = %f \n\n", rand_num, get_matrix_element(&Chaos, step, (int)rand_num), get_matrix_element(&Avg_chaos_val, 0, step));
			
			if (get_matrix_element(&Avg_chaos_val, 0, step) < 0.5)
			{			
				rand_num = ((fabs(get_matrix_element(&Avg_chaos_val, 0, step) - get_matrix_element(&Chaos, step, (int)rand_num))) / (1 - get_matrix_element(&Avg_chaos_val, 0, step))) * (MAX_NUM_MOVES - MIN_NUM_MOVES) + MIN_NUM_MOVES;
				num_perturbs = (int)rand_num;		
			}
			else
			{
				rand_num = ((fabs(get_matrix_element(&Avg_chaos_val, 0, step) - get_matrix_element(&Chaos, step, (int)rand_num))) / get_matrix_element(&Avg_chaos_val, 0, step)) * (MAX_NUM_MOVES - MIN_NUM_MOVES) + MIN_NUM_MOVES;
				num_perturbs = (int)rand_num;		
			}
		}
		else //only move one step because the last step was accepted (i.e. try to fins the local minimum)
		{			
			num_perturbs = 1;
		}	
		
		printf("\n num_perturbs = %d   step = %d \n", num_perturbs, step);  
		printf("\n The current processing period is: %f\n", curr_processing_period);	
		
		new_path = 0;
		
		copy_matrix(&Task_schedule, &New_task_schedule);
		
		fprintf(walk_file_ptr, "%d\n", num_perturbs);
		
		for (i = 0; i < num_perturbs; i++)
		{
			perturb_task_schedule(&New_task_schedule, num_tasks);
		}
		
		curr_processing_period = analyze_current_processing_period(num_procs, num_tasks, service_length, &Orig_task_durations, &New_task_schedule,
										 &Task_alphas, &S, &S0, &Beta, &S0_mult_Beta, &Proc_total_task_times);
			
		fprintf(result_file_ptr, "%f\n", curr_processing_period);
		if (curr_processing_period <= last_processing_period)
		{
			//if the new solution is the same or better than the last solution, then keep it
			copy_matrix(&New_task_schedule, &Task_schedule);
			last_processing_period = curr_processing_period;
			printf("Accepted the new solution\n");
			reject_count = 0;
			
			if (curr_processing_period <= optimal_processing_period)
			{
				//if the new solution is the best found so far, then copy it to the optimal solution
				optimal_processing_period = curr_processing_period;
				copy_matrix(&New_task_schedule, &Optimal_task_schedule);
			}
		}
		else //if the new solution is worse than the previous, than keep it with some probability related to the temperature
		{						
			prob_threshold = exp(-(curr_processing_period-last_processing_period)/((MAX_ITERATION-step+1)*K_const));
			rand_num = ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
		
			if (rand_num <= prob_threshold)
			{
				printf("Accepted the new solution with probability of acceptance of %f \n", prob_threshold);
				//keep the new solution (even though it is worse)
				copy_matrix(&New_task_schedule, &Task_schedule);
				last_processing_period = curr_processing_period;
				reject_count = 0;
			}	
			else
			{
				printf("Rejected the new solution with probability of acceptance of %f \n", prob_threshold);
				new_path = 1;//set the flag for a new path
				reject_count++;				
				
				if (reject_count >= RESTART_THRESHOLD)
				{
					//set the task schedule to the optimal task to restart searching
					copy_matrix(&Optimal_task_schedule, &Task_schedule);
					last_processing_period = optimal_processing_period;
					new_path = 0;
					reject_count = 0;
					printf("RESTARTED SEARCH\n");
				}
			}
		}
		step++; 
	}
	
	printf("\nThe very best processing period is: %f\n", optimal_processing_period);
	
	printf("Optimal Task schedule\n");
	for (i = 0; i < Optimal_task_schedule.rows; i++)
	{
		j = 0;
		
		current_task = (int)get_matrix_element(&Optimal_task_schedule, i, j);
		
		printf("Proc %d: ", i);
		//while the index is less than the size of the Task_schedule_ptr AND
		//the task is not the empty schedule slot 
		while ((j < Optimal_task_schedule.cols) && (current_task != EMPTY_SCHEDULE_SLOT))				   
		{			
			printf("%d ", current_task);
			j++;
			current_task = (int)get_matrix_element(&Optimal_task_schedule, i, j);				
		}
		printf("\n");	
	}
			
	free_matrix(&S);	
	free_matrix(&S0);
	free_matrix(&Beta);	
	free_matrix(&S0_mult_Beta);	
	free_matrix(&Task_alphas);
	free_matrix(&Task_schedule);
	free_matrix(&New_task_schedule);
	free_matrix(&Optimal_task_schedule);
	free_matrix(&Orig_task_durations);
	free_matrix(&Proc_total_task_times);
	free_matrix(&Chaos);
	free_matrix(&Avg_chaos_val);	
	
	fclose(walk_file_ptr);
	fclose(result_file_ptr);
			 		
	return 0;
}


