mpi_ga.c

MPI实现遗传算法

#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define POPSIZE		500
#define NVARS		6
#define MAXGENS		2000 
#define PXOVER		0.8             /* 杂交概率 */
#define PMUTATION	0.15            /* 变异概率 */

#define TASK_NUM(i)	((POPSIZE+pnum-1)/pnum)
//#define TASK_NUM(i)	((i)==pnum-1 && POPSIZE%((POPSIZE+pnum-1)/pnum)? POPSIZE%((POPSIZE+pnum-1)/pnum):(POPSIZE+pnum-1)/pnum)
#define EX_NUM		3

struct genotype
{
	double gene[NVARS]; 
	double fitness;            /* 适应度 */
	double rfitness;           /* 相对适应度 */
	double cfitness;           /* 累积适应度 */
} *population, *newpopulation;

double lowerbound[NVARS]; 
double upperbound[NVARS]; 
double local_best_individual[NVARS+1];	/* 局部最优解 */
double best_individual[NVARS+1];		/* 全局最优解 */

int generation;
double lower_fitness;
FILE *galog;               

int pid, pnum;

double randval(double low, double high)
{
	double val;
	val = ((double)(rand()%1000)/1000.0)*(high - low) + low;
	return val;
}

void initialize(void)
{
	FILE *infile;
	int i, j, r;
	double lbound, ubound;
	MPI_Status status;

	population = malloc(sizeof(struct genotype)*(TASK_NUM(pid)+1));
	newpopulation = malloc(sizeof(struct genotype)*(TASK_NUM(pid)+1));

	if (pid == pnum - 1) {
		if ((infile = fopen("ga_data.txt","r"))==NULL)
		{
			fprintf(galog,"\nCannot open input file!\n");
			exit(1);
		}

		srand(time(0));
		for (i=0; i<pnum; i++) {
			r = rand();
			MPI_Send(&r, 1, MPI_INT, i, 1, MPI_COMM_WORLD);
		}

		for (i = 0; i < NVARS; i++) {
			fscanf(infile, "%lf",&lowerbound[i]);
			fscanf(infile, "%lf",&upperbound[i]);
		}
	}
	else {
		MPI_Recv(&r, 1, MPI_INT, pnum-1, 1, MPI_COMM_WORLD, &status);
		srand(r);
	}

	MPI_Bcast(lowerbound, NVARS, MPI_DOUBLE, pnum-1, MPI_COMM_WORLD);
	MPI_Bcast(upperbound, NVARS, MPI_DOUBLE, pnum-1, MPI_COMM_WORLD);

	for (j = 0; j < TASK_NUM(pid); j++) {
		population[j].fitness  = 0;
		population[j].rfitness = 0;
		population[j].cfitness = 0;
		
		//printf("\nGene of member %d in %d is:\n", j, pid);
		for (i = 0; i < NVARS; i++) {
			population[j].gene[i] = randval(lowerbound[i], upperbound[i]);
			//printf("  %lf   ", population[j].gene[i]);
		}
	}
	if (pid == pnum - 1) fclose(infile);
}

void evaluate(void)
{
	int mem;
	int i;
	double x[NVARS+1];

	lower_fitness = 0.0;

	for (mem = 0; mem < TASK_NUM(pid); mem++) {
		for (i = 0; i < NVARS; i++)
			x[i+1] = population[mem].gene[i];

		population[mem].fitness = (x[1]*x[1]*x[1]*x[1]*x[1]*x[1])
			- (x[1]*x[2]*x[3]*x[4]*x[4]*x[6]) + (x[3]*x[3]*x[3]*x[3]*x[3]*x[5]*x[6])
			- (x[2]*x[4]*x[5]*x[5]*x[6]) - (x[1]*x[1]*x[3]*x[3]*x[5]*x[5])
			+ (x[2]*x[2]*x[3]*x[5]*x[6]) + (x[1]*x[2]*x[4]*x[4]*x[4]*x[5]);
		//printf(" The fitness of %d at %lf, %lf, %lf, %lf, %lf, %lf is %lf\n", mem, x[0], x[1], x[2], x[3], x[4], x[5], population[mem].fitness);
		if (population[mem].fitness < lower_fitness) lower_fitness = population[mem].fitness;
	}
}

void keep_the_best()
{
	int mem, i, cur_best = 0;

	for (mem = 0; mem < TASK_NUM(pid); mem++) {
		if (population[mem].fitness > population[TASK_NUM(pid)].fitness) {
			cur_best = mem;
			population[TASK_NUM(pid)].fitness = population[mem].fitness;
		}
	}
	//printf("\nG%4d: the best individual in Process %d is %d, the fitness is %lf.\n", generation, pid, cur_best, population[cur_best].fitness);
	for (i = 0; i < NVARS; i++) {
		population[TASK_NUM(pid)].gene[i] = population[cur_best].gene[i];
		local_best_individual[i+1] = population[cur_best].gene[i];
	}
	local_best_individual[0] = population[cur_best].fitness;
}

void pop_select(void)
{
	MPI_Status status;
	MPI_Request handle;
	int mem, i, j, k;
	double sum = 0;
	double p;
	static struct genotype ex_member[EX_NUM];

	for (mem = 0; mem < TASK_NUM(pid); mem++) {
		sum += (population[mem].fitness - lower_fitness);
		//printf(" Sum %d is %lf\n", mem, sum);
	}

	for (mem = 0; mem < TASK_NUM(pid); mem++) {
		population[mem].rfitness =  (population[mem].fitness - lower_fitness)/sum;
		//printf("The fitness of %d is %lf, the rfitness is %lf\n", mem, population[mem].fitness, population[mem].rfitness);
	}
	population[0].cfitness = population[0].rfitness;

	for (mem = 1; mem < TASK_NUM(pid); mem++) {
		population[mem].cfitness =  population[mem-1].cfitness + population[mem].rfitness;
	}

	for (i = 1; i <= EX_NUM; i++) {
		p = rand()%1000/1000.0;
		//printf("  p = %lf", p);
		//for (j=0; j<=TASK_NUM(pid); j++) printf(" %lf ", population[j].cfitness);
		//printf("\n");
		if (p < population[0].cfitness) {
			//printf("  %d send to %d: %lf, %lf, %lf, %lf, %lf, %lf: %lf\n", pid, (pid+i*generation)%pnum, 
			//	population[0].gene[0], population[0].gene[1], population[0].gene[2], population[0].gene[3],
			//	population[0].gene[4], population[0].gene[5], population[0].fitness);
			MPI_Isend(&population[0], sizeof(struct genotype)/sizeof(char), 
				MPI_CHAR, (pid+i*generation)%pnum, 0, MPI_COMM_WORLD, &handle);
		}
		else {
			for (j = 0; j < TASK_NUM(pid);j++)     
				if (p >= population[j].cfitness && p <  population[j+1].cfitness) {
					//printf("  %d send to %d: %lf, %lf, %lf, %lf, %lf, %lf: %lf\n", pid, (pid+i*generation)%pnum, 
					//	population[j+1].gene[0], population[j+1].gene[1], population[j+1].gene[2], population[j+1].gene[3],
					//	population[j+1].gene[4], population[j+1].gene[5], population[j+1].fitness);
					MPI_Isend(&population[j+1], sizeof(struct genotype)/sizeof(char), 
						MPI_CHAR, (pid+i*generation)%pnum, 0, MPI_COMM_WORLD, &handle);
					break;
				}
		}

		//printf("  %d receive from %d\n", pid, (pid+(pnum-i)*generation)%pnum);
		MPI_Recv (&ex_member[i-1], sizeof(struct genotype)/sizeof(char), MPI_CHAR, 
			(pid+(pnum-i)*generation)%pnum, 0, MPI_COMM_WORLD, &status);
		//printf("  %d receive from %d: %lf, %lf, %lf, %lf, %lf, %lf: %lf\n", pid, (pid+(pnum-i)*generation)%pnum, 
		//	ex_member[i-1].gene[0], ex_member[i-1].gene[1], ex_member[i-1].gene[2], ex_member[i-1].gene[3],
		//	ex_member[i-1].gene[4], ex_member[i-1].gene[5], ex_member[i-1].fitness);
	}

	for (i = 0; i < TASK_NUM(pid); i++)
	{
		p = rand()%1000/1000.0;
		if (p < population[0].cfitness)
			newpopulation[i] = population[0];     
		else {
			for (j = 0; j < TASK_NUM(pid); j++)     
				if (p >= population[j].cfitness &&
					p <  population[j+1].cfitness)
					newpopulation[i] = population[j+1];
		}
	}

	/* once a new population is created, copy it back */
	for (i = 0; i < TASK_NUM(pid); i++)
		population[i] = newpopulation[i];

	for (i = 0; i < EX_NUM; i++) {
		j = rand()%TASK_NUM(pid);
		if (population[j].fitness < ex_member[i].fitness) {
			for (k = 0; k < NVARS; k++) {
				population[j].gene[k] = ex_member[i].gene[k];
			}
			population[j].rfitness = 0;
			population[j].cfitness = 0;
			population[j].fitness = ex_member[i].fitness;
			//printf("  Update member %d: %lf, %lf, %lf, %lf, %lf, %lf: %lf\n", j, 
			//	population[j].gene[0], population[j].gene[1], population[j].gene[2], population[j].gene[3],
			//	population[j].gene[4], population[j].gene[5], population[j].fitness);
		}
	}
}

void swap(double *x, double *y)
{
	double temp;

	temp = *x;
	*x = *y;
	*y = temp;

}

void Xover(int one, int two)
{
	int i;
	int point; /* crossover point */

	if(NVARS > 1) {
		if(NVARS == 2)
			point = 1;
		else
			point = (rand() % (NVARS - 1)) + 1;

		for (i = 0; i < point; i++)
			swap(&population[one].gene[i], &population[two].gene[i]);
	}
}

void crossover(void)
{
	int i, mem, one;
	int first  =  0; /* count of the number of members chosen */
	double x;

	for (mem = 0; mem < TASK_NUM(pid); mem++) {
		x = rand()%1000/1000.0;
		if (x < PXOVER) {
			first++;
			if (first % 2 == 0)
				Xover(one, mem);
			else
				one = mem;
		}
	}
}

void mutate(void)
{
	int i, j;

	for (i = 0; i < TASK_NUM(pid); i++)
		for (j = 0; j < NVARS; j++)	{
			if (rand()%1000/1000.0 < PMUTATION) {
				/* find the bounds on the variable to be mutated */
				population[i].gene[j] = randval(lowerbound[j], upperbound[j]);
			}
		}
}

void elitist(void)
{
	int i;
	double best, worst; 
	int best_mem, worst_mem; 

	best = population[0].fitness;
	worst = population[0].fitness;
	for (i = 0; i < TASK_NUM(pid) - 1; ++i)
	{
		if(population[i].fitness > population[i+1].fitness)
		{     
			if (population[i].fitness >= best)
			{
				best = population[i].fitness;
				best_mem = i;
			}
			if (population[i+1].fitness <= worst)
			{
				worst = population[i+1].fitness;
				worst_mem = i + 1;
			}
		}
		else
		{
			if (population[i].fitness <= worst)
			{
				worst = population[i].fitness;
				worst_mem = i;
			}
			if (population[i+1].fitness >= best)
			{
				best = population[i+1].fitness;
				best_mem = i + 1;
			}
		}
	}

	if (best >= population[TASK_NUM(pid)].fitness)
	{
		for (i = 0; i < NVARS; i++) {
			population[TASK_NUM(pid)].gene[i] = population[best_mem].gene[i];
			local_best_individual[i+1] = population[best_mem].gene[i];
		}
		population[TASK_NUM(pid)].fitness = population[best_mem].fitness;
		local_best_individual[0] = population[best_mem].fitness;
	}
	else
	{
		for (i = 0; i < NVARS; i++)
			population[worst_mem].gene[i] = population[TASK_NUM(pid)].gene[i];
		population[worst_mem].fitness = population[TASK_NUM(pid)].fitness;
	}
	//printf(" Generation %d in Process %d: best = %lf, worst =%lf\n", generation, pid, local_best_individual[0], worst);
	//printf("\nG%4d: the best individual in Process %d is %d, the fitness is %lf.\n", generation, pid, best_mem, population[best_mem].fitness);
}

void report(void)
{
	int i, j;

	fprintf(galog, "\nGene data in %d for generation %d is:\n", pid, generation);
	for (j = 0; j < TASK_NUM(pid); j++) {
		fprintf(galog, "%4d  :", j);

		for (i = 0; i < NVARS; i++) {
			fprintf(galog, "\t%10.6lf  ", population[j].gene[i]);
		}
		fprintf(galog, ":   %20.6lf\n", population[j].fitness);
	}
}

void gene_max(double *in, double *inout, int *len, MPI_Datatype *dptr)
{
	int i;

	if (inout[0] < in[0]) {
		for (i=0; i < *len; ++i) {
			inout[i] = in[i];
		}
	}
}

int main(int argc, char **argv)
{
	double startwtime = 0.0, endwtime;
	int i, j, flag;
	int n = 0;
    MPI_Op my_op;

	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &pnum);
	MPI_Comm_rank(MPI_COMM_WORLD, &pid);

	fprintf(stdout, "Process %d of %d start ....\n",pid, pnum);

	//if ((galog = fopen("ga_log.txt","w"))==NULL)
	//{
	//	exit(1);
	//}
	galog = stdout;

	generation = 0;

	if (pid == pnum - 1) {
		srand(time(0));
	    startwtime = MPI_Wtime();
	}

	initialize();
	evaluate();
	keep_the_best();
	//report();
	while(generation<MAXGENS)
	{
		generation++;
		pop_select();
		crossover();
		mutate();
		//report();
		evaluate();
		elitist();
	}

    MPI_Op_create((MPI_User_function *)gene_max, 1, &my_op);
    MPI_Reduce(local_best_individual, best_individual, NVARS+1, MPI_DOUBLE, my_op, pnum-1, MPI_COMM_WORLD);

/*
	for (i = 0; i < NVARS; i++)
	{
		fprintf (galog,"\n var(%d) = %3.3f",i,population[POPSIZE].gene[i]);
	}
	fprintf(galog,"\n\n Best fitness = %3.3f",population[POPSIZE].fitness);
*/
	if (pid == pnum - 1) {
		fprintf(galog, "\n The solution is (population size = %d, # of generations = %d):\n", TASK_NUM(0)*pnum, generation);
		for (i=0; i<NVARS; i++) fprintf(galog, "    x[%d]   ", i);
		fprintf(galog, "        fitness   \n");
		for (i=0; i<NVARS; i++) fprintf(galog, "%10.6lf ", best_individual[i+1]);
		fprintf(galog, "%20.6lf\n", best_individual[0]);

		endwtime = MPI_Wtime();
		fprintf(galog, "\n Wall clock time =  %lf\n ", endwtime-startwtime);

		//fclose(galog);
	}

	//printf("\nProcess %d of %d Done ....\n", pid, pnum);
	fflush(stdout);

	MPI_Finalize();

	return 0;
}