fix: add random jitter to simluation

include/swarm.hpp

@@ -4,6 +4,8 @@
 #include <cstddef>
 
 #include <functional>
+#include <limits>
+#include <random>
 #include <utility>
 #include <vec.hpp>
 #include <vector>
@@ -94,6 +96,14 @@ struct Particle : public Agent<vec<N>> {
 
     void move(float dt=1) override { position = position + velocity * dt; }
 
+    void observe() {
+    float y = this->f(position);
+    if(y < pb) {
+        pb = y;
+        pb_pos = position;
+    }
+    }
+
     void step() override {
 	// Get the new parameter values
 	vec<3> params = alg(kviscosity, knostalgia, kpeer_pressure,
@@ -106,14 +116,8 @@ struct Particle : public Agent<vec<N>> {
 	curr_iter++;
 
     velocity = viscosity * velocity
-	    + nostalgia * (pb_pos - position)
-	    + peer_pressure * (peer.best().second - position);
-
-	float y = f(position);
-	if(y < pb) {
-	    pb = y;
-	    pb_pos = position;
-	}
+        + nostalgia * (random_coefficients() * (pb_pos - position))
+        + peer_pressure * (random_coefficients() * (peer.best().second - position));
     };
 
     const float kviscosity = 0.9f;
@@ -135,6 +139,17 @@ struct Particle : public Agent<vec<N>> {
     const vec<N> &get_velocity() const { return velocity; };
 
 private:
+    static vec<N> random_coefficients() {
+    static thread_local std::mt19937 generator(std::random_device{}());
+    static thread_local std::uniform_real_distribution<float> distribution(0.0f, 1.0f);
+
+    vec<N> coefficients;
+    for(std::size_t i = 0; i < N; ++i)
+        coefficients.c[i] = distribution(generator);
+
+    return coefficients;
+    }
+
     A alg;
     vec<N> position;
     vec<N> velocity;
@@ -160,8 +175,10 @@ struct Swarm : public Agent<vec<N>> {
     void step() override {
 	if(particles.empty()) return;
 
-	// find best before step() on each particle
-	// otherwise, if this were in best(), each p.step() would update it
+    for(auto &p : particles) {
+        p.observe();
+    }
+
 	const Particle<N> *b = &particles[0];
 	for(const auto &p : particles) {
 	    if(b->best().first > p.best().first) b = &p;
@@ -185,6 +202,12 @@ struct Swarm : public Agent<vec<N>> {
 	particles.push_back(
 	    Particle<N>(this->f, *this, pos, vel, max_iter)
 	);
+
+    const auto [particle_best_val, particle_best_pos] = particles.back().best();
+    if(particle_best_val < best_val) {
+        best_val = particle_best_val;
+        best_pos = particle_best_pos;
+    }
     }
 
     const std::vector<Particle<N>>& get_particles() { return particles; };
@@ -194,6 +217,6 @@ struct Swarm : public Agent<vec<N>> {
 private:
     unsigned max_iter;
     vec<N> best_pos;
-    float best_val;
+    float best_val = std::numeric_limits<float>::infinity();
     std::vector<Particle<N>> particles;
 };