refac: header consolidation + memory fix in de

diff --git a/include/cgo.h b/include/cgo.h
new file mode 100644 (file)
index 0000000..f81d509
--- /dev/null
+++ b/include/cgo.h
@@ -0,0 +1,3 @@
+#include "de.h"
+#include "pso.h"
+#include "vector.h"

diff --git a/include/pso.h b/include/pso.h
new file mode 100644 (file)
index 0000000..b888dbe
--- /dev/null
+++ b/include/pso.h
@@ -0,0 +1,36 @@
+#include "rand.h"
+
+typedef struct particle particle;
+
+typedef struct {
+  double *x;
+  double fitness;
+} particle_state;
+
+typedef struct {
+  int dim;
+  int particlec;
+  range *rrange;
+  double (*f)(double*);  // fitness function
+  double eps;            // convergence tolerance
+  double w;              // inertia
+  double c;              // cognitive factor
+  double s;              // social factor
+} pso_parameters;
+
+typedef struct {
+  particle *particles;
+  particle *global_best;
+} swarm;
+
+struct particle {
+  particle_state current;
+  particle_state best;
+  double *velocity;
+};
+
+void swarm_alloc(pso_parameters *par, swarm *sw);
+void swarm_free(pso_parameters *par, swarm *sw);
+void swarm_populate(pso_parameters *par, swarm *sw, const range *rrange);
+void swarm_velocity_update(pso_parameters *par, swarm *sw);
+void pso_optimize(pso_parameters *par, particle_state *res);

diff --git a/include/rand.h b/include/rand.h
index 7c0373dfd300c5c36dedc4bca81ef2d6b69bb302..172d6947654e4ccc9a04d260894f2d79409f1a34 100644 (file)
--- a/include/rand.h
+++ b/include/rand.h
@@ -1,3 +1,6 @@
+#ifndef RAND_H
+#define RAND_H
+
 typedef struct {
   double min;
   double max;
@@ -6,3 +9,5 @@ typedef struct {
 double rands();
 double randrange(double min, double max);
 void vec_pick(const int dim, const range *rrange, double *v);
+
+#endif

diff --git a/include/swarm.h b/include/swarm.h
deleted file mode 100644 (file)
index cfeb33f..0000000
--- a/include/swarm.h
+++ /dev/null
@@ -1,37 +0,0 @@
-#include "rand.h"
-
-typedef struct particle particle;
-
-typedef struct {
-  double *x;
-  double fitness;
-} particle_state;
-
-typedef struct {
-  int dim;
-  int particlec;
-  range *rrange;
-  double (*f)(double*);  // fitness function
-  double eps;            // convergence tolerance
-  double w;              // inertia
-  double c;              // cognitive factor
-  double s;              // social factor
-} swarm_parameters;
-
-typedef struct {
-  particle *particles;
-  particle *global_best;
-} swarm;
-
-struct particle {
-  particle_state current;
-  particle_state best;
-  double *velocity;
-};
-
-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/de.c b/src/de.c
index 441312b34a1846c7b70ad0502bf66b9b58ee8126..0e9f72c0c7900b173ecc2167525e540b74b65e3b 100644 (file)
--- a/src/de.c
+++ b/src/de.c
@@ -4,26 +4,25 @@
 
 void de_optimize(de_parameters *par, double *res) {
     int l = par->particlec * par->dim;
-    double *mem = malloc(sizeof(double) * l * 2);
+    double *mem = malloc(sizeof(double) * (l + par->dim));
     int conv = 0;
     double *xs = mem;
-    double *ys = mem + l;
-    for (double *x = xs; x < ys; x += par->dim)
+    double *y = mem + l;
+    for (double *x = xs; x < y; x += par->dim)
         vec_pick(par->dim, par->rrange, x);
     while (conv < 1000) {
-        for (double *x = xs; x < ys; x += par->dim) {
-            double *y = x + l;
+        for (double *x = xs; x < y; x += par->dim) {
             double *picks[4] = { x };  // agents { x, a, b, c }
-            for (int i = 1; i < 4; i++) {
-                double *pick = xs + (rand() % (par->particlec - i));
-                picks[i] = pick;
-                for (int k = 0; k < i; k++)
-                    if (pick >= picks[k]) picks[i]++;
+            for (int a = 1; a < 4; a++) {
+                double *pick = xs + (rand() % (par->particlec - a));
+                picks[a] = pick;
+                for (int ap = 0; ap < a; ap++)
+                    if (pick >= picks[ap]) picks[a]++;
             }
             int r = rand() % par->dim;
-            for (int i = 0; i < par->dim; i++) {
-                if (i == r || rands() < par->cr) y[i] = picks[1][i] + par->dw * (picks[2][i] - picks[3][i]);
-                else y[i] = x[i];
+            for (int d = 0; d < par->dim; d++) {
+                if (d == r || rands() < par->cr) y[d] = picks[1][d] + par->dw * (picks[2][d] - picks[3][d]);
+                else y[d] = x[d];
             }
             if (par->f(y) < par->f(x)) {
                 vec_copy(par->dim, y, x);
@@ -34,7 +33,7 @@ void de_optimize(de_parameters *par, double *res) {
     }
     double *best = xs;
     double bestf = par->f(best);
-    for (double *x = xs; x < ys; x += par->dim) {
+    for (double *x = xs; x < y; x += par->dim) {
         double cur = par->f(x);
         if (cur < bestf) {
             best = x;

diff --git a/src/main.c b/src/main.c
index b75439c040723e5ae8ae8cb128997d0095a34d12..f900b93f2385fc628b79bc5386fa48dd385e1f3d 100644 (file)
--- a/src/main.c
+++ b/src/main.c
@@ -1,5 +1,7 @@
 #include <argp.h>
 #include <string.h>
+#include <math.h>
+#include "cgo.h"
 
 enum algorithm {
     de,
@@ -39,11 +41,37 @@ static error_t parse_opt(int key, char *arg, struct argp_state *state) {
 
 static struct argp argp = { options, parse_opt, args_doc, doc };
 
+double b0(double *x) {
+  return -x[0];
+}
+
+double b1(double *x) {
+  return -x[1];
+}
+
+double b2(double *x) {
+  return x[0] - 10;
+}
+
+double b3(double *x) {
+  return x[1] - 15;
+}
+
+double f(double *x) {
+  return exp(b0(x)) + exp(b1(x)) + exp(b2(x)) + exp(b3(x)) + x[0] - x[1];
+}
+
 int main(int argc, char *argv[]) {
     struct arguments args;
     args.algorithm = de;
     args.outfile = "out.log";
     argp_parse(&argp, argc, argv, 0, 0, &args);
-    
+    range sample[2] = { (range) { -100.0, 100.0 }, { -100.0, 100.0 } };
+    de_parameters par = {2, 1000, sample, f};
+    double *res = vec_alloc(2);
+    de_optimize(&par, res);
+    printf("Converged at x=");
+    vec_print(2, res);
+    printf(" with fitness f(x)=%f\n", f(res));
     return 0;
 }

diff --git a/src/pso.c b/src/pso.c
new file mode 100644 (file)
index 0000000..e7efa80
--- /dev/null
+++ b/src/pso.c
@@ -0,0 +1,88 @@
+#include <math.h>
+#include <omp.h>
+#include <stdlib.h>
+#include "pso.h"
+#include "vector.h"
+
+void swarm_alloc(pso_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 swarm_free(pso_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(pso_parameters *par, swarm *sw, const range *rrange) {
+  sw->global_best = &sw->particles[0];
+  #pragma omp parallel for
+  for (int p = 0; p < par->particlec; p++) {
+    particle *pt = &sw->particles[p];
+    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 = par->f(pt->current.x);
+    pt->best.fitness = pt->current.fitness;
+    #pragma omp critical
+    {
+      if (pt->current.fitness < sw->global_best->current.fitness)
+        sw->global_best = pt;
+    }
+  }
+}
+
+void swarm_velocity_update(pso_parameters *par, swarm *sw) {
+  #pragma omp parallel for
+  for (int p = 0; p < par->particlec; p++) {
+    particle *pt = &sw->particles[p];
+    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 < par->particlec; p++) {
+    particle *pt = &sw->particles[p];
+    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)
+        sw->global_best = pt;
+    }
+  }
+}
+
+void pso_optimize(pso_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(par, &sw);
+    if (fabs(sw.global_best->current.fitness - best) < par->eps)
+      conv++;
+    else {
+      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);
+}

diff --git a/src/swarm.c b/src/swarm.c
deleted file mode 100644 (file)
index 028f51b..0000000
--- a/src/swarm.c
+++ /dev/null
@@ -1,88 +0,0 @@
-#include <math.h>
-#include <omp.h>
-#include <stdlib.h>
-#include "swarm.h"
-#include "vector.h"
-
-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 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_parameters *par, swarm *sw, const range *rrange) {
-  sw->global_best = &sw->particles[0];
-  #pragma omp parallel for
-  for (int p = 0; p < par->particlec; p++) {
-    particle *pt = &sw->particles[p];
-    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 = par->f(pt->current.x);
-    pt->best.fitness = pt->current.fitness;
-    #pragma omp critical
-    {
-      if (pt->current.fitness < sw->global_best->current.fitness)
-        sw->global_best = pt;
-    }
-  }
-}
-
-void swarm_velocity_update(swarm_parameters *par, swarm *sw) {
-  #pragma omp parallel for
-  for (int p = 0; p < par->particlec; p++) {
-    particle *pt = &sw->particles[p];
-    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 < par->particlec; p++) {
-    particle *pt = &sw->particles[p];
-    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)
-        sw->global_best = pt;
-    }
-  }
-}
-
-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(par, &sw);
-    if (fabs(sw.global_best->current.fitness - best) < par->eps)
-      conv++;
-    else {
-      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);
-}