Skip to content
Snippets Groups Projects
Commit b3429f94 authored by SIMONCINI DAVID's avatar SIMONCINI DAVID
Browse files

Initial commit

parents
No related branches found
No related tags found
No related merge requests found
Showing
with 1988 additions and 0 deletions
This diff is collapsed.
Compile:
mkdir build
cd build
cmake ../src/exe
make
Execute: Assuming 1aqt_0001.cfn is in data directory
./build/pils ./data/2erw_0001.cfn S "0 0.3 150" 500 60000 2erw-stats.txt >> 2erw-res.txt
S is an integer used as random seed
"0 0.3 150" are perturbation parameters: 0 means static perturbation, 0.3 is the size of the perturbation (proportional to sequence length) and 150 is the number of iterations allowed without improvement
500 is the maximum number of evaluations during each steepest descent
60000 is the maximum total number of evaluations.
1aqt-stats.txt contains statistics about the run
Our dataset is not included, feel free to contact us.
There is one sample instance in data subdirectory: unxz 2erw_0001.cfn.xz before running pils
You can generate your own instances using pompd :
https://forgemia.inra.fr/thomas.schiex/pompd
You will need to unzip the cfn.gz file generated with pompd
File added
cmake_minimum_required(VERSION 2.8)
######################################################################################
### 1) Set your application properties
######################################################################################
# Define your project name
PROJECT(fssp)
######################################################################################
# to have an optimised code
SET(CMAKE_CXX_FLAGS "-O3 -std=c++11")
#SET(CMAKE_CXX_FLAGS "-O3")
######################################################################################
### 2) Include the sources
######################################################################################
INCLUDE_DIRECTORIES(../../src)
######################################################################################
### 3) Define your target: just an executable here
#####################################################################################
ADD_EXECUTABLE(random_sampling random_sampling.cpp)
ADD_EXECUTABLE(firstImprovement_hc firstImprovement_hc.cpp)
ADD_EXECUTABLE(bestImprovement_hc bestImprovement_hc.cpp)
ADD_EXECUTABLE(doubleIncr_bestImprHC doubleIncr_bestImprHC.cpp)
ADD_EXECUTABLE(ils ils.cpp)
ADD_EXECUTABLE(readSolution readSolution.cpp)
ADD_EXECUTABLE(xsearch xsearch.h)
######################################################################################
### 4) Link the librairies for your executable
######################################################################################
#TARGET_LINK_LIBRARIES()
#ifndef __adaptivePerturb_h
#define __adaptivePerturb_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
//#include <algorithm>
#include <math.h>
#include <algo/perturbation.h>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
/*
Perturbation that modifies strength variables at random
strength is adaptive: increase when remains in the same local optima, decrease when new local optima
*/
class AdaptivePerturb : public Perturbation {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
AdaptivePerturb(std::mt19937 & _rng,
CostFunction & _eval,
NeighborEval & _neighborEval,
double _minStrength, double _maxStrength,
double _incrFactor, double _decrFactor) : Perturbation(_eval), rng(_rng), neighborEval(_neighborEval),
minStrength(_minStrength), maxStrength(_maxStrength),
incrFactor(_incrFactor), decrFactor(_decrFactor) {
ivar.resize(eval.size());
for(unsigned k = 0; k < ivar.size(); k++)
if (eval.n_values[k] > 1)
ivar[k] = k;
if (minStrength <= 0)
minStrength = 2;
if (maxStrength <= 0)
maxStrength = eval.size();
}
virtual void operator()(Solution & _solution) {
std::pair<int, int> neighbor;
double neighFit;
unsigned r;
unsigned n_strength;
// rounded the strength
double remains = strength - floor(strength);
n_strength = int (floor(strength));
if (remains > 0) {
if (runifd(rng) > remains)
n_strength++;
}
for(unsigned k = 0; k < n_strength; k++) {
// rnd neighbor
r = runif(rng) % (ivar.size() - k);
neighbor.first = ivar[r];
neighbor.second = (eval.n_values[ neighbor.first ]>1)?1 + runif(rng) % (eval.n_values[ neighbor.first ]-1):1;
// fitness of the neighbor
neighFit = neighborEval(_solution, neighbor);
nEval++;
// move
_solution.fitness(neighFit);
_solution[ neighbor.first ] = (_solution[ neighbor.first ] + neighbor.second) % eval.n_values[ neighbor.first ];
// "remove" the variable
std::swap(ivar[r], ivar[ivar.size() - 1 - k]);
}
}
virtual void init(Solution & _solution) {
strength = minStrength;
id = std::to_string(strength) + " ";
}
virtual void update(Solution & _previous, Solution & _solution) {
if (_previous.fitness() == _solution.fitness())
strength *= incrFactor;
else
strength *= decrFactor;
if (strength < minStrength)
strength = minStrength;
else
if (maxStrength < strength)
strength = maxStrength;
id = std::to_string(strength) + " ";
}
protected:
std::mt19937 & rng;
std::uniform_int_distribution<int> runif;
std::uniform_real_distribution<double> runifd;
// neighborhood
std::vector< unsigned > ivar;
// perturbation strength
double strength;
double minStrength, maxStrength;
double incrFactor, decrFactor;
NeighborEval & neighborEval;
};
#endif
#ifndef __bestImprovementHC_h
#define __bestImprovementHC_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
class BestImprovementHC : public LocalSearch {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
using LocalSearch::id;
BestImprovementHC(std::mt19937 & _rng,
CostFunction & _eval,
std::vector< std::pair<int, int> > & _neighborhood, NeighborEval & _neighborEval,
unsigned long long _nEvalMax,
std::ostream & _out) : LocalSearch(_eval), rng(_rng), neighborhood(_neighborhood), neighborEval(_neighborEval), nEvalMax(_nEvalMax), out(_out) {
}
virtual void operator()(Solution & _solution) {
unsigned long long nEvalLocal = 0;
unsigned long long nAccept = 0;
int r, n;
double neighFit, bestFit;
bool accept = true;
unsigned nBest;
std::vector<unsigned> iBest(neighborhood.size());
//out << id << nEval << " 0 " << _solution << std::endl;
while (nEvalLocal < nEvalMax && accept) {
out << id << nEvalLocal << " " << _solution << std::endl;
// find best neighbor
n = 0;
bestFit = neighborEval(_solution, neighborhood[n]);
nBest = 1;
iBest[0] = 0;
n++;
while (n < neighborhood.size() && nEvalLocal < nEvalMax) {
// test the neighbor
neighFit = neighborEval(_solution, neighborhood[n]);
// update best solution
if (neighFit < bestFit) {
bestFit = neighFit;
nBest = 1;
iBest[0] = n;
} else
if (neighFit == bestFit) { // when ties
nBest++;
iBest[nBest - 1] = n;
}
n++;
}
if (bestFit <= _solution.fitness()) {
accept = true;
if (nBest > 1) {
r = iBest[runif(rng) % nBest];
} else {
r = iBest[0];
}
// move
_solution.fitness(bestFit);
_solution[ neighborhood[r].first ] = (_solution[ neighborhood[r].first ] + neighborhood[r].second) % eval.n_values[ neighborhood[r].first ];
// out << id << nEval << " " << nAccept << " " << _solution << std::endl;
} else
accept = false;
nEvalLocal++;
}
nEval += nEvalLocal;
out << id << nEvalLocal << " " << _solution << std::endl;
}
protected:
std::mt19937 & rng;
std::uniform_int_distribution<int> runif;
unsigned long long nEvalMax;
// output
std::ostream & out;
// neighborhood
std::vector< std::pair<int, int> > neighborhood;
NeighborEval & neighborEval;
};
#endif
#ifndef __doubleIncr_biHC_h
#define __doubleIncr_biHC_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
#include <base/incrNeighborEval.h>
#include <base/bst.h>
class DoubleIncr_biHC : public LocalSearch {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
using LocalSearch::id;
DoubleIncr_biHC(std::mt19937 & _rng,
CostFunction & _eval,
IncrNeighborEval & _neighborEval,
unsigned long long _nEvalMax,
std::ostream & _out) : LocalSearch(_eval), rng(_rng), neighborEval(_neighborEval), nEvalMax(_nEvalMax), out(_out) {
delta.resize(eval.size());
neighborhoodSize = 0;
for(unsigned i = 0; i < delta.size(); i++) {
delta[i].resize(eval.n_values[i] - 1);
neighborhoodSize += eval.n_values[i] - 1;
}
}
void printDelta() {
for(unsigned i = 0; i < delta.size(); i++) {
for(unsigned k = 0; k < delta[i].size(); k++) {
std::cout << i << " " << (k+1) << " " << delta[i][k] << std::endl;
}
}
}
virtual void operator()(Solution & _solution) {
unsigned long long nEvalLocal = 0;
btree.makeEmpty();
bool accept = true;
unsigned nBest;
double bestDelta;
std::pair<int, int> bestNeigh;
std::vector< Node* > iBest(neighborhoodSize);
unsigned r, old_value;
while (nEvalLocal < nEvalMax && accept) {
if (nEvalLocal == 0) {
// first iteration:
// compute the initial deltas
std::pair<int, int> neighbor(0, 1);
for(unsigned i = 0; i < delta.size(); i++) {
for(unsigned k = 0; k < delta[i].size(); k++) {
neighbor.first = i;
neighbor.second = k + 1;
delta[i][k] = neighborEval(_solution, neighbor) - _solution.fitness();
// insert in the search tree
btree.insert(delta[i][k], neighbor);
}
}
} else {
// only update delta according to the accepted move: iBest[r]
unsigned i_move = iBest[r]->neighbor.first;
std::pair<int, int> neighbor(i_move, 1);
// update indice i_move with incremental evaluation (except back_move)
unsigned back_move = eval.n_values[i_move] - iBest[r]->neighbor.second - 1; // minus one because indice from 0
unsigned back_val = iBest[r]->neighbor.second - 1;
neighbor.second = back_move + 1;
// remove old value
btree.remove(delta[i_move][ back_move ], neighbor);
// compute new value
delta[i_move][ back_move ] = - delta[i_move][ back_val ];
// save new value
btree.insert(delta[i_move][ back_move ], neighbor);
unsigned k;
for(k = 0; k < back_move; k++) {
neighbor.second = k + 1;
// remove old value
btree.remove(delta[i_move][k], neighbor);
// compute new value
delta[i_move][k] = neighborEval(_solution, neighbor) - _solution.fitness();
// save new value
btree.insert(delta[i_move][k], neighbor);
}
for(k = back_move + 1; k < delta[i_move].size(); k++) {
neighbor.second = k + 1;
// remove old value
btree.remove(delta[i_move][k], neighbor);
// compute new value
delta[i_move][k] = neighborEval(_solution, neighbor) - _solution.fitness();
// save new value
btree.insert(delta[i_move][k], neighbor);
}
// update indice i linked to i_move with i_move < i
unsigned move_value, i;
for(unsigned kk = 0; kk < eval.links[i_move].size(); kk++) {
i = eval.links[i_move][ kk ];
neighbor.first = i;
for(k = 0; k < delta[ i ].size(); k++) {
move_value = (_solution[i] + k + 1) % eval.n_values[i];
neighbor.second = k + 1;
// remove old value
btree.remove(delta[i][k], neighbor);
// compute new value
delta[i][k] += eval.energy2[i_move][i][ old_value ][ _solution[i] ] - eval.energy2[i_move][i][ _solution[i_move] ][ _solution[i] ]
- eval.energy2[i_move][i][ old_value ][ move_value ] + eval.energy2[i_move][i][ _solution[i_move] ][ move_value ];
// save new value
btree.insert(delta[i][k], neighbor);
}
}
// update indice i linked to i_move with i < i_move
for(unsigned kk = 0; kk < neighborEval.interactions[i_move].size(); kk++) {
i = neighborEval.interactions[i_move][ kk ];
neighbor.first = i;
for(k = 0; k < delta[i].size(); k++) {
move_value = (_solution[i] + k + 1) % eval.n_values[i];
neighbor.second = k + 1;
// remove old value
btree.remove(delta[i][k], neighbor);
// compute new value
delta[i][k] += eval.energy2[i][i_move][ _solution[i] ][ old_value ] - eval.energy2[i][i_move][ _solution[i] ][ _solution[i_move] ]
+ eval.energy2[i][i_move][ move_value ][ _solution[i_move] ] - eval.energy2[i][i_move][ move_value ][ old_value ];
// save new value
btree.insert(delta[i][k], neighbor);
}
}
}
// find best neighbor
std::pair<int, int> bestNeigh;
btree.minimum(bestDelta, bestNeigh);
// find all best minima
btree.findall(bestDelta, nBest, iBest);
// Selection: select one of the best (if th(ere is ties)
if (bestDelta <= 0) {
accept = true;
if (nBest > 1) {
r = runif(rng) % nBest;
} else {
r = 0;
}
// save value before move (for double incr eval)
old_value = _solution[ iBest[r]->neighbor.first ];
// move
_solution.fitness(_solution.fitness() + bestDelta);
_solution[ iBest[r]->neighbor.first ] = (_solution[ iBest[r]->neighbor.first ] + iBest[r]->neighbor.second) % eval.n_values[ iBest[r]->neighbor.first ];
} else
accept = false;
nEvalLocal++;
}
nEval += nEvalLocal;
}
protected:
std::mt19937 & rng;
std::uniform_int_distribution<int> runif;
unsigned long long nEvalMax;
// output
std::ostream & out;
// incremental delta
unsigned neighborhoodSize;
std::vector< std::vector<double> > delta;
// neighborhood
IncrNeighborEval & neighborEval;
// to compute and update the minimum with BST
BST btree;
};
#endif
#ifndef __firstImprovementHC_h
#define __firstImprovementHC_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
class FirstImprovementHC : public LocalSearch {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
FirstImprovementHC(std::mt19937 & _rng,
CostFunction & _eval,
std::vector< std::pair<int, int> > & _neighborhood, NeighborEval & _neighborEval,
unsigned long long _nEvalMax,
std::ostream & _out, std::string & _id) : LocalSearch(_eval), rng(_rng), neighborhood(_neighborhood), neighborEval(_neighborEval), nEvalMax(_nEvalMax), out(_out), id(_id) {
}
virtual void operator()(Solution & _solution) {
unsigned long long nEvalLocal = 0;
unsigned long long nAccept = 0;
int r, n;
double neighFit;
bool accept = true;
//out << id << nEval << " 0 " << _solution << std::endl;
while (nEvalLocal < nEvalMax && accept) {
n = 0;
accept = false;
while (n < neighborhood.size() && !accept && nEval < nEvalMax) {
// random neighbor
r = runif(rng) % (neighborhood.size() - n);
// test the neighbor
neighFit = neighborEval(_solution, neighborhood[r]);
nEvalLocal++;
nEval++;
// acceptance
if (neighFit <= _solution.fitness()) {
accept = true;
nAccept++;
_solution.fitness(neighFit);
_solution[ neighborhood[r].first ] = (_solution[ neighborhood[r].first ] + neighborhood[r].second) % eval.n_values[ neighborhood[r].first ];
/*
if (neighFit < _solution.fitness()) {
out << id << nEval << " " << nAccept << " " << _solution << std::endl;
}
*/
} else {
std::swap(neighborhood[r], neighborhood[neighborhood.size() - 1 - n]);
}
n++;
}
}
}
protected:
std::mt19937 & rng;
std::uniform_int_distribution<int> runif;
unsigned long long nEvalMax;
// output
std::ostream & out;
// external id of the hc
std::string & id;
// neighborhood
std::vector< std::pair<int, int> > neighborhood;
NeighborEval & neighborEval;
};
#endif
#ifndef __ils_h
#define __ils_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
#include <base/incrNeighborEval.h>
#include <algo/perturbation.h>
using namespace std;
class ILS : public LocalSearch {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
using LocalSearch::id;
ILS(CostFunction & _eval,
IncrNeighborEval & _neighborEval,
LocalSearch & _ls,
Perturbation & _perturbation,
unsigned long long _nEvalMax,
unsigned _flatMax,
std::ostream & _out) : LocalSearch(_eval), neighborEval(_neighborEval),
ls(_ls), perturbation(_perturbation), nEvalMax(_nEvalMax), flatMax(_flatMax) , out(_out) {
}
virtual void operator()(Solution & _solution) {
unsigned i = 0;
id = std::to_string(i) + " " ;
perturbation.init(_solution);
// first ls
out << id << perturbation.id << nEval << " " << _solution << std::endl;
ls.id = std::to_string(i) + " ";
ls.nEval = 0;
ls(_solution);
nEval += ls.nEval;
//checkfit(_solution);
out << id << perturbation.id << nEval << " " << _solution << std::endl;
out << id << perturbation.id << nEval << " " << _solution << std::endl;
Solution previous;
i++;
unsigned flat = 0;
while (nEval < nEvalMax && flat < flatMax) {
if ( abs(_solution.fitness() - lastFit) < 10e-6 ){
flat++;
} else {
flat=0;
}
lastFit = _solution.fitness();
//std::cout << i << std::endl;
id = std::to_string(i) + " " ;
// copy the solution before perturbation
previous = _solution;
// perturbation
perturbation.nEval = 0;
perturbation(_solution);
nEval += perturbation.nEval;
//checkfit(_solution);
out << id << perturbation.id << nEval << " " << _solution << std::endl;
// local search
ls.nEval = 0;
ls.id = std::to_string(i) + " ";
ls(_solution);
nEval += ls.nEval;
// checkfit(_solution);
out << id << perturbation.id << nEval << " " << _solution << std::endl;
// update perturbation
perturbation.update(previous, _solution);
// selection
if (previous.fitness() < _solution.fitness()) {
// then previous solution is better, so restart from that solution
_solution = previous;
}
//checkfit(_solution);
out << id << perturbation.id << nEval << " " << _solution << std::endl;
i++;
}
}
void checkfit(Solution & _x) {
double f = _x.fitness();
eval(_x);
std::cout << f << " " << _x.fitness() << std::endl;
if (std::abs(f - _x.fitness()) > 1e-6) {
std::cout << "error fit" << std::endl;
//exit(1);
}
}
protected:
unsigned long long nEvalMax;
unsigned flatMax;
double lastFit = 0;
// output
std::ostream & out;
// local search (here HC)
LocalSearch & ls;
// perturbation
Perturbation & perturbation;
// neighborhood
IncrNeighborEval & neighborEval;
};
#endif
#ifndef _localSearch_h
#define _localSearch_h
#include <string.h>
#include <base/solution.h>
#include <base/costFunction.h>
class LocalSearch {
public:
LocalSearch(CostFunction & _eval) : eval(_eval) {
}
virtual void operator()(Solution & _solution) = 0;
// number of evaluation
unsigned long long nEval;
// id of the search (for output)
std::string id;
protected:
CostFunction & eval;
};
#endif
#ifndef __perturbation_h
#define __perturbation_h
#include <iostream>
#include <string>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
/*
Perturbation that modifies solution at random
*/
class Perturbation : public LocalSearch {
public:
Perturbation(CostFunction & _eval) : LocalSearch(_eval) {
}
virtual void init(Solution & _solution) {
}
virtual void update(Solution & _previous, Solution & _solution) {
}
};
#endif
/*
random_sampling.h
Author:
Sebastien Verel,
Univ. du Littoral Côte d'Opale, France.
David Simoncini
Univ. Toulouse 1 Capitole, France.
*/
#ifndef __randomSampling__h
#define __randomSampling__h
#include <vector>
#include <base/solution.h>
#include <init/init.h>
class RandomSampling {
public:
RandomSampling(CostFunction & _eval, Init & _init, unsigned _n) : eval(_eval), init(_init), n(_n) {
}
// run the sampling
void operator()() {
Solution x(eval.size());
for(unsigned i = 0; i < n; i++) {
init(x);
eval(x);
//x.print(); std::cout << std::endl;
sample.push_back( x.fitness() );
}
}
// fitness values of the sample
std::vector<double> sample;
protected:
// Automata
CostFunction & eval;
// initiatization
Init & init;
// size of the sample
unsigned n;
};
#endif
#ifndef __rndPerturb_h
#define __rndPerturb_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <math.h>
#include <algo/perturbation.h>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
/*
Perturbation that modifies strength variables at random
strength is random: between min and max values at each iteration
*/
class RandomPerturb : public Perturbation {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
RandomPerturb(std::mt19937 & _rng,
CostFunction & _eval,
NeighborEval & _neighborEval,
unsigned _minStrength, unsigned _maxStrength) : Perturbation(_eval), rng(_rng), neighborEval(_neighborEval),
minStrength(_minStrength), maxStrength(_maxStrength) {
ivar.resize(eval.size());
for(unsigned k = 0; k < ivar.size(); k++)
if (eval.n_values[k] > 1)
ivar[k] = k;
runifStrength.param( std::uniform_int_distribution<int>::param_type(minStrength, maxStrength) );
}
virtual void operator()(Solution & _solution) {
std::pair<int, int> neighbor;
double neighFit;
unsigned r;
for(unsigned k = 0; k < strength; k++) {
// rnd neighbor
r = runif(rng) % (ivar.size() - k);
neighbor.first = ivar[r];
neighbor.second = (eval.n_values[ neighbor.first ]>1)?1 + runif(rng) % (eval.n_values[ neighbor.first ]-1):1;
// fitness of the neighbor
neighFit = neighborEval(_solution, neighbor);
nEval++;
// move
_solution.fitness(neighFit);
_solution[ neighbor.first ] = (_solution[ neighbor.first ] + neighbor.second) % eval.n_values[ neighbor.first ];
// "remove" the variable
std::swap(ivar[r], ivar[ivar.size() - 1 - k]);
}
}
virtual void init(Solution & _solution) {
strength = runifStrength(rng);
id = std::to_string(strength) + " ";
}
virtual void update(Solution & _previous, Solution & _solution) {
strength = runifStrength(rng);
id = std::to_string(strength) + " ";
}
protected:
std::mt19937 & rng;
std::uniform_int_distribution<int> runif;
std::uniform_int_distribution<int> runifStrength;
// neighborhood
std::vector< unsigned > ivar;
// perturbation strength
unsigned strength;
unsigned minStrength, maxStrength;
NeighborEval & neighborEval;
};
#endif
import os
from pylab import *
import matplotlib.pyplot as plt
import numpy as numpy
statFile = open("../../build/stats.txt", 'r')
statData =statFile.read()
enP1=[]
enP2=[]
nEval=[]
nbCompCo=[]
nbIt=[]
It=0
for ligne in statFile:
enP1.append( int(ligne.split[0]) )
enP2.append( int(ligne.split[1]) )
nEval.append( int(ligne.split[2]) )
nbCompCo.append( int(ligne.split[3]) )
nbIT.append(It)
It+1
nbIt=linspace(0, len(enP1), len(enP1))
plt.figure()
plt.plot(enP1,nbIt,'r', label='parent 1')
plt.plot(enP2,nbIt,'b', label='parent 2')
plt.plot(nEval,nbIt,'g', label='number of evaluation')
plt.plot(nbCompCo,nbIt,'y', label='number of connexes components')
plt.legend()
plt.xlabel('number of iterations')
plt.show()
"""
plt.clf()
x=np.linspace(nbIter, 1000)
plt.show()
"""
stat.close()
#ifndef __staticPerturb_h
#define __staticPerturb_h
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <algo/perturbation.h>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
/*
Perturbation that modifies _strength variables at random
it uses incremental evaluation neighborEval at each modification
*/
class StaticPerturb : public Perturbation {
public:
using LocalSearch::nEval;
using LocalSearch::eval;
StaticPerturb(std::mt19937 & _rng,
CostFunction & _eval,
NeighborEval & _neighborEval,
unsigned _strength) : Perturbation(_eval), rng(_rng), neighborEval(_neighborEval), strength(_strength) {
ivar.resize(eval.size());
for(unsigned k = 0; k < ivar.size(); k++)
if (eval.n_values[k] > 1)
ivar[k] = k;
}
virtual void operator()(Solution & _solution) {
std::pair<int, int> neighbor;
double neighFit;
unsigned r;
for(unsigned k = 0; k < strength; k++) {
// rnd neighbor
r = runif(rng) % (ivar.size() - k);
neighbor.first = ivar[r];
neighbor.second = (eval.n_values[ neighbor.first ]>1)?1 + runif(rng) % (eval.n_values[ neighbor.first ]-1):1;
// fitness of the neighbor
neighFit = neighborEval(_solution, neighbor);
nEval++;
// move
_solution.fitness(neighFit);
_solution[ neighbor.first ] = (_solution[ neighbor.first ] + neighbor.second) % eval.n_values[ neighbor.first ];
// "remove" the variable
std::swap(ivar[r], ivar[ivar.size() - 1 - k]);
}
}
virtual void init(Solution & _solution) {
id = std::to_string(strength) + " ";
}
protected:
std::mt19937 & rng;
std::uniform_int_distribution<int> runif;
// neighborhood
std::vector< unsigned > ivar;
// perturbation strength
unsigned strength;
NeighborEval & neighborEval;
};
#endif
/*
#ifndef __ils_h
#define __ils_h
*/
#include <random>
#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <map>
#include <algo/localSearch.h>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/neighborEval.h>
#include <base/incrNeighborEval.h>
#include <algo/perturbation.h>
#include <base/Xover.h>
#include <init/init.h>
using namespace std;
class Xsearch : public LocalSearch {
public :
using LocalSearch::nEval;
using LocalSearch::eval;
using LocalSearch::id;
Xsearch(CostFunction & _eval,
IncrNeighborEval & _neighborEval,
LocalSearch & _ls1,
Perturbation & _pertu1,
unsigned long long _nEvalMax,
unsigned _flatMax,
ostream & statFlux) : LocalSearch(_eval), neighborEval(_neighborEval), ls1(_ls1), pertu1(_pertu1),
nEvalMax(_nEvalMax), flatMax(_flatMax), statF(statFlux) {
}
virtual void operator()(Solution & sol1){
unsigned i=0;
Xover xo(eval, statF);
pertu1.init(sol1);
ls1.nEval = 0;
ls1(sol1);
nEval += ls1.nEval;
unsigned flat = 0;
while(nEval < nEvalMax && flat < flatMax){
if ( abs(sol1.fitness() - lastFit) < 10e-6 ){
flat++;
} else {
flat=0;
}
lastFit = sol1.fitness();
sol2 = sol1;
pertu1.nEval =0;
pertu1(sol2);
nEval += pertu1.nEval;
ls1.nEval = 0;
ls1(sol2);
nEval += ls1.nEval;
statF << sol1.fitness() << " " << sol2.fitness() << " " << nEval << " ";
xo(sol1, sol2, sol1);
ls1.nEval = 0;
ls1(sol1);
nEval += ls1.nEval;
}
}
protected :
unsigned long long nEvalMax;
unsigned flatMax;
double lastFit = 0;
bool same;
bool growth;
Solution sol2;
ostream & statF;
LocalSearch & ls1;
Perturbation & pertu1;
IncrNeighborEval & neighborEval;
multimap<double, Solution> fitnesss;
multimap<double, Solution>::iterator it1;
};
/*
xover.h
*/
//#ifndef _Xover_h
//#define _Xover_h
#include <string>
#include <iostream>
#include <streambuf>
#include <stdlib.h>
#include <set>
#include <base/solution.h>
#include <base/costFunction.h>
#include <base/crossover.h>
using namespace std;
class Xover : public Crossover {
public :
Xover(CostFunction & _eval, ostream & _statF) : eval(_eval), statF(_statF) {
nComponents = 0;
component_id.resize(eval.n_variables);
components.resize(eval.n_variables + 1);
}
bool operator()(Solution & x1,Solution & x2, Solution & child) {
child.resize(x1.size());
compute_connected_components(x1, x2);
double g1, g2;
bool n1 = false;
bool n2 = false;
statF << nComponents << " ";
// first common component : copy x1 into child
double child_fit = 0;
child_fit = partial_func(0, x1) ;
for(unsigned i : components[0])
child[i] = x1[i];
// Others true components
for(unsigned c = 1; c <= nComponents; c++) {
g1 = partial_func(c, x1);
g2 = partial_func(c, x2);
if (g1 >= g2) {
// x2 is better for that component
n2 = true;
child_fit += g2;
for(unsigned i : components[c])
child[i] = x2[i];
} else {
// x1 is better for that component
n1 = true;
child_fit += g1;
for(unsigned i : components[c])
child[i] = x1[i];
}
}
child.fitness(child_fit);
statF << (n1 && n2) << endl ;
return (n1 && n2); // true when child is different from x1 and x2
}
void print_state() {
std::cout << endl << "ids:" ;
for(unsigned i : component_id) {
std::cout << " " << i ;
}
std::cout << std::endl;
std::cout << "nb components : " << nComponents << std::endl;
for(unsigned c = 0; c <= nComponents; c++) {
std::cout << c << ":" ;
for(unsigned i : components[c])
{
std::cout << " " << i << " links:" << endl ;
for (unsigned j : eval.links[i])
cout << j << " ";
cout << endl;
}
std::cout << std::endl;
}
std::cout << std::endl;
}
void afficheComponents() {
for(int i=0; i <= nComponents; i++){
cout << endl << "comp" << i << endl;
for(int j : components[i]){
cout << components[i][j] << " ";
}
}
cout << endl;
}
//number of connected components
unsigned nComponents;
//number of the connected components of each variables
std::vector<int> component_id;
//list of variables of each connected component
std::vector < std::vector<unsigned> > components;
//costfunction to evaluate the variables
CostFunction & eval;
ostream & statF;
protected:
void compute_connected_components (Solution & x1, Solution & x2) {
nComponents = 0;
std::fill(component_id.begin(), component_id.end(), -1);
for(unsigned i = 0; i < eval.n_variables; i++) // avoid to copy the references on x1 and x2 in set_component function
if (x1[i] == x2[i])
component_id[i] = 0; // in the share variable component
//for all variables
for (int i = 0; i < eval.n_variables; ++i)
{
if (component_id[i] == -1) {
nComponents++;
set_component(i);
}
}
for(unsigned c = 0; c <= nComponents; c++)
components[c].resize( 0 );
for(unsigned i = 0; i < eval.n_variables; i++)
components[ component_id[i] ].push_back(i);
}
void set_component(unsigned i) {
if (component_id[i] == -1) {
component_id[i] = nComponents;
for(unsigned j : eval.links[i]) // all neighbors in the variable links matrix
set_component(j);
for(unsigned j : eval.backlinks[i]) // check out backlinks too
set_component(j);
}
}
double partial_func(unsigned c, Solution & x){
double res= 0;
vector<bool> inthere;
inthere.resize(x.size());
fill(inthere.begin(), inthere.end(), false);
for(unsigned k : components[c])
{
inthere[k]=true;
}
for(unsigned k : components[c]) {
res+= eval.energy[k][x[k]]; // add the energy of each AA
for(unsigned i = 0; i < eval.links[k].size(); i++) {
unsigned l = eval.links[k][i];
if (inthere[l])
res += eval.energy2[k][l][x[k]][x[l]];
}
}
if (c!=0)
{
for(unsigned k : components[c]) {
for(unsigned i = 0; i < eval.links[k].size(); i++) {
unsigned l = eval.links[k][i];
if (!inthere[l])
res += eval.energy2[k][l][x[k]][x[l]];
}
for(unsigned i = 0; i < eval.backlinks[k].size(); i++) {
unsigned l = eval.backlinks[k][i];
if (!inthere[l])
res += eval.energy2[l][k][x[l]][x[k]];
}
}
}
return res;
}
};
#ifndef __BST_h
#define __BST_h
/*
Binary search Tree : AVL
*/
#include <iostream>
class Node {
public:
double fitness;
std::pair<int, int> neighbor;
int height;
Node * left;
Node * right;
Node(double _fitness, std::pair<int, int> & _neighbor) : fitness(_fitness), neighbor(_neighbor), height(1), left(NULL), right(NULL) {}
Node(Node * _source) : fitness(_source->fitness), neighbor(_source->neighbor), height(1), left(NULL), right(NULL) {}
~Node() {
}
void updateHeight() {
if (left == NULL) {
if (right == NULL)
height = 1;
else
height = right->height + 1;
} else
if (right == NULL)
height = left->height + 1;
else
height = std::max(left->height, right->height) + 1;
}
int balanceIndex() {
if (left == NULL) {
if (right == NULL)
return 0;
else
return - right->height;
} else
if (right == NULL)
return left->height ;
else
return left->height - right->height;
}
void print(std::ostream & _out) {
_out << "{" << height << ", " << fitness << " " << neighbor.first << " " << neighbor.second << "}";
}
};
class BST {
public:
BST() : root(NULL) { }
~BST() {
destroy(root);
}
void makeEmpty() {
destroy(root);
root = NULL;
}
void insert(double _fitness, std::pair<int, int> & _neighbor) {
root = insert(root, _fitness, _neighbor);
}
void remove(double _fitness, std::pair<int, int> & _neighbor) {
root = remove(root, _fitness, _neighbor);
}
void minimum(double & vmin, std::pair<int, int> & min) {
Node * n = minimum(root);
min.first = n->neighbor.first;
min.second = n->neighbor.second;
vmin = n->fitness;
}
void findall(double _fitness, unsigned & _nb, std::vector< Node* > & _res) {
_nb = 0;
findall(root, _fitness, _nb, _res);
}
void print(std::ostream & _out) {
prefixe(root, _out);
}
void copy(BST & _dest) {
_dest.root = copy(root);
}
void copy(Node * _source, BST & _dest) {
_dest.root = copy(_source);
}
bool checkHeight() {
if (root == NULL)
return true;
else
return (root->height == height(root));
}
bool checkBalance() {
if (root == NULL)
return true;
else
return checkBalance(root);
}
Node * root;
private:
void destroy(Node * n) {
if (n != NULL) {
destroy(n->left);
destroy(n->right);
delete n;
}
}
Node* copy(Node * n) {
if (n == NULL)
return n;
else {
Node * cn = new Node(n);
cn->left = copy(n->left);
cn->right = copy(n->right);
cn->height = n->height;
return cn;
}
}
unsigned height(Node * n) {
if (n == NULL)
return 0;
else {
unsigned l = height(n->left);
unsigned r = height(n->right);
return std::max(l, r) + 1;
}
}
bool checkBalance(Node * n) {
if (n == NULL)
return true;
else {
bool l = checkBalance(n->left);
bool r = checkBalance(n->right);
return l && r && (std::abs(n->balanceIndex()) < 2);
}
}
Node * minimum(Node * n) {
if (n == NULL)
return NULL;
else {
if (n->left == NULL)
return n;
else
return minimum(n->left);
}
}
void findall(Node * n, double fit, unsigned & nb, std::vector< Node* > & res) {
if (n != NULL) {
if (fit == n->fitness) {
res[nb] = n;
nb++;
findall(n->left, fit, nb, res);
findall(n->right, fit, nb, res);
} else
if (fit < n->fitness)
findall(n->left, fit, nb, res);
else
findall(n->right, fit, nb, res);
}
}
Node* insert(Node * n, double _fitness, std::pair<int, int> & _neighbor) {
if (n == NULL) {
n = new Node(_fitness, _neighbor);
} else {
if (_fitness < n->fitness ||
(_fitness == n->fitness && _neighbor.first < n->neighbor.first) ||
(_fitness == n->fitness && _neighbor.first == n->neighbor.first && _neighbor.second < n->neighbor.second)) {
n->left = insert(n->left, _fitness, _neighbor);
n->updateHeight();
n = balance(n);
} else {
n->right = insert(n->right, _fitness, _neighbor);
n->updateHeight();
n = balance(n);
}
}
return n;
}
Node* remove(Node * n, double _fitness, std::pair<int, int> & _neighbor) {
if (n != NULL) {
if (n->fitness == _fitness && n->neighbor.first == _neighbor.first && n->neighbor.second == _neighbor.second) {
// remove this node
if (n->left == NULL) {
Node * tmp = n;
n = balance(n->right);
//n = balance(n);
delete tmp;
} else {
std::pair<Node*, Node*> m = cutMax(n->left);
m.second->left = m.first;
m.second->right = n->right;
m.second->updateHeight();
delete n;
n = m.second;
}
} else
if (_fitness < n->fitness ||
(_fitness == n->fitness && _neighbor.first < n->neighbor.first) ||
(_fitness == n->fitness && _neighbor.first == n->neighbor.first && _neighbor.second < n->neighbor.second)) {
n->left = remove(n->left, _fitness, _neighbor);
n->updateHeight();
} else {
n->right = remove(n->right, _fitness, _neighbor);
n->updateHeight();
}
n = balance(n);
}
return n;
}
std::pair<Node*, Node*> cutMax(Node * n) {
if (n->right == NULL) {
return std::pair<Node*, Node*>(n->left, n);
} else {
std::pair<Node*, Node*> res = cutMax(n->right);
n->right = res.first;
n->updateHeight();
n = balance(n);
return std::pair<Node*, Node*>(n, res.second);
}
}
void prefixe(Node * n, std::ostream & _out) {
_out << "(";
if (n != NULL) {
n->print(_out);
prefixe(n->left, _out);
prefixe(n->right, _out);
}
_out << ")";
}
int heightNode(Node * n) {
return n ? n->height : 0;
}
Node* balance(Node * n) {
if (n != NULL) {
int b = n->balanceIndex();
if (std::abs(b) > 1) {
if (b == 2) {
if (heightNode(n->left->left) >= heightNode(n->left->right)) {
n = rotateRight(n);
} else {
n->left = rotateLeft(n->left);
n = rotateRight(n);
}
} else {
if (heightNode(n->right->right) >= heightNode(n->right->left)) {
n = rotateLeft(n);
} else {
n->right = rotateRight(n->right);
n = rotateLeft(n);
}
}
}
}
return n;
}
Node* rotateLeft(Node * p) {
// right
Node * rt = p->right;
// left of the right
p->right = rt->left;
p->updateHeight();
rt->left = p;
rt->updateHeight();
return rt;
}
Node* rotateRight(Node * p) {
// left
Node * lt = p->left;
p->left = lt->right;
p->updateHeight();
lt->right = p;
lt->updateHeight();
return lt;
}
};
#endif
\ No newline at end of file
/*
costFunction.h
Author:
Sebastien Verel,
Univ. du Littoral Côte d'Opale, France.
David Simoncini
Univ. Toulouse 1 Capitole, France.
*/
#ifndef _costFunction_h
#define _costFunction_h
#include <string>
#include <iostream>
#include <streambuf>
#include <base/solution.h>
#include "rapidjson/document.h"
#include "rapidjson/filereadstream.h"
#include <cstdio>
using namespace std;
/*
Basic Energy Function
*/
class CostFunction {
public:
CostFunction(const char * _instance_fileName) {
n_variables = 0;
rapidjson::Document document;
FILE* fp = fopen(_instance_fileName, "rb"); // non-Windows use "r"
if (!fp)
std::cerr << "Impossible to open " << _instance_fileName << std::endl;
// read first line which is a comment
int i, ch;
while ((ch = fgetc(fp)) != '\n' && ch != EOF) { i++; }
char readBuffer[65536];
rapidjson::FileReadStream is(fp, readBuffer, sizeof(readBuffer));
document.ParseStream(is);
fclose(fp);
if (document.IsObject()) {
// variables
if (document.HasMember("variables") && document["variables"].IsObject()) {
for (rapidjson::Value::ConstMemberIterator itr = document["variables"].MemberBegin(); itr != document["variables"].MemberEnd(); ++itr) {
std::string v = itr->name.GetString();
n_variables++;
// read possibles values
const rapidjson::Value& a = document["variables"][v.c_str()];
// number of possible values
n_values.push_back( a.Size() );
// save the name of the values
std::vector<std::string> names;
for (rapidjson::SizeType k = 0; k < a.Size(); k++) {
names.push_back( a[k].GetString() );
}
value_name.push_back(names);
}
}
// energy functions
if (document.HasMember("functions") && document["functions"].IsObject()) {
// linear part
energy.resize(n_variables);
char e[256];
for(unsigned i = 0; i < n_variables; i++) {
sprintf(e, "E%d", i+1);
const rapidjson::Value& a = document["functions"][e]["costs"];
for (rapidjson::SizeType k = 0; k < a.Size(); k++) {
energy[i].push_back(a[k].GetDouble());
}
}
// quadratic part
energy2.resize(n_variables);
links.resize(n_variables);
backlinks.resize(n_variables);
for(unsigned i = 0; i < n_variables; i++) {
energy2[i].resize(n_variables);
for(unsigned j = i + 1; j < n_variables; j++) {
sprintf(e, "E%d_%d", i+1, j+1);
rapidjson::Value::ConstMemberIterator itr = document["functions"].FindMember(e);
rapidjson::Value::ConstMemberIterator itr2;
if (itr != document["functions"].MemberEnd()) {
links[i].push_back(j);
backlinks[j].push_back(i);
// two possibilities
itr2 = itr->value.FindMember("defaultcost");
if (itr2 == itr->value.MemberEnd()) {
// without defaultcost:
const rapidjson::Value& costs = itr->value["costs"];
// quadratic term : energy2[i][j][a1][a2] : position i and j, and then acide/angle a1 in position i, and a2 in position j
energy2[i][j].resize(n_values[i]);
for (rapidjson::SizeType k = 0; k < costs.Size(); k++) {
if (k / n_values[j] < n_values[i])
energy2[i][j][k / n_values[j]].push_back(costs[k].GetDouble());
}
} else {
// with defaultcost
double defaultcost = itr2->value.GetDouble();
// without defaultcost:
const rapidjson::Value& costs = itr->value["costs"];
// quadratic term : energy2[i][j][a1][a2] : position i and j, and then acide/angle a1 in position i, and a2 in position j
energy2[i][j].resize(n_values[i]);
for(unsigned a1 = 0; a1 < n_values[i]; a1++)
energy2[i][j][a1].resize(n_values[j], defaultcost);
rapidjson::SizeType k = 0;
while (k < costs.Size()) {
energy2[i][j][ costs[k].GetInt() ][ costs[k + 1].GetInt() ] = costs[k + 2].GetDouble();
k += 3;
}
}
// check size
if (energy2[i][j].size() != n_values[i]) {
std::cout << "error instance " << i << " " << j << " " << n_values[i] << " " << energy2[i][j].size() << std::endl;
}
for(unsigned a = 0; a < n_values[i]; a++)
if (energy2[i][j][a].size() != n_values[j]) {
std::cout << "error instance " << i << " " << n_values[i] << " " << a << " " << energy2[i][j][a].size() << std::endl;
}
}
}
}
}
} else
std::cerr << "Impossible to parse the json file " << _instance_fileName << std::endl;
}
unsigned int size() const {
return n_variables;
};
unsigned int variable_size(unsigned k) const {
return n_values[k];
};
void operator()(Solution & _x) {
double fit = 0;
unsigned j;
for(unsigned i = 0; i < _x.size(); i++) {
fit += energy[i][ _x[i] ];
for(unsigned k = 0; k < links[i].size(); k++) {
j = links[i][k];
if (_x[i] < energy2[i][j].size() && _x[j] < energy2[i][j][_x[i]].size())
fit += energy2[i][j][ _x[i] ][ _x[j] ];
else {
if (_x[i] >= energy2[i][j].size())
std::cerr << "error i=" << i << " " << _x[i] << " " << energy2[i][j].size() << std::endl;
else
std::cerr << "error j=" << j << " " << _x[j] << " " << energy2[i][j][_x[i]].size() << std::endl;
}
}
}
_x.fitness(fit);
}
void print() {
std::cout << n_variables << std::endl;
for(unsigned i = 0; i < n_variables; i++) {
std::cout << i << " " << n_values[i] ;
for(unsigned j = 0; j < n_values[i]; j++) {
std::cout << " " << value_name[i][j];
}
std::cout << std::endl;
}
for(unsigned i = 0; i < n_variables; i++) {
std::cout << i ;
for(unsigned k = 0; k < links[i].size(); k++) {
std::cout << " " << links[i][k];
}
std::cout << std::endl;
}
for(unsigned i = 0; i < n_variables; i++) {
std::cout << i << " " << n_values[i] ;
for(unsigned j = 0; j < n_values[i]; j++) {
std::cout << " " << energy[i][j];
}
std::cout << std::endl;
}
for(unsigned i = 0; i < n_variables; i++) {
for(unsigned j = i+1; j < n_variables; j++) {
if (energy2[i][j].size() > 0) {
for(unsigned a1 = 0; a1 < energy2[i][j].size(); a1++)
for(unsigned a2 = 0; a2 < energy2[i][j][a1].size(); a2++)
std::cout << i << " " << j << " " << a1 << " " << a2 << " " << energy2[i][j][a1][a2] << std::endl;
}
}
}
}
void printOnShort(std::ostream& _os) const {
_os << n_variables << std::endl;
_os << n_values[0] ;
for(unsigned i = 1; i < n_variables; i++) {
_os << " " << n_values[i] ;
}
_os << std::endl;
// l_i
for(unsigned i = 0; i < n_variables; i++) {
_os << i ;
for(unsigned k = 0; k < links[i].size(); k++) {
if (i < links[i][k])
_os << " " << links[i][k] ;
}
_os << std::endl ;
}
}
void printlinks() {
for(int i; i<links.size() ; i ++)
{
cout << " links " << i << ": ";
for (unsigned j : links[i]){
cout << j << " ";
}
cout << endl;
}
cout << endl;
}
void printEnergy2() {
for(unsigned i = 0; i< 2; i ++)
{
cout << i << ": " << endl;
for (unsigned j =0; j< n_variables; j++){
cout << " " << j << ": " << endl ;
for(unsigned k = 0; k< energy2[i][j].size(); k++) {
for (unsigned l = 0; l< energy2[i][j][k].size(); l++){
cout << " " << k << " " << l << " " << energy2[i][j][k][l] << endl ;
}
}
cout << endl;
}
cout << endl;
}
}
//protected:
unsigned int n_variables;
// number of values for each variable
std::vector<unsigned> n_values;
// value names
std::vector< std::vector< std::string > > value_name;
// linear term
std::vector< std::vector<double> > energy;
// interaction between variables with i < j
std::vector< std::vector<unsigned> > links; // who am I connecting to ?
std::vector< std::vector<unsigned> > backlinks; // who is connecting to me ?
// quadratic term : energy2[i][j][a1][a2] : position i and j, and then acide/angle a1 in position i, and a2 in position j
std::vector< std::vector< std::vector< std::vector<double> > > > energy2;
};
#endif
#ifndef _crossover_h
#define _crossover_h
class Crossover {
public:
/*
return true when the child solution is different to both parent solutions (improvement)
*/
virtual bool operator()(Solution & x1, Solution & x2, Solution & child) = 0;
};
#endif
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment