/* Swissknife PSO For future Standard PSO 2007, to replace PSO 2006 (Particle Swarm Central http://particleswarm.info) CFC: Call For Contribution Contact for remarks, suggestions etc.: MauriceClerc@WriteMe.com Contributors (last update 2007-02) James Kennedy and Russel Eberhart Maurice Clerc Manfred Stickel Abhi DattaSharma Tim Blackwell Renato Krohling Riccardo Poli William Langdon Hongbo Liu Vladimiro Miranda -------------------------------- Motivation Quite often some authors say they compare their PSO versions to the "standard one" ... which is never the same! So the idea is to define a real standard at least for one year. We propose here a PSO with a lot of options. You can suggest some others At the end of the process a lot of them will be removed, fo the final version should be a compromise between effectiveness and parsimony. Note that it does not intend to be the best one on the market (in particular there is no adaptation of the swarm size nor of the coefficients). Also this future "standard" will probably still accept a few options, like "clamping or not", or "fixed or random topology", if it appears that there is no consensus, and if results are equivalent. However this version will be kept as a "swiss knife" PSO, for future research. Last updates -------------------------------- 2007-02-27 new option for _when_ random topology is modified (topology=3) 2007-02-24 re-code the way random topology is designed. It is now faster. 2007-02-18 moves 21 and 22 2007-02-12 "jump+local search" model (move 20) 2007-01-31 Added some "adaptive" inertia weighs (moves 18 and 19) 2006-09-29 Partly rewritten so that it should be compatible with more compilers 2006-07-01 Option "no clamping, penalty value as fitness". Should be equivalent to "no clamping and no evaluation" 2006-06-20 Just for fun, added a (bad) fitness function that looks for n consecutive primes 2006-06-13 Cleaned up the code (slightly!) and added some comments 2006-05-31 Adaptive Bare Bones 2006-05-29 Gaussian kernels functions 2006-05-29 stop criterion option 2006-05-28 Option search space= list of values 2006-05-27 move option 16 (classical but with Gaussian distribution) 2006-04-23 Probability of evaluation 2006-04-18 Adaptive quantisation (?) 2006-04-15 Distorted Search Space option (?) 2006-04-12 move options 13,14,15 (?, ?, ?) 2006-04-11 Central Force option (?) 2006-04-11 DE-like (Differential Evolution) 2006-04-10 move option 9. Not really an improvement, even with a big K 2006-04-10 move option 8. (?) 2006-03-15 ---------------------------------- Compared to Standard PSO 2006 this code is more modular. It has been entirely rewritten so that it would be easier to translate it into Java, or any object language. On the other hand it is a bit slower, also because there are so many options, i.e. a lot of "if ..." -------------------------------- Parameters S := swarm size K := maximum number of particles _informed_ by a given one T := topology of the information links (info-network) w := first cognitive/confidence coefficient c := second cognitive/confidence coefficient R := random distribution of c B := rule "to keep the particle in the box" (clamping). It may be "do nothing" -------------------------------- Equations A) Classical For each particle and each dimension v(t+1) = w*v(t) + R(c)*(p(t)-x(t)) + R(c)*(g(t)-x(t)) x(t+1) = x(t) + v(t+1) where v(t) := velocity at time t x(t) := position at time t p(t) := best previous position of the particle g(t) := best previous position of the informants of the particle (not necessarily the whole swarm) B) Fully informed. Quite similar, but a R(c)*(p(t)-x(t)) term for each informant, not only for the best one. C) Bare-bones (without velocity) x(t+1)=D(p(t), g(t)) with D := random distribution "around" p(t) and g(t), for example Normal, Cauchy, or Levy. -------------------------------- Initialisation Initial positions are chosen at random inside the search space (which is supposed to be a hyperparallelepid, and even often a hypercube), according to an uniform distribution. This is not the best way, but the one of the original PSO, and quite simple. Some options for initial velocity. For example it can be defined as the difference of two random positions. The resulting distribution is a "triangular" one. Or it can be just set to zero. -------------------------------- Information links topology (info-network) A lot of works have been done about this topic. The main result is there is no best "fixed" topology, hence the random approach used in Standard PSO 2006. It seems it is still the best compromise, but some others have nevertheless been added. -------------------------------- Some results IMPORTANT As this is a beta version here are some results on classical test functions. So you can suggest any modification that improve them. The challenge is to find a algorithm that is good for the five functions below (0 to 4), and for any offset. This last point is very important. Too often authors test their algorithm almost only on functions whose optimum is near of the centre of the search space. In particular, "good" for function 4 (Tripod) would mean a success rate greater than 50% for any offset value. Search space Max. number of eval. Admissible error 0 Parabola/Sphere [-100,100]^30 9000 0.0001 1 Griewank [-100,100]^30 9000 0.0001 2 Rosenbrock [-10,10]^30 40000 0.0001 3 Rastrigin [-10,10]^30 40000 0.0001 4 Tripod [-100,100]^2 10000 0.0001 Offsets: 0 (0%), 0.5 (50%), 0.9 (90%), 1 (100%) Some good combinations A = best over the five functions and the three first offsets (not 100%) B = best over the five functions and the four offsets ......... A1 A2 A3 B1 C1 C2 A4 C3 velInit 0 0 0 0 0 0 0 0 AISS 0 0 0 2 0 0 0 0 clamping 0 0 0 0 1 1 1 1 move 0 0 0 0 18 19 0 21 topology 0 0 0 0 0 0 3 0 itself 1 1 1 1 1 1 1 1 K 3 S 3 3 3 3 3 1 probEval 1 1 0.5 1 1 1 1 1 Examples Success rate (or best mean value if success rate=0) over 50 runs. Note that the variance of the success rate is usually quite high, so you may easily find statistically equivalent other values With combination A1 offset 0% 50% 90% 100% Sphere 10% 8% 4% 0% Griewank 42% 26% 30% 24% Rosenbrock 31.1 26.94 25.04 22.98 Rastrigin 61.7 81.11 156.54 240.06 Tripod 44% 78% 100% 100% With A2 offset 0% 50% 90% 100% Sphere 92% 82% 72% 28% Griewank 16% 30% 24% 16% Rosenbrock 18.49 13.11 13.05 15.40 Rastrigin 86.58 159.35 377.1 494.28 Tripod 42% 64% 100% 100% With A3 offset 0% 50% 90% 100% Parabola 70% 58% 30% 0% Griewank 40% 34% 42% 16% Rosenbrock 41.96 34.69 26.89 25.22 Rastrigin 63.82 82.27 148.37 215.13 Tripod 40% 81% 100% 100% With C1 offset 0% 50% 90% 100% Parabola 82% 88% 60% 68% Griewank 36% 14% 24% 36% Rosenbrock 36.7 27.53 25.21 25.47 Rastrigin 59.50 68.31 113.86 143.33 Tripod 52% 88% 100% 100% With C2 offset 0% 50% 90% 100% Parabola 90% 84% 54% 6% Griewank 46% 44% 32% 24% Rosenbrock 34.9 28.9 23.8 24.3 Rastrigin 60.4 81.0 149.9 214.9 Tripod 44% 86% 100% 100% With A4 offset 0% 50% 90% 100% Parabola 30% Griewank 30% Rosenbrock 28.4 Rastrigin 61.1 Tripod 40% With C3 offset 0% 50% 90% 100% Parabola 88% 78% 78% 100% Griewank 52% 44% 48% 60% Rosenbrock 27.8 26.7 29 22.6 Rastrigin 182.1 180.3 116.4 32.9 Tripod 50% 24% 100% 100% ----------------------------------- Use See PROBLEM and PARAMETERS sections. Define the problem (you may add your own one in perf()) Choose your options Run and enjoy! ---------------------------------------------------------------------- */ #include "stdio.h" #include "math.h" #include #include #define D_max 110 // Max number of dimensions of the search space #define R_max 40000 // Max number of runs #define S_max 100 // Max swarm size #define AV_max 20 // Max of admissible positions on each dimension // when the search space is a list of such positions #define zero 0 // 1.0e-30 // To avoid numerical instabilities // Structures struct quantum { double q[D_max]; int size; }; struct SS { int D; double max[D_max]; double min[D_max]; struct quantum q; // Quantisation step size. 0 => continuous problem int avSize[D_max]; // If > 0, the search space is defined by a list // of admissible values, at most AV_max for each dimension double av[D_max][AV_max]; // Admissible values }; struct param { int AISS; // Kind of Adaptive Local Search Space to use double alpha; // For central option, for chaotic option double beta; // For central force option double c; // Confidence coefficient int CI; // Flag for Constant Informed PSO int clamping; // Flag for Position clamping double distort; // For distorted search space option double Dt; // For chaotic option int itself; // For neighbourhood double jump; // Parameter to define the probability of a jump int K; // Max number of particles informed by a given one int local1; // To estimate the number of normal searches (phase 1) int local2; // To estimate the number of local searches (phase 2) double localRadius; // Relative radius for local search int mode; // Parallel or sequential int move; // Move option (standard, deterministic, etc) double p; // Probability threshold for random topology double probEval; // Probability of evaluation double qScale1; // For adaptive quantisation double qScale2; // For adaptive quantisation int quantis; // Flag for adaptive quantisation int randChoice; // Flag for random choice or not of particles int S; // Swarm size int saveEnergy; // Flag for saving swarm energy on a file int stop; // Flag for stop criterion int topology; // Neighbourhood topology int traj; // For testing beta version: rank of the // particle whose trajectory is saved int velInit; // Kind of velocity initialisation int velReInit; // Flag for velocity reinitialisation double w; // Confidence coefficient double w1; // For adaptive versions double w2; // For adaptive versions }; struct position { double f; int improved; int size; double x[D_max]; }; struct velocity { int size; double v[D_max]; }; struct problem { double epsilon; int evalMax; int function; double objective; double offset[D_max]; struct SS SS; // Search space }; struct swarm { int best; struct position P[S_max]; // Previous best positions found by each particle struct position PP[S_max]; // Previous previous bests int S; struct velocity V[S_max]; // Velocities struct position X[S_max + 1]; // Positions + queen struct position XP[S_max + 1]; // Previous positions. // If used, V is not necessary, for V=X-XP int worst; }; struct result { double nEval; struct swarm SW; }; // Sub-programs double alea (double a, double b); double alea_cauchy (double mean, double gamma, double tMin, double tMax); int alea_integer (int a, int b); double alea_levy (double mean, double scale, double alpha); double alea_normal (double mean, double stdev); struct velocity alea_sphere(int D, double radius); double alea_triangle (double mean, double width); struct position clampToList (struct position x, struct SS SS); double gauss (double x, double mu, double var); double energy (struct swarm SW); double energyS (struct position x, struct position xp); double distanceL (struct position x1, struct position x2, double L); struct position distort (struct position x, struct SS, int way, double alpha); double distortX (double x, double min, double max, double alpha); double max(double a,double b); double min(double a, double b); double normL (struct velocity,double L); double perf (struct position x, int function, double offset[], double objective, struct SS SS, double alpha); // Fitness evaluation struct quantum quantAdapt (struct quantum quant, struct quantum q, double scale); struct position quantis (struct position x, struct SS SS, struct quantum quants, int option); struct position queen (struct swarm SW, struct problem pb, struct SS SS, double distort); // Just for tests struct result PSO (struct SS SS, struct param param, struct problem problem); int sign (double x); void traj (struct result R, int rank); // Global variables long double E; // exp(1). Useful for some test functions long double pi; // Useful for some test functions // File(s); FILE * f_run; FILE * f_stag; FILE * f_synth; FILE * f_traj; FILE * f_energy; // ================================================= int main () { struct position bestResult; int d; // Current dimension double error; // Current error double errorMean; // Average error double errorMin; // Best result through several runs double errorMeanBest[R_max]; double evalMean; // Mean number of evaluations int nFailure; // Number of unsuccessful runs struct param param; struct problem pb; int run, runMax; struct result result; double shift, shiftMin, shiftMax, shiftDelta; int shiftType; struct SS SS; // Search space double successRate; double variance; f_energy = fopen ("f_energy.txt", "w"); f_run = fopen ("f_run.txt", "w"); f_stag = fopen ("f_stag.txt", "w"); f_synth = fopen ("f_synth.txt", "w"); f_traj = fopen ("f_traj.txt", "w"); E = exp ((long double) 1); pi = acos ((long double) -1); // ----------------------------------------------- PROBLEM pb.function = 2; // Function code /* 0 Parabola (Sphere) 1 Griewank 2 Rosenbrock (Banana) 3 Rastrigin 4 Tripod (dimension 2) 99 Tests (see perf() ) */ // ------------ pb.epsilon = 0.0001; // Acceptable error pb.objective = 0; // Objective value // Define the loop on offset values shiftType = 0; // Kind of offset sequence: // 0 => deterministic (diagonal) // 1 => random in the search range shiftMin = 0; // Relative offset [0,1]. Difficult value: 1 / sqrt( 3 ); shiftMax = shiftMin; shiftDelta = 0.1; // Keep it positive /* If shiftType=0, shiftMin and shiftMax are relative offsets in [-1,1]. The real offset will be shift*(xmax-xmin)/2 If shiftType=1, it just mean there are L loops of at most runMax runs (see below) with L=(shiftMax-shiftMin)/shiftDelta */ runMax = 10; // Numbers of runs if (runMax > R_max) runMax = R_max; // ------------------ Search space switch (pb.function) { case 0: // Parabola SS.D = 30; // Dimension for (d = 0; d < SS.D; d++) SS.avSize[d] = 0; // NOT a list of admissible // values for (d = 0; d < SS.D; d++) { SS.min[d] = -100; SS.max[d] = 100; SS.q.q[d] = 0; // E / 10; // Relative quantisation, in [0,1]. // Worst case for tests: transcendant number // The real one is q[d]*(max[d]-min[d])/2 } pb.evalMax = 9000; // Max number of evaluations for each run break; case 1: // Griewank SS.D = 30; // Dimension for (d = 0; d < SS.D; d++) SS.avSize[d] = 0; // Boundaries for (d = 0; d < SS.D; d++) { SS.min[d] = -100; SS.max[d] = 100; SS.q.q[d] = 0; // SS.min[d] = -600; SS.max[d] = 600; SS.q.q[d] = 0; } pb.evalMax = 9000; break; case 2: // Rosenbrock SS.D = 30; // Dimension for (d = 0; d < SS.D; d++) SS.avSize[d] = 0; // Boundaries for (d = 0; d < SS.D; d++) { SS.min[d] = -10; SS.max[d] = 10; SS.q.q[d] = 0; } pb.evalMax = 40000 ; break; case 3: // Rastrigin SS.D = 30; // Dimension for (d = 0; d < SS.D; d++) SS.avSize[d] = 0; // Boundaries for (d = 0; d < SS.D; d++) { SS.min[d] = -10; SS.max[d] = 10; SS.q.q[d] = 0; // E / 10; } pb.evalMax = 40000; break; case 4: // Tripod SS.D = 2; // Dimension for (d = 0; d < SS.D; d++) SS.avSize[d] = 0; // Boundaries for (d = 0; d < SS.D; d++) { SS.min[d] = -100; SS.max[d] = 100; SS.q.q[d] = 0; // E/10; } pb.evalMax = 10000; //10000; break; } SS.q.size = SS.D; // ----------------------------------------------------- // PARAMETERS // * means "seems quite good" param.AISS = 0; // AISS for current iteration // 0 => no AISS. ** // 1 => possible for any particle // 2 => possible just once by neighbourhood. * param.alpha=2; // For Central Force move option param.beta = 1; // For Central Force move option: // distance^beta // Note 1: you should use full-connected topology, or, better FIPS // Note 2: for "classical" central force, use a negative beta //param.c; // See param.w param.CI = 0; // For Constant Informed PSO // 0 => don't use CI option * // 1 => use PP if the particle has improved its position, P else // 2 => use P, and PP param.clamping =1; // -1 => no clamping. WARNING: the function has to be defined everywhere // 0 => no clamping AND no evaluation. WARNING: the program // may never stop (in particular with option move 20 (jumps)) * // 1 => classical. Set to bounds, and velocity to zero // 2 => 0 or 1 with probability 0.5 // 3 no clamping, but set a penalty as fitness value param.distort = 1; // Distorted Search Space // 1 => no distortion * // <1 => contraction // >1 => expansion param.Dt = 1; // Only for Central Force move option. Time step param.itself = 1; // 0 => a particle doesn't belong to its neighbourhood. It // can't be its own best neighbour . // 1 => a particle belongs to its neighbourhood. To be // consistent, you should choose this value when using FIPS. * //param.K // See after param.S // Only for move option(s) using local search (20) param.local1= 1; // Number of consecutive normal searches (phase 1) param.local2= 20; // Number of consecutive local searches (phase 2) param.localRadius=0.2; // Relative radius for local search /* local1 local2 localRadius => (for offset 0%) Sphere 20 10 0.3 92% Griewank 20 10 0.3 60% Rosenbrock 1 2 0.2 27.0 Rastrigin 20 10 0.3 59.2 Tripod 1 2 0.2 82% */ param.mode = 0; // 0 => sequential mode ("classical" PSO). * // 1 => parallel mode *** TO COMPLETE *** ... or maybe not // Up to now any parallel version is always a bit less effective // than the same in sequential mode param.move = 0; // * means "seems quite good" // ? means "really questionable" // 0 *=> classical (uniform distribution). // 1 => using the previous previous best // 2 => classical, but deterministic // 3 *=> Gauss Bare-bones PSO (weighted). // 4 => FIPS (Fully informed PSO), deterministic // 5 => FIPS, random (uniform). // 6 => Cauchy Bare-bones PSO // 7 ?=> Lévy Bare-bones PSO // 8 ?=> two nearest informants on each dimension + random // 9 => two best informants // 10 => DE-like. No velocity. K is not used. S should be quite high // (at least 5*D) // 11 ? => FIPS Central force, deterministic. TO MODIFY // 12 ? => Central force. TO MODIFY // 13 => Stupid Gauss Bare-bones PSO // 14 => Towards two promising places + velocity // 15 => Towards a promising place + velocity // 16 => classical (gaussian distribution) // 17 => adaptive Bare Bones // 18 *=> adaptive inertia weight formula w'=w*f(gbest)/f(pbest) // and adaptive c=(w'+1)^2/2 // Note: assuming f > 0 // 19 => adaptive inertia weight formula // and adaptive c= formula coming from Stagnation Analysis // 20 => normal + local searches // 21 => FIPS according to fitnesses // 22 => Normal + imprecise information //param.p; // See after param.S param.probEval = 1; // Probability of evaluation (1 for standard PSO) // 0.5 seems also interesting param.qScale1 = 2; // Scaling quantisation step if // improvement. >1 => bigger step // Active only if "quantis">0 param.qScale2 = pi / 4; // Scaling quantisation step if no improvement. // <1 <=> smaller step // Active only if "quantis">0 param.quantis = 0; // How to take q(d) into account in search space definition // 0 => constant quantisation (see search space definition of the // current problem) // 1 => global adaptive quantisation, according to global best status // 2 => individual adaptive quantisation * // 3 => global, according to all best statuses, and to stag/prog sequences param.randChoice=0; // 0* => at each iteration, particles are modified // always according to the same order 0..S-1 // 1 => random order // 2 => TO DO. at each iteration, sort the particles by fitness // and perform the next loop according to this order param.S = (int) (10 + 2 * sqrt ((double) SS.D)); // Swarm size if (param.S > S_max) param.S = S_max; param.K=3; // 3*." // For some move options (9), you need to increase it, however // (say K=S) // Probability threshold for informants param.p=1-pow(1-(double)1/(param.S),param.K); //p=min(param.K,param.S)/param.S; // Approximation param.saveEnergy = 0; // Just for information // 0 => do nothing special // 1 => save swarm energy on a file param.stop = 0; // Stop criterion // 0 => error < pb.epsilon // 1 => eval < pb.evalMax param.topology = 0; // 0 => random, modified if there has been no improvement. * // 1 => fixed, Star (full-connected topology) // 2 => fixed, Ring // 3 => random, modified after a number of iterations // depending on the current topology param.traj = -1; // For test. Rank of the particle whose trajectory is saved // -1 => no one param.velInit = 0; // 0 => initial velocity set to zero // 1 => "triangular" random distribution param.velReInit = 0; // 0 => no velocity reinitialisation after re-quantisation // >0 => reinitialisation after re-quantisation param.w = 1 / (2 * log ((double) 2)); param.c = 0.5 + log ((double) 2); param.w1 = 0.74; param.w2 = 0.74; // Note that for Weighted Gaussian Bare-Bones there are good theoretical reasons // to choose w1=0.74 // (cf. http://clerc.maurice.free.fr/pso/Ideas.pdf) // Note also that the same kind of analysis for Cauchy Bare-Bones suggests // that the best "weight" should be just 1. // ------------------------------------------------------- // Some information printf ("\n Function %i ", pb.function); printf ("\n Swarm size %i", param.S); // =========================================================== // RUNs // --------- Loop on offset for (shift = shiftMin; shift <= shiftMax;shift = shift + shiftDelta) { printf ("\n-----------"); printf ("\nShift %f", shift); switch (shiftType) { case 0: // Deterministic, on diagonal for (d = 0; d < SS.D; d++) { pb.offset[d] = shift * (SS.max[d] - SS.min[d]) / 2; } break; case 1: // Random for (d = 0; d < SS.D; d++) { pb.offset[d] = alea (-1, 1) * (SS.max[d] - SS.min[d]) / 2; } break; } errorMean = 0; evalMean = 0; nFailure = 0; for (run = 0; run < runMax; run++) { srand (clock () / 100); // May improve pseudo-randomness result = PSO (SS, param, pb); error = result.SW.P[result.SW.best].f; if (error > pb.epsilon) nFailure = nFailure + 1; // Result display printf ("\nRun %i. Eval %f. Error %f ", run, result.nEval, error); bestResult = distort (result.SW.P[result.SW.best], SS, -1, param.distort); printf ("\n Position :\n"); for (d = 0; d < SS.D; d++) printf (" %f", bestResult.x[d]); // Save result // fprintf( f_run, "\n%i %.0f %f ", run, result.nEval, error ); // for ( d = 0; d < SS.D; d++ ) fprintf( f_run, " %f", bestResult.x[d] ); // Compute/store some statistical information if (run == 0) errorMin = error; else if (error < errorMin) errorMin = error; evalMean = evalMean + result.nEval; errorMean = errorMean + error; errorMeanBest[run] = error; } // End loop on "run" // END. Display some statistical information evalMean = evalMean / (double) run; errorMean = errorMean / (double) run; printf ("\n Eval. (mean)= %f", evalMean); printf ("\n Error (mean) = %f", errorMean); // Variance variance = 0; for (run = 0; run < runMax; run++) variance = variance + pow (errorMeanBest[run] - errorMean, 2); variance = sqrt (variance / runMax); printf ("\n Std. dev. %f", variance); // Success rate and minimum value successRate = 100 * (1 - nFailure / (double) run); printf ("\n Success rate = %.2f%%", successRate); if (runMax > 1) printf ("\n Best min value = %f", errorMin); // Save /* fprintf(f_synth,"\n"); for (d=0;d normal search to begin nbLocalSearchMax=0; nbLocalSearch=param.local2; phase=1; // Phase to begin with // Take quantisation into account R.SW.X[s] = quantis (R.SW.X[s], ISS, quant[s], 0); // If list, clamp to the nearest admissible value R.SW.X[s] = clampToList (R.SW.X[s], ISS); } // First evaluations for (s = 0; s < R.SW.S; s++) { R.SW.X[s].f = perf (R.SW.X[s], pb.function, pb.offset, pb.objective, SS, param.distort); R.SW.P[s] = R.SW.X[s]; // Best position = current one R.SW.P[s].improved = 0; // No improvement R.SW.PP[s] = R.SW.P[s]; } R.nEval = R.SW.S; // Find the best and the worst R.SW.best = 0; R.SW.worst = 0; for (s = 1; s < R.SW.S; s++) { if (R.SW.P[s].f < R.SW.P[R.SW.best].f) R.SW.best = s; if (R.SW.P[s].f > R.SW.P[R.SW.worst].f) R.SW.worst = s; } errorPrev = R.SW.P[R.SW.best].f; penalty = R.SW.P[R.SW.worst].f; // Display the best /* printf( " Best value after init. %f ", errorPrev ); printf( "\n Position :\n" ); for ( d = 0; d < SS.D; d++ ) printf( " %f", R.SW.P[R.SW.best].x[d] ); */ initLinks = 1; // So that information links will beinitialized // Note: It is also a flag saying "No improvement" noStop = 0; stagSeq = 0; // Just for information. Number ofstagnation phases progSeq = 0; // Just for information. Number of xs.xphases progNb = 0; stagNb = 0; w=param.w; switch (param.move) { case 3: // Gauss Bare Bones w = param.w1; break; } switch (param.topology) { case 0: // Random break; case 1: // Fully connected for (s = 0; s < R.SW.S; s++) { for (m = 0; m < R.SW.S; m++) LINKS[m][s] = 1; // Init to "no link" if (param.itself == 0) LINKS[s][s] = 0; } break; case 2: // Ring printf ("\n Ring topology, K=%i", param.K); for (s = 0; s < R.SW.S; s++) { for (m = 0; m < R.SW.S; m++) LINKS[m][s] = 0; // Init "no link" } for (s = 0; s < R.SW.S; s++) { for (i = 0; i < param.K; i++) { m = s - param.K / 2 + i; if (m < 0) m = m + param.S; if (m >= param.S) m = m - param.S; LINKS[m][s] = 1; } if (param.itself == 0) LINKS[s][s] = 0; } break; } LinksTot=0; // Total number of links for (s = 0; s < R.SW.S; s++) { for (m = 0; m < R.SW.S; m++) if (LINKS[m][s] == 1) LinksTot=LinksTot+1; } // ---------------------------------------------- ITERATIONS iter=0; iterBegin=0; while (noStop == 0) { iter=iter+1; // Queen (for tests) R.SW.XP[queenRank] = R.SW.X[queenRank]; R.SW.X[queenRank] = queen (R.SW, pb, SS, param.distort); if (param.saveEnergy > 0) { // temp=energy(R.SW); // Swarm temp = energyS (R.SW.X[queenRank], R.SW.XP[queenRank]) * R.SW.S; // Queen fprintf (f_energy, "%.0f %f \n", R.nEval, temp); } if (param.traj >= 0) // Save trajectory of a particle { traj (R, param.traj); } randTopology=0; switch (param.topology) { case 0: if (initLinks == 1) randTopology=1; break; case 3: if (iter-iterBegin>LinksTot/2) { iterBegin=0; randTopology=1; } } if (randTopology==1) // Random topology { // Who informs who, at random for (s = 0; s < R.SW.S; s++) { for (m = 0; m < R.SW.S; m++) { temp = alea (0, 1); if (temp= R.SW.S) s1 = s; switch (param.move) // Informants { default: // Find the best informant g = s1; // according to best previous gp = s1; // according to best previous previous for (m = s1; m < R.SW.S; m++) { if (LINKS[m][s] == 1 && R.SW.P[m].f < R.SW.P[g].f) g = m; if (LINKS[m][s] == 1 && R.SW.PP[m].f < R.SW.PP[gp].f) gp = m; } break; case 8: // The two nearest ones, on each dimension for (k = 0; k < 2; k++) { for (d = 0; d < SS.D; d++) { distN = SS.max[d] - SS.min[d]; for (s1 = 0; s1 < R.SW.S; s1++) { if (k == 1 && s1 == informant[d][0]) continue; dist = fabs (R.SW.X[s].x[d] - R.SW.P[s1].x[d]); if (dist < distN) { sN = s1; distN = dist; } } informant[d][k] = sN; } // printf("\nk=%i, %i %i ",k,informant[0][k],informant[1] [k]); } break; case 9: // The two best informants g1 = s1; for (m = s1; m < R.SW.S; m++) { if (LINKS[m][s] == 1 && R.SW.P[m].f < R.SW.P[g1].f) g1 = m; } g2 = s1; for (m = s1; m < R.SW.S; m++) { if (m == g1) continue; if (LINKS[m][s] == 1 && R.SW.P[m].f < R.SW.P[g2].f) g2 = m; } break; case 4: // FIPS case 5: // FIPS case 11: // FIPS case 13: // Stupid Bare-Bones. No need of g case 14: // FIPS case 21: // FIPS break; } // ... Adaptive Individual Search Space switch (param.AISS) { case 0: // No AISS AISS = 1 == 0; // Always false break; case 1: // True, if no improvement AISS = R.SW.P[s].improved <= 0; break; case 2: // True if no improvement AND the first // time in g's neighbourhood AISS = R.SW.P[s].improved <= 0 && nAISS[g] == 0; if (AISS) nAISS[g]++; break; } for (d = 0; d < SS.D; d++) { if (AISS) { ps = R.SW.P[g].x[d]; xM = 2 * ps - R.SW.X[s].x[d]; if (xM >= ps) { ISS.min[d] = R.SW.X[s].x[d]; ISS.max[d] = xM; if (ISS.max[d] > SS.max[d]) ISS.max[d] = SS.max[d]; } else { ISS.min[d] = xM; if (ISS.min[d] < SS.min[d]) ISS.min[d] = SS.min[d]; ISS.max[d] = R.SW.X[s].x[d]; } } // else { xMin[d] = ISS.min[d]; xMax[d] = ISS.max[d]; } } for (d = 0; d < SS.D; d++) { Vmax.v[d]=ISS.max[d]-ISS.min[d]; } diameter=normL(Vmax,param.alpha); // ... compute the new velocity, and move noMove = 0; wTemp=w; cTemp=param.c; switch (param.move) // Define the two informants { default: // Classical switch (param.CI) // Constant Informed option { default: // classical p and g ( => NOT constant informed) P1 = R.SW.P[s]; P2 = R.SW.P[g]; break; case 1: // p ang if improvement, previous_p and g else if (R.SW.P[s].improved < 0) { P1 = R.SW.P[s]; P2 = R.SW.P[g]; } else { P1 = R.SW.PP[s]; P2 = R.SW.PP[g]; // P1 = R.SW.PP[s]; P2 = R.SW.PP[gp]; } break; } } switch (param.move) // Some adaptations for some move options { default: break; case 18: wTemp=w*P2.f/P1.f; cTemp=pow(wTemp+1,2)/2; break; case 19: wTemp=w*P2.f/P1.f; cTemp=(beta*wTemp+delta)*0.5+ sqrt(wTemp*(alpha*wTemp+gammac)+ eta+pow(beta*wTemp+delta,2)*0.25); cTemp=cTemp*0.5; break; } switch (param.move) // Move { default: switch (param.CI) { default: temp = 1; // <=1 for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d]=wTemp *R.SW.V[s].v[d] + alea(0, temp * cTemp) * (P1.x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d]=R.SW.V[s].v[d] + alea(0,(2 - temp) * cTemp) * (P2.x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 2: temp = -0.5 + sqrt ((double) 5) / 2; for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d]=param.w *R.SW.V[s].v[d] + alea(0, temp * param.c) * (P1.x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d]=R.SW.V[s].v[d] + alea (0, param.c) * (P2.x[d] - R.SW.X[s].x[d]); /* R.SW.V[s].v[d] = R.SW.V[s].v[d] + alea( 0,temp ) * ( R.SW.PP[s].x[d] - R.SW.X[s].x[d] ); */ R.SW.V[s].v[d]=R.SW.V[s].v[d] +alea(0, temp*temp* param.c) * (R.SW.PP[s].x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; } break; case 1: // Using previous previous best for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = param.w *R.SW.V[s].v[d] + alea (0,param.c) *(R.SW.P[s].x[d] - R.SW.X[s].x[d]); temp = alea (0, 1) * (R.SW.P[g].x[d] + R.SW.PP[g].x[d]); R.SW.V[s].v[d] =R.SW.V[s].v[d] + alea (0, param.c) * (temp - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 2: // Deterministic for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] =param.w *R.SW.V[s].v[d] + 0.5 * param.c * (R.SW.P[s].x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d] = R.SW.V[s].v[d] + 0.5 * param.c * (R.SW.P[g].x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 17: // Adaptive Gaussian Bare Bones w = param.w1 + (param.w2 -param.w1) * (0.5 -atan (0.5 *R.SW.P[s].improved / param.S) / pi); // printf("\n%f",w); // NO break here case 3: // Gauss Bare Bones for (d = 0; d < SS.D; d++) { mean = (R.SW.P[s].x[d] + R.SW.P[g].x[d]) / 2; stdev = w * fabs (R.SW.P[s].x[d] - R.SW.P[g].x[d]); R.SW.X[s].x[d] = alea_normal (mean, stdev); } break; case 4: // FIPS, deterministic neigh = 0; // Number of neighbours for (m = 0; m < param.S; m++) if (LINKS[m][s] == 1) neigh++; temp = 2 * param.c / neigh; for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = param.w * R.SW.V[s].v[d]; for (m = 0; m < param.S; m++) { if (LINKS[m][s] == 1) R.SW.V[s].v[d] = R.SW.V[s].v[d] + temp * (R.SW.P[m].x[d] - R.SW.X[s].x[d]); } R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 5: // FIPS, random neigh = 0; // Number of neighbours for (m = 0; m < param.S; m++) { if (LINKS[m][s] == 1) neigh++; } temp = 2 * param.c / neigh; for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = param.w * R.SW.V[s].v[d]; for (m = 0; m < param.S; m++) { if (LINKS[m][s] == 1) R.SW.V[s].v[d] = R.SW.V[s].v[d] + alea (0,temp) * (R.SW.P[m].x[d] -R.SW.X[s].x[d]); } R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 6: // Cauchy Bare Bones for (d = 0; d < SS.D; d++) { mean = (R.SW.P[s].x[d] + R.SW.P[g].x[d]) / 2; gamma = fabs (R.SW.P[s].x[d] - R.SW.P[g].x[d]) / 2; // So that the new position will be inside the search space tMin = -atan (fabs (mean - gamma) / gamma); temp = -atan (fabs (mean - SS.min[d]) / gamma); if (temp > tMin) tMin = temp; tMax = atan (fabs (mean + gamma) / gamma); temp = atan (fabs (mean - SS.max[d]) / gamma); if (temp < tMax) tMax = temp; R.SW.X[s].x[d] = alea_cauchy (mean, gamma, tMin, tMax); } break; case 7: // Lévy Bare Bones for (d = 0; d < SS.D; d++) { mean = (R.SW.P[s].x[d] + R.SW.P[g].x[d]) / 2; stdev = fabs (R.SW.P[s].x[d] - R.SW.P[g].x[d]) / 2; R.SW.X[s].x[d] = alea_levy (mean, stdev, 1); } break; case 8: // two nearest informants on each dimension for (d = 0; d < SS.D; d++) { x1 = R.SW.P[informant[d][0]].x[d]; x2 = R.SW.P[informant[d][1]].x[d]; // temp = alea_triangle((x1+x2)/2 , fabs(x1-x2)); // R.SW.X[s].x[d]= temp+ // param.w*(R.SW.X[s].x[d]-R.SW.X[s].x[d]); // temp=alea((x1+x2)/2, fabs(x1-x2)); temp = alea (0, param.c) * (x1 -R.SW.X[s].x[d]) + alea (0, param.c) * (x2 - R.SW.X[s].x[d]); // temp=alea(0,param.c)*(x1-R.SW.X[s].x[d]); R.SW.V[s].v[d] = param.w * R.SW.V[s].v[d] + temp; R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 9: // Using the two best informants for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = param.w * R.SW.V[s].v[d] + alea (0,param.c) *(R.SW.P[g1].x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d] = R.SW.V[s].v[d] + alea (0, param.c) * (R.SW.P[g2].x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 10: // DE-like. No velocity g = alea_integer (0, param.S - 1); g1 = alea_integer (0, param.S - 1); g2 = alea_integer (0, param.S - 1); for (d = 0; d < SS.D; d++) { /* temp = alea( 0, 1 ); R.SW.X[s].x[d] = temp * R.SW.P[s].x[d] + ( 1 -temp ) * R.SW.P[g].x[d] + alea( 0, 2 ) * ( R.SW.P[g1].x[d] - R.SW.P[g2].x[d] ); */ temp = R.SW.P[g].x[d] + alea (0,2) * (R.SW.P[g1].x[d] - R.SW.P[g2].x[d]); if (alea (0, 1) < alea (0, 2)) R.SW.X[s].x[d] = temp; } break; case 11: // FIPS + Central force neigh = 0; // Number of neighbours centralTot=0; // Compute the weights for (m = 0; m < param.S; m++) { if (LINKS[m][s] == 1 && R.SW.P[m].f<=R.SW.X[s].f ) { neigh++; temp = pow(distanceL (R.SW.X[s], R.SW.P[m],2),param.beta); central[m]=temp*exp(-R.SW.P[m].f); centralTot=centralTot+central[m]; } else central[m]=0; } // Normalize coefficients for (m = 0; m < param.S; m++) { central[m]=central[m]/centralTot; } deltaT = param.Dt; for (d = 0; d < SS.D; d++) { //R.SW.V[s].v[d]=0; R.SW.V[s].v[d]=param.w*R.SW.V[s].v[d]; for (m = 0; m < param.S; m++) { temp=2*param.c*central[m]; if (LINKS[m][s] == 1 && R.SW.P[m].f<=R.SW.X[s].f ) { if (R.SW.P[m].x[d] - R.SW.X[s].x[d]>=0) sign=1; else sign=-1; R.SW.V[s].v[d] = R.SW.V[s].v[d] + temp * deltaT * sign; } /* if (LINKS[m][s] == 1 && R.SW.X[m].f<=R.SW.X[s].f ) { //temp=2*(R.SW.X[s].f-R.SW.X[m].f)/(R.SW.X[s].f+R.SW.X[m].f); temp=R.SW.X[s].f-R.SW.X[m].f; temp=pow(temp,param.alpha); R.SW.V[s].v[d] = R.SW.V[s].v[d] + central[m] * deltaT * temp * (R.SW.X[m].x[d] - R.SW.X[s].x[d]); } */ } R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d] * deltaT; } break; case 12: // Central force temp=distanceL (R.SW.X[s], R.SW.P[s],param.alpha)/diameter; /* if (temp>2) central[s] = 1/temp; else */ central[s]=1; central[s]=central[s]*pow(param.w+1,2)/2; temp=distanceL (R.SW.X[s], R.SW.P[g],param.alpha)/diameter; /* if (temp>2) central[g] = 1/temp; else */ central[g]=1; central[g]=central[g]*pow(param.w+1,2)/2; central[s]=cTemp; central[g]=cTemp; printf("\n%f %f %f",wTemp,central[s],central[g]); deltaT = param.Dt; for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = wTemp * R.SW.V[s].v[d] + alea(0,central[s]) * deltaT *(P1.x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d] = R.SW.V[s].v[d] + alea(0,central[g]) * deltaT * (P2.x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d] * deltaT; break; } case 13: // Pure normal distribution (no velocity). // Weighted Bare-Bones for (d = 0; d < SS.D; d++) { stdev = param.w * fabs (R.SW.P[s].x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = alea_normal (R.SW.P[s].x[d], stdev); } break; case 14: // Towards two promising places + velocity g = alea_integer (0, param.S - 1); if (R.SW.P[g].f > R.SW.P[s].f) { noMove = 1; break; } for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = param.w *R.SW.V[s].v[d] + alea (0,param.c) * (R.SW.P[s].x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d] = R.SW.V[s].v[d] + alea (0,param.c) * (R.SW.P[g].x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 15: // Towards a promising place + velocity g = alea_integer (0, param.S - 1); if (R.SW.P[g].f > R.SW.P[s].f) { noMove = 1; break; } for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = param.w * R.SW.V[s].v[d] + alea (0,param.c) * (R.SW.P[g].x[d] - R.SW.X[s].x[d]); /* stdev = param.w*fabs( R.SW.P[g].x[d] - R.SW.X[s].x[d] ); R.SW.V[s].v[d] = param.w * R.SW.V[s].v[d] + alea_normal(R.SW.P[g].x[d] , stdev ); */ R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 16: // With Gaussian distribution (cf. Renato Krohling) for (d = 0; d < SS.D; d++) { R.SW.V[s].v[d] = fabs(alea_normal (0, 1)) * (R.SW.P[s].x[d] - R.SW.X[s].x[d]); R.SW.V[s].v[d] = R.SW.V[s].v[d] +fabs(alea_normal(0, 1)) * (R.SW.P[g].x[d] - R.SW.X[s].x[d]); R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } break; case 20: // If too many iterations without improvement, then "jump" //if (-R.SW.P[s].improved>param.jump) if (phase==1 && nbNormalSearch ISS.max[d]) outside++; } // printf("\n d%i outside %i",d,outside); if (outside == 0) // If inside, the position is evaluated { if (noEval == 0) { R.SW.X[s].f = perf (R.SW.X[s], pb.function, pb.offset, pb.objective, SS, param.distort); R.nEval = R.nEval + 1; } } break; case 1: // Set to the bounds, and v to zero clamping1: for (d = 0; d < SS.D; d++) { if (R.SW.X[s].x[d] < ISS.min[d]) { R.SW.X[s].x[d] = ISS.min[d]; R.SW.V[s].v[d] = 0; } if (R.SW.X[s].x[d] > ISS.max[d]) { R.SW.X[s].x[d] = ISS.max[d]; R.SW.V[s].v[d] = 0; } } if (noEval == 0) { R.SW.X[s].f = perf(R.SW.X[s],pb.function, pb.offset, pb.objective, SS, param.distort); R.nEval = R.nEval + 1; } break; case 2: // 0 or 1 with probability 0.5 if (alea (0, 1) < 0.5) goto clamping0; else goto clamping1; case 3: // no clamping and use a penalty as fitness value outside = 0; for (d = 0; d < SS.D; d++) { if (R.SW.X[s].x[d] < ISS.min[d] || R.SW.X[s].x[d] > ISS.max[d]) outside++; } if (outside > 0) // If outside, the position is artifically evaluated { R.SW.X[s].f = penalty; } else if (noEval == 0) { R.SW.X[s].f = perf(R.SW.X[s],pb.function, pb.offset, pb.objective, SS, param.distort); R.nEval = R.nEval + 1; } break; } } else // AISS (no need of clamping) { if (noEval == 0) { R.SW.X[s].f = perf (R.SW.X[s],pb.function, pb.offset, pb.objective, SS, param.distort); R.nEval = R.nEval + 1; } } // ... update the best previous position if (R.SW.X[s].f < R.SW.P[s].f) // Improvement { R.SW.PP[s] = R.SW.P[s]; // Update previous previous best R.SW.P[s] = R.SW.X[s]; // Again improvement if (R.SW.P[s].improved > 0) R.SW.P[s].improved = R.SW.P[s].improved + 1; else // Improvement after stagnation { R.SW.P[s].improved = 1; } nAISS[s] = 0; // Re-initialise the number of particles // that have use it for AISS // ... update the best of the bests if (R.SW.P[s].f < R.SW.P[R.SW.best].f) { R.SW.best = s; } } else // (again) no improvement { if (R.SW.P[s].improved < 0) R.SW.P[s].improved = R.SW.P[s].improved - 1; else { R.SW.P[s].improved = -1; // Stagnation after improvement } } } // End of "for (s=0 ... " // Check if finished error = R.SW.P[R.SW.best].f; if (error < errorPrev) // Improvement { initLinks = 0; progNb++; // Progress length if (progNb == 1) // If previous phase was stagnation { // Save stagnation info, for future tests // stag[stagSeq] = -stagNb; stagSeq++; // fprintf( f_stag, "%i \n", -stagNb ); } stagNb = 0; } else // No improvement { initLinks = 1; // Information links will be reinitialized stagNb++; // Stagnation length if (stagNb == 1) // If previous phase was improvement { // Save stagnation info, for future tests // prog[progSeq] = progNb; progSeq++; // fprintf( f_stag, "%i \n", progNb ); } progNb = 0; } // Adapt quantisation for (s = 0; s < R.SW.S; s++) quantUpdate[s] = 0; switch (param.quantis) { case 1: // Global, according to global best status if (initLinks == 1) // If no improvement { ISS.q = quantAdapt (ISS.q, ISS.q, param.qScale2); for (s = 0; s < R.SW.S; s++) quantUpdate[s] = 1; } else // If improvement { ISS.q = quantAdapt (ISS.q, ISS.q, param.qScale1); for (s = 0; s < R.SW.S; s++) quantUpdate[s] = 1; } break; case 2: // Individual for (s = 0; s < R.SW.S; s++) { if (R.SW.P[s].improved > 0) // If progression { quant[s] = quantAdapt (quant[s], ISS.q, param.qScale1); quantUpdate[s] = 1; } if (R.SW.P[s].improved < 0) // If stagnation { quant[s] = quantAdapt (quant[s], ISS.q, param.qScale2); quantUpdate[s] = 1; } } break; case 3: // Global, according to all particle statuses prog = 0; stag = 0; for (s = 0; s < R.SW.S; s++) // Loop on the swarm { if (R.SW.P[s].improved > 0) prog++; // If progression if (R.SW.P[s].improved <= 0) stag++; // If stagnation } // if ( prog > 0 and initLinks==0) if (prog > R.SW.S / 2 && initLinks == 0) { ISS.q = quantAdapt (ISS.q, SS.q, param.qScale1); for (s = 0; s < R.SW.S; s++) quantUpdate[s] = 1; } // if ( stag > 0 and initLinks==1 ) // if ( stag > R.SW.S/2) // and initLinks==1 ) // if(stag==R.SW.S) if (initLinks == 1) // If no improvement { ISS.q = quantAdapt (ISS.q, SS.q, param.qScale2); for (s = 0; s < R.SW.S; s++) quantUpdate[s] = 1; // printf("\n %f",ISS.q.q[0]); } break; } // Velocity reinitialisations if (param.velReInit > 0) { for (s = 0; s < R.SW.S; s++) // Loop on the swarm { if (quantUpdate[s] == 1) { switch (param.quantis) { case 1: // Global case 3: // Global for (d = 0; d < SS.D; d++) { /* temp = ( alea( ISS.min[d], ISS.max[d] ) - alea( ISS.min[d], ISS.max[d] ) ) / 2; temp= quant[s].q[d]*(alea( 0, 4)-2 ); temp = (alea( 0, ISS.q.q[d] )+ alea( 0, ISS.q.q[d]))*(alea(0,2)-1); */ temp = alea_normal (0, param.w * ISS.q.q[d]); temp = temp * (ISS.max[d] - ISS.min[d]); if (fabs (temp) > fabs (R.SW.V[s].v[d])) R.SW.V[s].v[d] = temp; } break; case 2: // Individual for (d = 0; d < SS.D; d++) { /* temp = (alea( 0, quant[s].q[d] )+ alea(0, quant[s].q[d]))*(alea(0,2)-1); temp = alea_normal( 0, param.w *quant[s].q[d] ); */ temp = (alea(0, quant[s].q[d]) - alea (0, quant[s].q[d])); temp = temp * (ISS.max[d] - ISS.min[d]); if (fabs (temp) > fabs (R.SW.V[s].v[d])) R.SW.V[s].v[d] = temp; } break; // printf( "\n reinit velocity %i %f", // s,R.SW.V[s].v[0] ); } } } } errorPrev = error; switch (param.stop) { case 0: if (error > pb.epsilon && R.nEval < pb.evalMax) noStop = 0; // Won't stop else noStop = 1; // Will stop break; case 1: if (R.nEval < pb.evalMax) noStop = 0; // Won't stop else noStop = 1; // Will stop break; } } // End of loop on iterations // printf( "\n and the winner is ... %i", R.SW.best ); if (param.traj >= 0) // Save trajectory of a particle { traj (R, param.traj); } // fprintf( f_stag, "\nEND" ); return R; } // =========================================================== double alea (double a, double b) { // random number (uniform distribution) in [a b] double r; r = (double) rand (); r = r / RAND_MAX; return a + r * (b - a); } // =========================================================== double alea_cauchy (double mean, double gamma, double thetaMin, double thetaMax) { return (mean + gamma * tan (alea (thetaMin, thetaMax))); } // =========================================================== int alea_integer (int a, int b) { // Integer random number in [a b] int ir; double r; r = alea (0, 1); ir = (int) (a + r * (b + 1 - a)); if (ir > b) ir = b; return ir; } // =========================================================== double alea_levy (double mean, double scale, double alpha) { /* alpha in ]0,2] */ double r, s, t, u, v; u = pi * (alea (0, 1) - 0.5); do { v = -log (alea (0, 1)); } while (v == 0); t = sin (alpha * u) / pow (cos (u), 1 / alpha); s = pow (cos ((1 - alpha) * u) / v, (1 - alpha) / alpha); r = scale * t * s + mean; return r; } // =========================================================== double alea_normal (double mean, double std_dev) { /* Use the polar form of the Box-Muller transformation to obtain a pseudo random number from a Gaussian distribution */ double x1, x2, w, y1; // double y2; do { x1 = 2.0 * alea (0, 1) - 1.0; x2 = 2.0 * alea (0, 1) - 1.0; w = x1 * x1 + x2 * x2; } while (w >= 1.0); w = sqrt (-2.0 * log (w) / w); y1 = x1 * w; // y2 = x2 * w; y1 = y1 * std_dev + mean; return y1; } //====================================================== struct velocity alea_sphere(int D,double radius) { int j; double l; double pw; double r; struct velocity v; v.size=D; pw=1/(double)D; // ----------------------------------- Step 1. Random direction // This is a theorem ... l=0; for (j=0;j %f",d,x.x[d],y); } return xc; } // =========================================================== struct position distort (struct position x, struct SS SS, int way, double alpha) { double alphaT; int d; struct position y; if (alpha <= 0) return x; y.f = x.f; y.improved = x.improved; y.size = x.size; if (way == 1) alphaT = alpha; else alphaT = 1.0 / alpha; for (d = 0; d < x.size; d++) y.x[d] = distortX (x.x[d], SS.min[d], SS.max[d], alphaT); return y; } double distortX (double x, double min, double max, double alpha) { /* Distort a position */ double Dx = (max - min) / 2; double y; y = sign (x) * Dx * (1 - pow (1 - fabs (x) / Dx, alpha)); return y; } // =========================================================== double energy (struct swarm SW) { // Kinetic energy of the whole swarm double e = 0; int s; for (s = 0; s < SW.S; s++) e = e + energyS (SW.X[s], SW.XP[s]); return e; } double energyS (struct position x, struct position xp) { // Kinetic energy of a particle (mass=1/2) int d; double kinetic = 0; for (d = 0; d < x.size; d++) // Kinetic { kinetic = kinetic + pow (x.x[d] - xp.x[d], 2); } return kinetic; } // =========================================================== double gauss (double x, double mu, double var) { return exp (-0.5 * (x - mu) * (x - mu) / var) / sqrt (2 * pi * var); } // =========================================================== double normL (struct velocity v,double L) { // L-norm of a vector int d; double n; n = 0; for (d = 0; d < v.size; d++) n = n + pow(fabs(v.v[d]),L); n = pow (n, 1/L); return n; } // =========================================================== double distanceL (struct position x1, struct position x2,double L) { // Euclidean distance between two positions int d; double n; n = 0; for (d = 0; d < x1.size; d++) n = n + pow (fabs(x1.x[d] - x2.x[d]), L); n = pow (n, 1/L); return n; } // =========================================================== double max(double a,double b) { if (a>b) return a; return b; } double min(double a, double b) { if (a q.q[d]) quantA.q[d] = q.q[d]; } return quantA; /* //Random if ( scale < 1 ) for ( d = 0; d < quant.size; d++ ) { quantA.q[d] = quant.q[d] * alea( 0, scale ); } else for ( d = 0; d < quant.size; d++ ) { quantA.q[d] = quant.q[d] * alea( 1, scale ); if ( quantA.q[d] > q.q[d] ) quantA.q[d] = q.q[d]; } //printf("\n scale %f, q[0] %f, quant %f quantA[0] %.30f",scale, q[0],quant.q[0],quantA.q[0]); return quantA; */ } // =========================================================== struct position quantis (struct position x, struct SS SS, struct quantum quants, int option) { /* Quantisatition of a position Only values like x+k*q (k integer) are admissible */ int d; double qd; struct position quantx; quantx = x; for (d = 0; d < x.size; d++) { if (option == 2) qd = quants.q[d]; // Individual, may be different for each particle else qd = SS.q.q[d]; // Global if (qd > zero) // Note that qd can't be < 0 { qd = qd * (SS.max[d] - SS.min[d]) / 2; quantx.x[d] = qd * floor (0.5 + x.x[d] / qd); } } return quantx; } // =========================================================== struct position queen (struct swarm SW, struct problem pb, struct SS SS, double distort) { int d; struct position q; struct quantum dummy; // ={0}; int s; q.size = SW.X[0].size; for (d = 0; d < q.size; d++) { q.x[d] = SW.P[0].x[d]; for (s = 1; s < SW.S; s++) { q.x[d] = q.x[d] + SW.P[s].x[d]; } q.x[d] = q.x[d] / SW.S; } // Evaluation // It takes quantisation into account q = quantis (q, SS, dummy, 0); // If list, clamp to the nearest admissible value q = clampToList (q, SS); // printf("\nqueen eval"); q.f = perf (q, pb.function, pb.offset, pb.objective, SS, distort); return q; } // =========================================================== void traj (struct result R, int rank) { // Save trajectory of particle "rank" int d; fprintf (f_traj, "%i ", (int) R.nEval); for (d = 0; d < R.SW.X[0].size; d++) { fprintf (f_traj, "%f ", R.SW.X[rank].x[d]); // Current // position // fprintf( f_traj, "%f ", R.SW.P[rank].x[d] );// Best previous } fprintf (f_traj, "%i ", R.SW.best); fprintf (f_traj, "\n"); } // =========================================================== double perf (struct position x, int function, double offset[], double objective, struct SS SS, double alpha) { // Evaluate the fitness value for the particle of rank s int d; int k; double f, p, xd, x1, x2; int out; double s11, s12, s21, s22; double t0, tt, t1; struct position xs; if (alpha < 1 || alpha > 1) xs = distort (x, SS, -1, alpha); else xs = x; switch (function) { case 0: // Parabola (Sphere) f = 0; for (d = 0; d < xs.size; d++) { xd = xs.x[d] - offset[d]; f = f + xd * xd; } break; case 1: // Griewank f = 0; p = 1; for (d = 0; d < xs.size; d++) { xd = xs.x[d] - offset[d]; f = f + xd * xd; p = p * cos (xd / sqrt ((double) (d + 1))); } f = f / 4000 - p + 1; break; case 2: // Rosenbrock f = 0; t0 = xs.x[0] - offset[0] + 1; // Solution on (0,...0) when // offset=0 for (d = 1; d < xs.size; d++) { t1 = xs.x[d] - offset[d] + 1; tt = 1 - t0; f += tt * tt; tt = t1 - t0 * t0; f += 100 * tt * tt; t0 = t1; } break; case 3: // Rastrigin k = 10; f = 0; for (d = 0; d < xs.size; d++) { xd = xs.x[d] - offset[d]; f =f+ xd * xd - k * cos (2 * pi * xd); } f =f+ xs.size * k; break; case 4: // 2D Tripod function // Note that there is a big discontinuity right on the solution // point. x1 = xs.x[0] - offset[0]; x2 = xs.x[1] - offset[1]; s11 = (1.0 - sign (x1)) / 2; s12 = (1.0 + sign (x1)) / 2; s21 = (1.0 - sign (x2)) / 2; s22 = (1.0 + sign (x2)) / 2; //f = s21 * (fabs (x1) - x2); // Solution on (0,0) f = s21 * (fabs (x1) +fabs(x2+50)); // Solution on (0,-50) f = f + s22 * (s11 * (1 + fabs (x1 + 50) + fabs (x2 - 50)) + s12 * (2 + fabs (x1 - 50) + fabs (x2 - 50))); break; case 99: // Test // --------------- Bill Langdon out = 0; for (d = 0; d < SS.D; d++) { if (xs.x[d] < SS.min[d] || xs.x[d] > SS.max[d]) out = 1; } if (out == 1) f = 0; else f = 2 * (xs.x[0] - xs.x[1]); // on [-10,10] f = 40 - f; // we are looking for a minimum break; // -----------------Bill Langdon out = 0; for (d = 0; d < SS.D; d++) { if (xs.x[d] < SS.min[d] || xs.x[d] > SS.max[d]) out = 1; } if (out == 1) f = 0; else f = 0.063 * xs.x[0]; // on [-10,10] f = 0.63 - f; break; // ---------------------------Bill Langdon out = 0; for (d = 0; d < SS.D; d++) { if (xs.x[d] < SS.min[d] || xs.x[d] > SS.max[d]) out = 1; } if (out == 1) f = 0; else f = 0.00124 * xs.x[0] * xs.x[0] * xs.x[1]; // on [-10,10] f = 1.24 - f; break; // --------------------------- A 1D example x1 = xs.x[0]; f = x1 * fabs (sin (1 / x1)) / (x1 * x1 + 1) + 0.4667; break; } return fabs (f - objective); }