arpso.c

一个遗传粒子群算法原代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <time.h>

#include <functions.c>
#define pi 3.14159265358979
#define e  2.71828182845904
#define TICKSIZE 5
#define bool int
#define true 1
#define false 0
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))

#define DBL_MAX pow(2,1023)
int outputDeltaEvals = 20;

int t;                              // timestep
int tMax;                           // max timestep
int swarmsize;                      // number of particles
int nParticles;						// actual number of particles (this number can vary because of fusion and fission)
int dim;                            // dimensionality of the function
int topo;                           // neighbourhood topology
int runs;                           // Runs
int runno;							// Run number
int procent;						//......
long evaluations=0;
int notImproved=0;
int particleLowerBound = 4;
int modeFusion = 1;
int tLastFission = -1;

int evalImproved;
int remGbest = -1;
double dir;
double dLow, dHigh;
double diversity;
double diagonal;

double w,wstart,wend;               // inertia weight
double theta, theta_start, theta_end, ignorance;// theta & ignorance
double phi1,phi2;                   // ...
double vmax;                        // ...
double minlimit,maxlimit;           // limits of the problem
double mininitlimit,maxinitlimit;   // limits of initialization for particles
int probID;                         //
double bfactor;                     // bounce factor
double radius;                      // particle radius
int colllistlength = 25;
int colllistPtr;
int *partvelchanged;
int nAllowed = 0;
int maxAllowed;
int windowSize=10;
int windowPos;
double fitnessThreshold = 0.1;
double *fusionPoint;

typedef struct statistics {
	int nEvals;
	double nParticles, avgbestfitness, avgavgfitness, stddev, posStddev;
} statistics;

statistics *logdata, *logdataRuns;

FILE *positionfile;
FILE *numberparticlesfile,*avgbestfitnessfile,*avgavgfitnessfile,*stddevfile,*posstddevfile,*besteverfile;

typedef struct par {
	double *pos;
	double *vel;
	double fitness;
	
	// -- bruges til m錶 for improvement over fitwindowsize timestep
	double *fitnesslist;
	double fitnessImprove;
	// ---
	double *best;
	double *bestvel;
	double bestfit;
	// ---
	double *gbest;
	double gbestfit;
	int *nlist;					    // indexlist of neighbours
	int changedvel;
	// --- collision
	int *colllist;					// list of number of collisions over the last colllistlength timesteps
	int colllistsum;					// sum of the colllist
	
	int oldVal;
	// ---
	int behavior;						// 0 is global; 1 is local ; ...
	int tImproved;					// timesteps since last fitness improvement
	int tChangedBehavior;				// timesteps since last behavior change
	// ---
	bool allowed;   // .. to be put to the gbest
	
	//-- Nuklear Stuff
	int parSize;
	int fusedWith;
} par;

par *swarm;
par *bestever;

double randval(double low, double high);
void allocspace();
void initvalues();
void updvel();
//void collcontrol();
void newpos();
void assignfitn();
void updneigh();
void report();
void updparams();
void logging();
void writestats();
void logstats();
void writepositions();
int overlap(int i, int j);
double dist(int i, int j);
double dist2(double *i, double *j);
double decreaseTheta(double a, double b, int time);
void randomizevectordir(int i);
void collectbestever();
void updColllist();
void fitImprovement();
double distanceToMeanPoint(int index);
int findFissionNumber();
void meanPoint(double *avgpoint);
void printpositions();
void printbestfitness();
int findFissionParticle(int ModelNr);

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

int mainParticle(int i) {
     return (swarm[i].fusedWith == i);
}

void allocspace() {
	int i;
	
	if (runs == 1 && dim == 2) {
		if ((positionfile = fopen("positions.dat","w")) == NULL) {
			printf("The file positions.dat couldn't be opened\n");
			exit(-1);
		}
	}
	
	if ((numberparticlesfile = fopen("numberofparticles.dat","w")) == NULL) {
		printf("The file numberofparticles.dat couldn't be opened\n");
		exit(-1);
	}
	
	if ((avgbestfitnessfile = fopen("avgbestfitness.dat","w")) == NULL) {
		printf("The file avgbestfitness.dat couldn't be opened\n");
		exit(-1);
	}
	
	if ((avgavgfitnessfile = fopen("avgavgfitness.dat","w")) == NULL) {
		printf("The file avgavgfitness.dat couldn't be opened\n");
		exit(-1);
	}
	
	if ((stddevfile = fopen("stddev.dat","w")) == NULL) {
		printf("The file stddev.dat couldn't be opened\n");
		exit(-1);
	}
	
	if ((posstddevfile = fopen("posstddev.dat","w")) == NULL) {
		printf("The file posstddev.dat couldn't be opened\n");
		exit(-1);
	}
	
	if ((besteverfile = fopen("bestever.dat","w")) == NULL) {
		printf("The file bestever.dat couldn't be opened\n");
		exit(-1);
	}
	
	logdata = malloc((tMax+1) * sizeof(statistics));
	// max number of evaluations is at most tmax+swarmsize
	// (if nEvals is tmax-1 before a "timestep" starts, then a whole new swarm can be evaluated)
	// we will only output each outputDeltaEvals evaluation-value, hence the formula for the size of logdataRuns
	logdataRuns = malloc(((tMax+swarmsize)/outputDeltaEvals+1) * sizeof(statistics));
	partvelchanged = malloc(swarmsize * sizeof(int));
	
	bestever = malloc((runs+1) * sizeof(par));
	
	for(i=0;i<runs+1;i++)
		bestever[i].pos = malloc(dim * sizeof(double));
	
	swarm = malloc((swarmsize+1) * sizeof(par));
	fusionPoint = malloc(dim * sizeof(double));
	
	for(i=0;i<swarmsize+1;i++) {
		swarm[i].pos = malloc(dim * sizeof(double));
		swarm[i].vel = malloc(dim * sizeof(double));
		// ---
		swarm[i].best = malloc(dim * sizeof(double));
		swarm[i].bestvel = malloc(dim * sizeof(double));
		// ---
		swarm[i].gbest = malloc(dim * sizeof(double));
		// ---
		swarm[i].colllist = malloc(colllistlength*sizeof(int));
		swarm[i].fitnesslist = malloc(windowSize * sizeof(double));
		
		switch(topo){
			
		case 0:
			// fully connected
			swarm[i].nlist = malloc(swarmsize * sizeof(int));
			break;
			/*
		case 1:
			// circle k=2
			swarm[i].nlist = malloc(3 * sizeof(int));
			break;
		case 2:
			if (i==0) swarm[i].nlist = malloc(swarmsize * sizeof(int));
			else swarm[i].nlist = malloc(2 * sizeof(int));
			break;
			*/
		}
	}
}

void initvalues() {
	int i,j;

	dir = 1;
	//dLow = 0.01;
	//dHigh = 0.1;
	diversity = DBL_MAX;
  diagonal = sqrt(dim) * (maxlimit - minlimit); // = sqrt(dim * pow(maxlimit - minlimit,2));

	evalImproved = 0;
	modeFusion = 1;
	tLastFission = -1;
	nParticles = swarmsize/1;// nParticles = swarmsize/initparticlessize
	w= wstart;
	if (runno==0) {
		for(i=0;i<swarmsize;i++)
			partvelchanged[i] = 0;
	}

	for(i=0;i<dim;i++)
	  fusionPoint[i] = 0;
	
	// stats
	for(i=0;i<tMax;i++) {
		logdata[i].nEvals = -1;
	}
	for(i=0;i<((tMax+swarmsize)/outputDeltaEvals+1);i++) {
		logdataRuns[i].nEvals = -1;
	}
	
	// position and velocity
	for(i=0;i<swarmsize+1;i++) {
		for(j=0;j<dim;j++) {
			swarm[i].pos[j] = randval(mininitlimit,maxinitlimit);
			swarm[i].vel[j] = randval(-vmax,vmax); 
		}
		swarm[i].fitness = DBL_MAX;
		swarm[i].bestfit = DBL_MAX;
		swarm[i].gbestfit = DBL_MAX;
		swarm[i].changedvel = 0;
		swarm[i].tImproved = 0;
		swarm[i].parSize = 1;    
		swarm[i].fitnessImprove = 0;
		// Should only be done in fusion model 1!!!
		if (i!=(swarmsize))
		  swarm[i].fusedWith = i ;//% (swarmsize/2);
		else
		  swarm[i].fusedWith = i; 
	}
	
	
	
	// initializing collisionlist
	colllistPtr = 0;
	for(i=0;i<swarmsize;i++) {
		swarm[i].colllistsum = 0;
		for(j=0;j<colllistlength;j++) {
			swarm[i].colllist[j] = 0;
		}
	}
	
	switch(topo){
	case 0:
		// fully connected
		for(i=0;i<swarmsize;i++){
			for(j=0;j<swarmsize;j++)
				swarm[i].nlist[j]=j;
		}
		break;
		/*
	case 1:
		// circle k=2
		for(i=0;i<swarmsize;i++){
			swarm[i].nlist[0] = (i-1) % swarmsize;
			swarm[i].nlist[1] = i;
			swarm[i].nlist[2] = (i+1) % swarmsize;
			if (swarm[i].nlist[0]<0)
				swarm[i].nlist[0] += swarmsize;
			if (swarm[i].nlist[1]<0)
				swarm[i].nlist[1] += swarmsize;
		}
		break;
	case 2:
		// particle 0 connected to all
		for(i=0;i<swarmsize;i++)
			swarm[0].nlist[i]=i;
		// all connected to particle 0
		for(i=1;i<swarmsize;i++){
			swarm[i].nlist[0]=0;
			swarm[i].nlist[1]=i;
		}
		break;
		*/
	}
	
	// best position and fitness
	windowPos = 0;
	assignfitn();
	windowPos = 0;	
	// initializing fitnesslist
	for(i=0;i<swarmsize+1;i++) {
		for(j=0;j<windowSize;j++) {
			swarm[i].fitnesslist[j] = DBL_MAX;
		} 
	}
	
	// global best position and fitness
	for(i=0;i<swarmsize;i++) {
		for(j=0;j<dim;j++)
			swarm[i].gbest[j] = swarm[i].pos[j];
		swarm[i].gbestfit = swarm[i].fitness;
	}
}

void updvel() {
	int i,j;
	w = (wstart-wend) * (tMax - evaluations) /tMax +wend;    //time variant weight, linear from weight to 0.4
	for(i=0;i<swarmsize;i++) {
		if (mainParticle(i)) {
			for(j=0;j<dim;j++) {
				swarm[i].vel[j] = w * swarm[i].vel[j] + dir * randval(0,phi1) * (swarm[i].best[j] - swarm[i].pos[j]) + dir * randval(0,phi2) * (swarm[i].gbest[j] - swarm[i].pos[j]);
				// checking for maxvelocity
				if (swarm[i].vel[j] > vmax)
					swarm[i].vel[j] = vmax;
				if (swarm[i].vel[j] < -vmax)
					swarm[i].vel[j] = -vmax;
			}
		}
	}
}

double getDiversity() {
  int i,j;
  double posStddev;
  double *avgpoint = malloc(dim * sizeof(double));
  // spredning paa partiklernes positioner
  for(i=0;i<dim;i++) {
    for(j=0;j<swarmsize;j++) {
      if (mainParticle(j))
	avgpoint[i] += swarm[j].pos[i];
    }
  }
  for(i=0;i<dim;i++)
    avgpoint[i] /= nParticles;
  for(i=0;i<swarmsize;i++) {
    if (mainParticle(i))
      posStddev += pow(dist2(avgpoint, swarm[i].pos),2);
  }
  posStddev /= MAX(nParticles - 1,1);
  posStddev = sqrt(posStddev);

  free(avgpoint);
  return posStddev;
}

void setPosAndVel(int a, int method) {
  int i;
  double fac3 = getDiversity();
  double fac2 = 1/(sqrt(t+1));
  double fac = (double)1-(double)evaluations/(double)tMax;
  for(i=0;i<dim;i++) {
    switch (method) {
    case 0:
      {