Add genetic_search

.gitignore

@@ -28,4 +28,5 @@ main.pdf
 /src/dump
 /src/random_local_search
 /src/sa_local_search_*
+/src/genetic_search_*
 /src/discr_calc/*_discr

src/halton.cpp

@@ -1,6 +1,7 @@
 #include "halton.h"
 #include <cassert>
 
+using namespace std;
 
 // compute \phi_p^\pi(i)
 double phi(int p, permutation &pi, int i) {
@@ -26,6 +27,13 @@ permutation& halton::get_pi(int i) {
   return pi[i];
 }
 
+void halton::set_pis(vector<permutation> &pis) {
+  assert((int)pis.size() == dim);
+  for(int i = 0; i < dim; i++)
+    assert(pis[i].get_size() == p[i]);
+  pi = pis;
+}
+
 void halton::compute(int i, double* res) {
   for(int d = 0; d < dim; d++) {
     res[d] = phi(p[d], pi[d], i);

src/halton.h

@@ -10,6 +10,7 @@ class halton {
   public:
     halton(int dim, int *p); // p is not deep-copied
     permutation& get_pi(int i);
+    void set_pis(std::vector<permutation> &pis);
     int get_p(int i);
     void compute(int i, double* res);
   protected:

src/main.cpp

@@ -21,7 +21,9 @@ int main(int argc, char **argv) {
 
   // random_search rs(dim, npoints, p);
   // rs.run(iterations);
-  sa_local_search s(dim, npoints, p, 0.992, 0.01/dim);
+  // sa_local_search s(dim, npoints, p, 0.992, 0.01/dim);
+  // s.run(iterations);
+  genetic_search s(dim, npoints, p, 10, 10, 0.5);
   s.run(iterations);
 
   return EXIT_SUCCESS;

src/permutation.cpp

@@ -43,6 +43,10 @@ bool permutation::check() {
   return true;
 }
 
+int permutation::get_size() {
+  return size;
+}
+
 // Algorithm:
 // for each i (randomly : we use a permutation :P (it always starts with 0 but
 // it's OK since a[0]=b[0]=0 and also we want ret[0]=0))

src/permutation.h

@@ -11,6 +11,7 @@ class permutation {
     void random();
     int operator[](int i); // return the image of i (i and its image 0-based)
     bool check(); // check invariant : sigma[0] = 0 and it is a permutation
+    int get_size();
     friend permutation permutation_crossover(permutation &a, permutation &b);
   protected:
     int size;

src/search.cpp

@@ -1,7 +1,10 @@
 #include "search.h"
 #include <cmath>
 #include <cassert>
+#include <algorithm> // sort
+#include <tuple> // tie, ignore
 
+using namespace std;
 
 /**********************    search    ***********************/
 
@@ -69,9 +72,12 @@ void local_search::init() {
   check();
 }
 
-
 void local_search::random_neighbour() {
   undoable = true;
+  _random_neighbour();
+}
+
+void local_search::_random_neighbour() {
   i = rand()%dim;
   int p = ha.get_p(i);
   t1 = 1+rand()%(p-1);
@@ -85,6 +91,10 @@ void local_search::random_neighbour() {
 void local_search::undo() {
   assert(undoable);
   undoable = false;
+  _undo();
+}
+
+void local_search::_undo() {
   ha.get_pi(i).transpose(t1, t2);
 }
 
@@ -115,7 +125,7 @@ bool random_local_search::accept(double previous, double current) {
 }
 
 
-/**********************    sa_search    ***********************/
+/**********************    sa_local_search    ***********************/
 
 sa_local_search::sa_local_search(int dim, int npoints, int *p, double lambda, double temp)
     :local_search(dim, npoints, p), lambda(lambda), temp(temp) {
@@ -130,3 +140,81 @@ bool sa_local_search::accept(double previous, double current) {
   //printf("%lf\t%lf\n", temp, current);
   return current < previous || (double)rand()/RAND_MAX < exp((previous - current)/temp);
 }
+
+
+/**********************    genetic_search    ***********************/
+
+genetic_search::genetic_search(int dim, int npoints, int *p, int mu, int lambda, double c) :
+    search(dim, npoints, p), mu(mu), lambda(lambda), c(c) {
+  assert(mu > 0);
+  assert(lambda > 0);
+  assert(0 <= c && c <= 1);
+
+  filename = std::string("genetic_search_") + std::to_string(mu)
+    + "_" + std::to_string(lambda)
+    + "_" + std::to_string(c);
+
+  //init(); // not necessary here
+
+  genes.reserve(mu+lambda);
+  vector<permutation> base;
+
+  // basic vector of permutations
+  base.reserve(dim);
+  for(int i = 0; i < dim; i++)
+    base.emplace_back(p[i]);
+
+  // generate mu random genes
+  for(int i = 0; i < mu; i++) {
+    genes.push_back(make_pair(1., base));
+    for(int j = 0; j < dim; j++)
+      genes[i].second[j].random();
+  }
+
+  // "reserve" lambda other genes (later we'll only change permutations already
+  // created here)
+  for(int i = 0; i < lambda; i++)
+    genes.push_back(make_pair(1., base));
+}
+
+bool genes_ord(pair<double, vector<permutation>> &a, pair<double, vector<permutation>> &b) {
+  return a.first < b.first;
+}
+
+void genetic_search::_run(int iterations) {
+  for(int t = 0; t < iterations; t++) {
+    // generate lambda new genes
+    for(int i = 0; i < lambda; i++) {
+      vector<permutation> &gene =  genes[mu+i].second;
+      if((double)rand()/RAND_MAX < c) {
+        // crossover
+        int g1 = rand()%mu, g2 = rand()%mu;
+        for(int d = 0; d < dim; d++)
+          gene[d] = permutation_crossover(genes[g1].second[d], genes[g2].second[d]);
+      }
+      else {
+        // mutation from a gene
+        gene = genes[rand()%mu].second;
+        int id = rand()%dim;
+        int p = ha.get_p(id);
+        int t1 = 1+rand()%(p-1), t2 = 1+rand()%(p-2);
+        // ensure t1 != t2
+        if(t2 >= t1)
+          t2++;
+        gene[id].transpose(t1, t2);
+      }
+    }
+    // for all gene
+    for(int i = 0; i < mu+lambda; i++) {
+      // replace permutations
+      // compute points
+      ha.set_pis(genes[i].second);
+      compute();
+      // compute discrepancy and check if it is good
+      check();
+      genes[i].first = current;
+    }
+    // sort genes : we want the mu firsts at the beginning
+    sort(genes.begin(), genes.end(), genes_ord);
+  }
+}

src/search.h

@@ -4,6 +4,8 @@
 #include "halton.h"
 #include "pointset.h"
 #include <string>
+#include <vector>
+#include <utility> // pair
 
 class search {
   public:
@@ -37,8 +39,10 @@ class local_search: public search {
     virtual bool accept(double previous, double current) = 0; // return true if we must accept this solution
   protected:
     void init(); // random initial configuration
-    virtual void random_neighbour();
-    virtual void undo(); // (no history : only one modification is remembered)
+    void random_neighbour();
+    virtual void _random_neighbour();
+    void undo(); // (no history : only one modification is remembered)
+    virtual void _undo();
   private:
     double previous;
     // for undo :
@@ -62,4 +66,16 @@ class sa_local_search: public local_search {
     double temp;
 };
 
+// mu+lambda genetic algorithm
+// c : proba of making a crossover
+class genetic_search: public search {
+  public:
+    genetic_search(int dim, int npoints, int *p, int mu, int lamda, double c);
+    virtual void _run(int iteartions);
+  protected:
+    int mu, lambda;
+    double c;
+    std::vector<std::pair<double, std::vector<permutation>>> genes; // discrepancy / permutations
+};
+
 #endif // SEARCH_H