initial commit

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
index 0000000..a806caf
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,2 @@
+opt: src/particle_swarm.c src/vector.c
+       gcc -Iinclude -lm -fopenmp src/vector.c src/particle_swarm.c -o opt

diff --git a/include/particle_swarm.h b/include/particle_swarm.h
new file mode 100644 (file)
index 0000000..a9fdedb
--- /dev/null
+++ b/include/particle_swarm.h
@@ -0,0 +1,27 @@
+typedef struct particle particle;
+
+typedef struct {
+  double min;
+  double max;
+} range;
+
+typedef struct {
+  double *x;
+  double fitness;
+} particle_state;
+
+typedef struct {
+  int particlec;
+  particle *particles;
+  particle *global_best;
+} swarm;
+
+struct particle {
+  particle_state current;
+  particle_state best;
+  double *velocity;
+};
+
+void pick(const int dim, const range *vrange, double *v);
+void swarm_init(swarm *sw, const int dim, const range *vrange);
+void swarm_velocity_update(swarm *sw, const int dim);

diff --git a/include/vector.h b/include/vector.h
new file mode 100644 (file)
index 0000000..7e7160a
--- /dev/null
+++ b/include/vector.h
@@ -0,0 +1,8 @@
+double *vec_alloc(const int dim);
+void vec_print(const int dim, const double *v);
+void vec_copy(const int dim, const double *v, double *vp);
+void vec_scale(const int dim, const double s, double *v);
+void vec_inv(const int dim, double *v);
+void vec_add(const int dim, const double *b, double *a);  // a = a + b
+void vec_sub(const int dim, const double *b, double *a);  // a = a - b
+double vec_dotprod(const int dim, const double *a, const double *b);

diff --git a/src/particle_swarm.c b/src/particle_swarm.c
new file mode 100644 (file)
index 0000000..f1e2b41
--- /dev/null
+++ b/src/particle_swarm.c
@@ -0,0 +1,102 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "particle_swarm.h"
+#include "vector.h"
+
+double c = 2.05;
+double s = 2.05;
+double w = 0.8;
+
+double rands() {
+  return (double) rand() / RAND_MAX;
+}
+
+double randrange(double min, double max) {
+  return rands() * (max - min) + min;
+}
+
+void pick(const int dim, const range *vrange, double *v) {
+  for (int d = 0; d < dim; d++) v[d] = randrange(vrange[d].min, vrange[d].max);
+}
+
+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];
+}
+
+void swarm_init(swarm *sw, const int dim, const range *vrange) {
+  sw->particles = malloc(sizeof(particle) * sw->particlec);
+  sw->global_best = &sw->particles[0];
+  for (int p = 0; p < sw->particlec; p++) {
+    particle *pt = &sw->particles[p];
+    pt->current.x = vec_alloc(dim);
+    pt->best.x = vec_alloc(dim);
+    pt->velocity = vec_alloc(dim);
+    for (int d = 0; d < dim; d++) {
+      pick(dim, vrange, pt->current.x);
+      vec_copy(dim, pt->current.x, pt->best.x);
+    }
+    pt->current.fitness = f(pt->current.x);
+    pt->best.fitness = pt->current.fitness;
+    if (pt->current.fitness < sw->global_best->current.fitness) sw->global_best = pt;
+  }
+}
+
+void swarm_velocity_update(swarm *sw, const int dim) {
+  for (int p = 0; p < sw->particlec; p++) {
+    particle *pt = &sw->particles[p];
+    for (int d = 0; d < dim; d++)
+      pt->velocity[d] = w * pt->velocity[d] + c * rands() * (pt->best.x[d] - pt->current.x[d]) + s * rands() * (sw->global_best->current.x[d] - pt->current.x[d]);
+  }
+  for (int p = 0; p < sw->particlec; p++) {
+    particle *pt = &sw->particles[p];
+    vec_add(dim, pt->velocity, pt->current.x);
+    pt->current.fitness = f(pt->current.x);
+    if (pt->current.fitness < pt->best.fitness) pt->best = pt->current;
+    if (pt->current.fitness < sw->global_best->current.fitness) sw->global_best = &sw->particles[p];
+  }
+}
+
+int main() {
+  double eps = pow(2, 40);
+  int dim = 2;
+  int particlec = 100;
+  range vrange[2] = {
+    (range) { -100.0, 100.0 },
+    (range) { -100.0, 100.0 }
+  };
+  swarm sw;
+  sw.particlec = particlec;
+  swarm_init(&sw, dim, vrange);
+  int conv = 0;
+  double best = 0.0;
+  while (conv < 10000) {
+    swarm_velocity_update(&sw, dim);
+    if (fabs(sw.global_best->current.fitness - best) < eps)
+      conv++;
+    else {
+      best = sw.global_best->current.fitness;
+      conv = 0;
+    }
+  }
+  printf("PSO: converged in x=");
+  vec_print(dim, sw.global_best->current.x);
+  printf(" with fitness f(x)=%f\n", sw.global_best->current.fitness);
+  return 0;
+}

diff --git a/src/vector.c b/src/vector.c
new file mode 100644 (file)
index 0000000..fa39447
--- /dev/null
+++ b/src/vector.c
@@ -0,0 +1,46 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+double *vec_alloc(const int dim)
+{
+  return malloc(sizeof(double) * dim);
+}
+
+void vec_print(const int dim, const double *v)
+{
+  printf("[");
+  for (int d = 0; d < dim - 1; d++) printf("%f, ", v[d]);
+  printf("%f]", v[dim - 1]);
+}
+
+void vec_copy(const int dim, const double *v, double *vp)
+{
+  for (int d = 0; d < dim; d++) vp[d] = v[d];
+}
+
+void vec_scale(const int dim, const double s, double *v)
+{
+  for (int d = 0; d < dim; d++) v[d] *= s;
+}
+
+void vec_inv(const int dim, double *v)
+{
+  vec_scale(dim, -1.0, v);
+}
+
+void vec_add(const int dim, const double *b, double *a)
+{
+  for (int d = 0; d < dim; d++) a[d] += b[d];
+}
+
+void vec_sub(const int dim, const double *b, double *a)
+{
+  for (int d = 0; d < dim; d++) a[d] -= b[d];
+}
+
+double vec_dotprod(const int dim, const double *a, const double *b)
+{
+  double sum = 0.0;
+  for (int d = 0; d < dim; d++) sum += a[d] * b[d];
+  return sum;
+}