refactor: split off randomization

diff --git a/include/rand.h b/include/rand.h
new file mode 100644 (file)
index 0000000..7c0373d
--- /dev/null
+++ b/include/rand.h
@@ -0,0 +1,8 @@
+typedef struct {
+  double min;
+  double max;
+} range;
+
+double rands();
+double randrange(double min, double max);
+void vec_pick(const int dim, const range *rrange, double *v);

diff --git a/include/swarm.h b/include/swarm.h
index b3f58e95fbffdb865b51e9b79f263323a891e0d6..cfeb33f91675443fb2414b6405d0684a4f1d05cd 100644 (file)
--- a/include/swarm.h
+++ b/include/swarm.h
@@ -1,9 +1,6 @@
-typedef struct particle particle;
+#include "rand.h"
 
-typedef struct {
-  double min;
-  double max;
-} range;
+typedef struct particle particle;
 
 typedef struct {
   double *x;
@@ -13,6 +10,7 @@ typedef struct {
 typedef struct {
   int dim;
   int particlec;
+  range *rrange;
   double (*f)(double*);  // fitness function
   double eps;            // convergence tolerance
   double w;              // inertia
@@ -21,7 +19,6 @@ typedef struct {
 } swarm_parameters;
 
 typedef struct {
-  swarm_parameters params;
   particle *particles;
   particle *global_best;
 } swarm;
@@ -32,7 +29,9 @@ struct particle {
   double *velocity;
 };
 
-void pick(const int dim, const range *rrange, double *v);
-void swarm_populate(swarm *sw, const range *rrange);
-void swarm_velocity_update(swarm *sw);
-void swarm_optimize(swarm *sw);
+void swarm_alloc(swarm_parameters *par, swarm *sw);
+void swarm_free(swarm_parameters *par, swarm *sw);
+void swarm_populate(swarm_parameters *par, swarm *sw, const range *rrange);
+void swarm_velocity_update(swarm_parameters *par, swarm *sw);
+
+void swarm_optimize(swarm_parameters *par, particle_state *res);

diff --git a/src/main.c b/src/main.c
index 94ddcfad2e6860320365354963e90490b638d0f5..430874c586d4551f7cf0d9fa24ffff8038cf57f6 100644 (file)
--- a/src/main.c
+++ b/src/main.c
@@ -24,24 +24,25 @@ double f(double *x) {
 }
 
 int main(int argc, char *argv[]) {
-  swarm sw = {
+  range rrange[2] = {
+    (range) { -100.0, 100.0 },
+    (range) { -100.0, 100.0 }
+  };
+  swarm_parameters par = {
     2,
     100,
+    rrange,
     f,
     pow(2, 40),
     0.8,
     2.05,
     2.05
   };
-  range rrange[2] = {
-    (range) { -100.0, 100.0 },
-    (range) { -100.0, 100.0 }
-  };
-  swarm_populate(&sw, rrange);
-  swarm_optimize(&sw);
+  particle_state best = { vec_alloc(par.dim), 0.0 };
+  swarm_optimize(&par, &best);
   printf("apso: converged in x=");
-  vec_print(sw.params.dim, sw.global_best->current.x);
-  printf(" with fitness f(x)=%f\n", sw.global_best->current.fitness);
+  vec_print(par.dim, best.x);
+  printf(" with fitness f(x)=%f\n", best.fitness);
   // GtkApplication *app = gui_init();
   // return g_application_run(G_APPLICATION(app), argc, argv);
   return 0;

diff --git a/src/rand.c b/src/rand.c
new file mode 100644 (file)
index 0000000..323834e
--- /dev/null
+++ b/src/rand.c
@@ -0,0 +1,14 @@
+#include <stdlib.h>
+#include "rand.h"
+
+double rands() {
+  return (double) rand() / RAND_MAX;
+}
+
+double randrange(double min, double max) {
+  return rands() * (max - min) + min;
+}
+
+void vec_pick(const int dim, const range *rrange, double *v) {
+  for (int d = 0; d < dim; d++) v[d] = randrange(rrange[d].min, rrange[d].max);
+}

diff --git a/src/swarm.c b/src/swarm.c
index c4074df91e7bca599217d35397e33c936162a523..028f51bc92dc08d7e9c92e036730bfd58cdb02bc 100644 (file)
--- a/src/swarm.c
+++ b/src/swarm.c
@@ -1,35 +1,39 @@
 #include <math.h>
-#include <stdlib.h>
 #include <omp.h>
+#include <stdlib.h>
 #include "swarm.h"
 #include "vector.h"
 
-double rands() {
-  return (double) rand() / RAND_MAX;
-}
-
-double randrange(double min, double max) {
-  return rands() * (max - min) + min;
+void swarm_alloc(swarm_parameters *par, swarm *sw) {
+  sw->particles = malloc(sizeof(particle) * par->particlec);
+  #pragma omp parallel for
+  for (int p = 0; p < par->particlec; p++) {
+    sw->particles[p].current.x = vec_alloc(par->dim);
+    sw->particles[p].best.x = vec_alloc(par->dim);
+    sw->particles[p].velocity = vec_alloc(par->dim);
+  }
 }
 
-void pick(const int dim, const range *rrange, double *v) {
-  for (int d = 0; d < dim; d++) v[d] = randrange(rrange[d].min, rrange[d].max);
+void swarm_free(swarm_parameters *par, swarm *sw) {
+  #pragma omp parallel for
+  for (int p = 0; p < par->particlec; p++) {
+    free(sw->particles[p].current.x);
+    free(sw->particles[p].best.x);
+    free(sw->particles[p].velocity);
+  }
+  free(sw->particles);
 }
 
-void swarm_populate(swarm *sw, const range *rrange) {
-  sw->particles = malloc(sizeof(particle) * sw->params.particlec);
+void swarm_populate(swarm_parameters *par, swarm *sw, const range *rrange) {
   sw->global_best = &sw->particles[0];
   #pragma omp parallel for
-  for (int p = 0; p < sw->params.particlec; p++) {
+  for (int p = 0; p < par->particlec; p++) {
     particle *pt = &sw->particles[p];
-    pt->current.x = vec_alloc(sw->params.dim);
-    pt->best.x = vec_alloc(sw->params.dim);
-    pt->velocity = vec_alloc(sw->params.dim);
-    for (int d = 0; d < sw->params.dim; d++) {
-      pick(sw->params.dim, rrange, pt->current.x);
-      vec_copy(sw->params.dim, pt->current.x, pt->best.x);
+    for (int d = 0; d < par->dim; d++) {
+      vec_pick(par->dim, rrange, pt->current.x);
+      vec_copy(par->dim, pt->current.x, pt->best.x);
     }
-    pt->current.fitness = sw->params.f(pt->current.x);
+    pt->current.fitness = par->f(pt->current.x);
     pt->best.fitness = pt->current.fitness;
     #pragma omp critical
     {
@@ -39,19 +43,22 @@ void swarm_populate(swarm *sw, const range *rrange) {
   }
 }
 
-void swarm_velocity_update(swarm *sw) {
+void swarm_velocity_update(swarm_parameters *par, swarm *sw) {
   #pragma omp parallel for
-  for (int p = 0; p < sw->params.particlec; p++) {
+  for (int p = 0; p < par->particlec; p++) {
     particle *pt = &sw->particles[p];
-    for (int d = 0; d < sw->params.dim; d++)
-      pt->velocity[d] = sw->params.w * pt->velocity[d] + sw->params.c * rands() * (pt->best.x[d] - pt->current.x[d]) + sw->params.s * rands() * (sw->global_best->current.x[d] - pt->current.x[d]);
+    for (int d = 0; d < par->dim; d++)
+      pt->velocity[d] = par->w * pt->velocity[d] + par->c * rands() * (pt->best.x[d] - pt->current.x[d]) + par->s * rands() * (sw->global_best->current.x[d] - pt->current.x[d]);
   }
   #pragma omp parallel for
-  for (int p = 0; p < sw->params.particlec; p++) {
+  for (int p = 0; p < par->particlec; p++) {
     particle *pt = &sw->particles[p];
-    vec_add(sw->params.dim, pt->velocity, pt->current.x);
-    pt->current.fitness = sw->params.f(pt->current.x);
-    if (pt->current.fitness < pt->best.fitness) pt->best = pt->current;
+    vec_add(par->dim, pt->velocity, pt->current.x);
+    pt->current.fitness = par->f(pt->current.x);
+    if (pt->current.fitness < pt->best.fitness) {
+      vec_copy(par->dim, pt->current.x, pt->best.x);
+      pt->best.fitness = pt->current.fitness;
+    }
     #pragma omp critical
     {
       if (pt->current.fitness < sw->global_best->current.fitness)
@@ -60,16 +67,22 @@ void swarm_velocity_update(swarm *sw) {
   }
 }
 
-void swarm_optimize(swarm *sw) {
+void swarm_optimize(swarm_parameters *par, particle_state *res) {
+  swarm sw;
+  swarm_alloc(par, &sw);
+  swarm_populate(par, &sw, par->rrange);
   int conv = 0;
   double best = INFINITY;
   while (conv < 1000) {
-    swarm_velocity_update(sw);
-    if (fabs(sw->global_best->current.fitness - best) < sw->params.eps)
+    swarm_velocity_update(par, &sw);
+    if (fabs(sw.global_best->current.fitness - best) < par->eps)
       conv++;
     else {
-      best = sw->global_best->current.fitness;
+      best = sw.global_best->current.fitness;
       conv = 0;
     }
   }
+  vec_copy(par->dim, sw.global_best->best.x, res->x);
+  res->fitness = sw.global_best->best.fitness;
+  swarm_free(par, &sw);
 }