projet_optim_m2/src/search.cpp

#include "search.h"
#include <cmath>
#include <cassert>
#include <algorithm> // sort
#include <tuple> // tie, ignore

using namespace std;

/**********************    search    ***********************/

search::search(int dim, int npoints, int *p) : dim(dim), npoints(npoints),
    best(1.), ha(dim, p), ps(dim, npoints) {
  assert(dim > 0);
  assert(npoints > 0);
}

void search::check() {
  current = ps.discrepancy();
  if(current < best) {
    best = current;
    FILE *f = fopen(filename.c_str(), "w");
    ps.dump(f);
    fclose(f);
  }
}

double search::run(int iterations) {
  _run(iterations);
  fprintf(stderr,"Best discrepancy found after %d iterations : %g\nResult saved in %s\n", iterations, best, filename.c_str());
  printf("%g\n", best);
  return best;
}

// compute points
void search::compute() {
  for(int i = 0; i < npoints; i++)
    ha.compute(i+1, ps.point(i));
}


/**********************    random_search    ***********************/

random_search::random_search(int dim, int npoints, int *p) : search(dim, npoints, p) {
  filename = std::string("random_search");
}

void random_search::_run(int iterations) {
  for(int t = 0; t < iterations; t++) {
    // shuffle permutations
    for(int i = 0; i < dim; i++)
      ha.get_pi(i).random();

    // compute points
    compute();
    // compute discrepancy and check if it is good
    check();
  }
}


/**********************    local_search    ***********************/

local_search::local_search(int dim, int npoints, int *p) : search(dim, npoints, p),
    undoable(false) {
}

void local_search::init() {
  // random initial configuration
  for(int i = 0; i < dim; i++)
    ha.get_pi(i).random();
  compute();
  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);
  t2 = 1+rand()%(p-2);
  // ensure t1 != t2
  if(t2 >= t1)
    t2++;
  ha.get_pi(i).transpose(t1, t2);
}

void local_search::undo() {
  assert(undoable);
  undoable = false;
  _undo();
}

void local_search::_undo() {
  ha.get_pi(i).transpose(t1, t2);
}

void local_search::_run(int iterations) {
  for(int t = 0; t < iterations; t++) {
    previous = current;
    random_neighbour();
    compute();
    check();
    if(!accept(previous, current)) {
      current = previous;
      undo();
    }
  }
}


/**********************    random_local_search    ***********************/

random_local_search::random_local_search(int dim, int npoints, int *p)
    :local_search(dim, npoints, p) {
  filename = std::string("random_local_search");
  init();
}

bool random_local_search::accept(double previous, double current) {
  return current < previous;
}


/**********************    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) {
  assert(0 <= lambda && lambda < 1);
  filename = std::string("sa_local_search_") + std::to_string(lambda);
  init();
}

bool sa_local_search::accept(double previous, double current) {
  if(current < previous)
    temp *= lambda;
  //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);
  }
}