Various changes and optimizations

This commit is contained in:
maartenkeijzer 2007-09-10 10:30:53 +00:00
commit 36c14fed9e
8 changed files with 366 additions and 145 deletions

View file

@ -3,6 +3,8 @@
#include <eoAlgo.h>
#include <moo/eoEpsilonArchive.h>
#include <utils/eoStat.h>
template <class EOT>
class eoEpsMOEA : public eoAlgo<EOT> {
@ -13,14 +15,16 @@ class eoEpsMOEA : public eoAlgo<EOT> {
eoContinue<EOT>& _continuator,
eoEvalFunc<EOT>& _eval,
eoGenOp<EOT>& _op,
const std::vector<double>& eps
const std::vector<double>& eps,
unsigned max_archive_size
) : continuator(_continuator),
eval (_eval),
loopEval(_eval),
popEval(loopEval),
op(_op),
archive(eps)
{}
archive(eps, max_archive_size)
{
}
void operator()(eoPop<EOT>& pop) {
@ -48,7 +52,11 @@ class eoEpsMOEA : public eoAlgo<EOT> {
}
}
} while (continuator(pop));
// quite expensive to copy the entire archive time and time again, but this will make it work more seamlessly with the rest of EO
offspring.clear();
archive.appendTo(offspring);
} while (continuator(offspring)); // check archive
// return archive
pop.clear();
@ -58,13 +66,23 @@ class eoEpsMOEA : public eoAlgo<EOT> {
private :
void update_pop(eoPop<EOT>& pop, const EOT& offspring) {
std::vector<unsigned> dominated;
for (unsigned i = 0; i < pop.size(); ++i) {
int dom = offspring.fitness().check_dominance(pop[i].fitness());
if (dom == 1) {
pop[i] = offspring;
return;
switch (dom) {
case 1 : // indy dominates this
dominated.push_back(i);
break;
case -1 : // is dominated, do not insert
return;
case 0: // incomparable
break;
}
if (dom == -1) return; //
}
if (dominated.size()) {
pop[ dominated[ rng.random(dominated.size()) ] ] = offspring;
}
// non-dominated everywhere, overwrite random one
@ -111,7 +129,7 @@ class eoEpsMOEA : public eoAlgo<EOT> {
eoGenOp<EOT>& op;
eoEpsilonArchive<EOT> archive;
};
#endif

View file

@ -11,15 +11,21 @@ class eoEpsilonArchive {
struct Node {
EOT element;
std::vector<double> discretized;
std::vector<long> discretized;
Node(const EOT& eo) : element(eo), discretized(Traits::nObjectives()) {}
dominance::dominance_result check(const Node& other) const {
return dominance::check_discrete(discretized, other.discretized);
}
};
typedef std::list<Node> archive_t;
typedef std::vector<Node> archive_t;
archive_t archive;
std::vector<double> eps;
std::vector<double> inv_eps;
unsigned max_size;
public:
@ -27,8 +33,12 @@ class eoEpsilonArchive {
static double tol() { return 1e-6; }
eoEpsilonArchive(const std::vector<double>& eps_) : eps(eps_) {
if (eps.size() != Traits::nObjectives()) throw std::logic_error("eoEpsilonArchive: need one epsilon for each objective");
eoEpsilonArchive(const std::vector<double>& eps_, unsigned max_size_ = std::numeric_limits<unsigned>::max()) : max_size(max_size_) {
inv_eps.resize(eps_.size());
for (unsigned i = 0; i < inv_eps.size(); ++i) {
inv_eps[i] = 1.0 / eps_[i];
}
if (inv_eps.size() != Traits::nObjectives()) throw std::logic_error("eoEpsilonArchive: need one epsilon for each objective");
}
bool empty() { return archive.size() == 0; }
@ -43,29 +53,31 @@ class eoEpsilonArchive {
for (unsigned i = 0; i < eo.fitness().size(); ++i) {
double val = Traits::direction(i) * eo.fitness()[i];
node.discretized[i] = floor(val/eps[i]);
node.discretized[i] = (long) floor(val*inv_eps[i]);
}
using namespace dominance;
typename archive_t::iterator boxIt = archive.end();
unsigned box = archive.size();
// update archive
for (typename archive_t::iterator it = archive.begin(); it != archive.end(); ++it) {
dominance_result result = check<eoEpsilonArchive<EOT> >(node.discretized, it->discretized);
for (unsigned i = 0; i != archive.size(); ++i) {
dominance_result result = node.check(archive[i]); //check<eoEpsilonArchive<EOT> >(node.discretized, archive[i].discretized);
switch (result) {
case first : { // remove dominated archive member
it = archive.erase(it);
std::swap( archive[i], archive.back());
archive.pop_back();
--i;
break;
}
case second : {
//cout << it->discretized[0] << ' ' << it->discretized[1] << " dominates " << node.discretized[0] << ' ' << node.discretized[1] << endl;
return; // new one does not belong in archive
}
case non_dominated : break; // both non-dominated, put in different boxes
case non_dominated_equal : { // in same box
boxIt = it;
box = i;
goto exitLoop; // break
}
}
@ -74,33 +86,60 @@ class eoEpsilonArchive {
exitLoop:
// insert
if (boxIt == archive.end()) { // non-dominated, new box
if (box >= archive.size()) { // non-dominated, new box
archive.push_back(node);
} else { // fight it out
int dom = node.element.fitness().check_dominance( boxIt->element.fitness() );
int dom = node.element.fitness().check_dominance( archive[box].element.fitness() );
switch (dom) {
case 1: *boxIt = node; break;
case 1: archive[box] = node; break;
case -1: break;
case 0: {
double d1 = 0.0;
double d2 = 0.0;
for (unsigned i = 0; i < node.element.fitness().size(); ++i) {
double a = Traits::direction(i) * node.element.fitness()[i];
double b = Traits::direction(i) * boxIt->element.fitness()[i];
d1 += pow( (a - node.discretized[i]) / eps[i], 2.0);
d2 += pow( (b - node.discretized[i]) / eps[i], 2.0);
double d1 = 0.0;
double d2 = 0.0;
for (unsigned i = 0; i < node.element.fitness().size(); ++i) {
double a = Traits::direction(i) * node.element.fitness()[i] * inv_eps[i];
double b = Traits::direction(i) * archive[box].element.fitness()[i] * inv_eps[i];
d1 += pow( a - node.discretized[i], 2.0);
d2 += pow( b - node.discretized[i], 2.0);
}
if (d1 > d2) {
archive[box] = node;
}
break;
}
if (d1 > d2) {
*boxIt = node;
}
break;
}
}
}
/*
static unsigned counter = 0;
if (++counter % 500 == 0) {
std::vector<long> mins(archive[0].discretized.size(), std::numeric_limits<long>::max());
std::vector<long> maxs(archive[0].discretized.size(), std::numeric_limits<long>::min());
for (unsigned i = 0; i < archive.size(); ++i) {
for (unsigned dim = 0; dim < archive[i].discretized.size(); ++dim) {
mins[dim] = std::min( mins[dim], archive[i].discretized[dim] );
maxs[dim] = std::max( maxs[dim], archive[i].discretized[dim] );
}
}
std::cout << "Range ";
for (unsigned dim = 0; dim < mins.size(); ++dim) {
std::cout << (maxs[dim] - mins[dim]) << ' ';
}
std::cout << archive.size() << std::endl;
}*/
if (archive.size() > max_size) {
unsigned idx = rng.random(archive.size());
if (idx != archive.size()-1) std::swap(archive[idx], archive.back());
archive.pop_back();
}
}
void appendTo(eoPop<EOT>& pop) const {
@ -111,9 +150,7 @@ exitLoop:
const EOT& selectRandom() const {
int i = rng.random(archive.size());
typename archive_t::const_iterator it = archive.begin();
while (i-- > 0) it++;
return it->element;
return archive[i].element;
}
};

View file

@ -4,20 +4,6 @@ using namespace std;
namespace detail {
namespace {
struct CompareOn {
unsigned dim;
double tol;
CompareOn(unsigned d, double t) : dim(d), tol(t) {}
bool operator()(const FitnessInfo& a, const FitnessInfo& b) {
return a.fitness[dim] > b.fitness[dim] && fabs(a.fitness[dim] - b.fitness[dim]) > tol;
}
};
} // end anonymous namespace
void one_objective(std::vector<FitnessInfo>& fitness, std::vector< std::vector<FitnessInfo> >& front, double tol)
{
std::sort(fitness.begin(), fitness.end(), CompareOn(0, tol));
@ -136,6 +122,7 @@ void m_objectives(std::vector<FitnessInfo>& fitness, std::vector< std::vector<Fi
}
void front_sorter_impl(std::vector<FitnessInfo>& fitness, std::vector< std::vector<FitnessInfo> >& front_indices, double tol) {
switch (fitness[0].fitness.size())
{
case 1:

View file

@ -43,6 +43,18 @@ namespace detail {
FitnessInfo(const std::vector<double>& fitness_, unsigned index_) : fitness(fitness_), index(index_) {}
};
struct CompareOn {
unsigned dim;
double tol;
CompareOn(unsigned d, double t) : dim(d), tol(t) {}
bool operator()(const FitnessInfo& a, const FitnessInfo& b) {
return a.fitness[dim] > b.fitness[dim] && fabs(a.fitness[dim] - b.fitness[dim]) > tol;
}
};
extern void front_sorter_impl(std::vector<FitnessInfo>& fitness, std::vector< std::vector<FitnessInfo> >& fronts, double tol);
} // namespace detail

View file

@ -105,12 +105,44 @@ namespace dominance {
// else b dominates a (check for non-dominance done above
return second;
}
template <class IntType>
inline dominance_result check_discrete(const std::vector<IntType>& a, const std::vector<IntType>& b) {
bool all_equal = true;
bool a_better_in_one = false;
bool b_better_in_one = false;
for (unsigned i = 0; i < a.size(); ++i) {
if ( a[i] != b[i] ) {
all_equal = false;
if (a[i] > b[i]) {
a_better_in_one = true;
} else {
b_better_in_one = true;
}
// check if we need to go on
if (a_better_in_one && b_better_in_one) return non_dominated;
}
}
if (all_equal) return non_dominated_equal;
if (a_better_in_one) return first;
// else b dominates a (check for non-dominance done above
return second;
}
template <class FitnessTraits>
inline dominance_result check(const std::vector<double>& p1, const std::vector<double>& p2) {
return check<FitnessTraits>(p1, p2, FitnessTraits::tol());
}
} // namespace dominance
@ -134,10 +166,22 @@ public :
eoMOFitness() : std::vector<double>(FitnessTraits::nObjectives()), worthValid(false) { reset(); }
eoMOFitness(double def) : std::vector<double>(FitnessTraits::nObjectives(),def), worthValid(false) {}
explicit eoMOFitness(double def) : std::vector<double>(FitnessTraits::nObjectives(),def), worthValid(false) {}
// Ctr from a std::vector<double>
eoMOFitness(std::vector<double> & _v) : std::vector<double>(_v), worthValid(false) {}
explicit eoMOFitness(std::vector<double> & _v) : std::vector<double>(_v), worthValid(false) {}
virtual ~eoMOFitness() {} // in case someone wants to subclass
eoMOFitness(const eoMOFitness<FitnessTraits>& org) { operator=(org); }
eoMOFitness<FitnessTraits>& operator=(const eoMOFitness<FitnessTraits>& org) {
std::vector<double>::operator=(org);
worth = org.worth;
worthValid = org.worthValid;
return *this;
}
void reset() {
@ -177,7 +221,9 @@ public :
std::vector<double> f;
for (unsigned j = 0; j < FitnessTraits::nObjectives(); ++j) {
if (FitnessTraits::direction(j) != 0) f.push_back( FitnessTraits::direction(j) * this->operator[](j));
if (FitnessTraits::direction(j) != 0) {
f.push_back( FitnessTraits::direction(j) * this->operator[](j));
}
}
return f;
@ -192,17 +238,17 @@ public :
/// compare *not* on dominance, but on worth
bool operator<(const eoMOFitness<FitnessTraits>& _other) const
{
return getWorth() > _other.getWorth();
return getWorth() < _other.getWorth();
}
bool operator>(const eoMOFitness<FitnessTraits>& _other) const
{
return _other < *this;
return _other > *this;
}
bool operator<=(const eoMOFitness<FitnessTraits>& _other) const
{
return getWorth() >= _other.getWorth();
return getWorth() <= _other.getWorth();
}
bool operator>=(const eoMOFitness<FitnessTraits>& _other) const
@ -242,4 +288,11 @@ public :
};
namespace dominance {
template <class FitnessTraits>
inline dominance_result check(const eoMOFitness<FitnessTraits>& p1, const eoMOFitness<FitnessTraits>& p2) {
return check<FitnessTraits>(p1, p2, FitnessTraits::tol());
}
}
#endif

View file

@ -1,6 +1,8 @@
#include <eoNSGA_II_Eval.h>
namespace nsga2 {
using namespace std;
typedef std::pair<double, unsigned> double_index_pair;
@ -15,8 +17,6 @@ namespace nsga2 {
void assign_worths(const std::vector<detail::FitnessInfo>& front, unsigned rank, std::vector<double>& worths) {
unsigned i;
unsigned nObjectives = front[0].fitness.size(); //traits::nObjectives(); //_pop[_cf[0]].fitness().size();
std::vector<double> niche_distance(front.size());
@ -25,7 +25,7 @@ namespace nsga2 {
{
std::vector<std::pair<double, unsigned> > performance(front.size());
for (i =0; i < front.size(); ++i)
for (unsigned i =0; i < front.size(); ++i)
{
performance[i].first = front[i].fitness[o];
performance[i].second = i;
@ -35,16 +35,18 @@ namespace nsga2 {
std::vector<double> nc(front.size(), 0.0);
for (i = 1; i < front.size()-1; ++i)
for (unsigned i = 1; i < front.size()-1; ++i)
{ // and yet another level of indirection
nc[performance[i].second] = performance[i+1].first - performance[i-1].first;
}
// set boundary at max_dist + 1 (so it will get chosen over all the others
//nc[performance[0].second] += 0;
nc[performance.back().second] += std::numeric_limits<double>::infinity(); // best on objective
nc[performance.back().second] = std::numeric_limits<double>::infinity(); // best on objective
for (i = 0; i < nc.size(); ++i)
for (unsigned i = 0; i < nc.size(); ++i)
{
niche_distance[i] += nc[i];
}
@ -53,12 +55,12 @@ namespace nsga2 {
// now we've got niche_distances, scale them between (0, 1), making sure that infinities get maximum rank
double max = 0;
for (unsigned i = 0; i < niche_distance[i]; ++i) {
for (unsigned i = 0; i < niche_distance.size(); ++i) {
if (niche_distance[i] != std::numeric_limits<double>::infinity()) {
max = std::max(max, niche_distance[i]);
}
}
for (unsigned i = 0; i < front.size(); ++i) {
double dist = niche_distance[i];
if (dist == std::numeric_limits<double>::infinity()) {

View file

@ -34,15 +34,15 @@ class eoNSGA_II_Eval : public eoMOEval<EOT>
typename eoFrontSorter<EOT>::front_t front = sorter(pop);
// calculate worths
std::vector<double> worths(pop.size());
for (unsigned i = 0; i < front.size(); ++i) {
nsga2::assign_worths(front[i], front.size() - i, worths);
}
// set worths
for (unsigned i = 0; i < pop.size(); ++i) {
typename EOT::Fitness f = pop[i]->fitness();
f.setWorth(worths[i]);
pop[i]->fitness(f);
pop[i]->fitnessReference().setWorth( worths[i]);
}
}

View file

@ -7,109 +7,221 @@ double calc_distance(const std::vector<double>& f1, const std::vector<double>& f
double dist = 0;
for (unsigned i = 0; i < f1.size(); ++i) {
double d = (f1[i] - f2[i]);
dist += d*d;
dist += pow(fabs(d), 1./2);
//dist += pow(fabs(d), 2.0);
}
return dist;
}
unsigned assign_worths(std::vector<detail::FitnessInfo> front, unsigned rank, std::vector<double>& worths) {
unsigned nDim = front[0].fitness.size();
unsigned nDim = front[0].fitness.size();
// find boundary points
std::vector<unsigned> processed(nDim);
// find boundary points
std::vector<unsigned> processed(nDim);
for (unsigned i = 1; i < front.size(); ++i) {
for (unsigned dim = 0; dim < nDim; ++dim) {
if (front[i].fitness[dim] > front[processed[dim]].fitness[dim]) {
processed[dim] = i;
for (unsigned i = 1; i < front.size(); ++i) {
for (unsigned dim = 0; dim < nDim; ++dim) {
if (front[i].fitness[dim] > front[processed[dim]].fitness[dim]) {
processed[dim] = i;
}
}
}
}
// assign fitness to processed and store in boundaries
std::vector<detail::FitnessInfo> boundaries;
for (unsigned i = 0; i < processed.size(); ++i) {
worths[ front[ processed[i] ].index] = rank;
//cout << "Boundary " << i << ' ' << front[processed[i]].index << ' ' << parents[ front[ processed[i] ].index]->fitness() << endl;
boundaries.push_back( front[ processed[i] ] );
}
rank--;
// clean up processed (make unique)
sort(processed.begin(), processed.end()); // swap out last first
for (unsigned i = 1; i < processed.size(); ++i) {
if (processed[i] == processed[i-1]) {
std::swap(processed[i], processed.back());
processed.pop_back();
--i;
}
}
// remove boundaries from front
while (processed.size()) {
unsigned idx = processed.back();
//std::cout << idx << ' ' ;
processed.pop_back();
std::swap(front.back(), front[idx]);
front.pop_back();
}
//std::cout << std::endl;
// calculate distances
std::vector<double> distances(front.size(), std::numeric_limits<double>::infinity());
unsigned selected = 0;
// select based on maximum distance to nearest processed point
for (unsigned i = 0; i < front.size(); ++i) {
for (unsigned k = 0; k < boundaries.size(); ++k) {
double d = calc_distance( front[i].fitness, boundaries[k].fitness);
if (d < distances[i]) {
distances[i] = d;
}
}
if (distances[i] > distances[selected]) {
selected = i;
}
}
while (!front.empty()) {
detail::FitnessInfo last = front[selected];
std::swap(front[selected], front.back());
front.pop_back();
std::swap(distances[selected], distances.back());
distances.pop_back();
// set worth
worths[last.index] = rank--;
if (front.empty()) break;
selected = 0;
for (unsigned i = 0; i < front.size(); ++i) {
double d = calc_distance(front[i].fitness, last.fitness);
// assign fitness to processed and store in boundaries
std::vector<detail::FitnessInfo> boundaries;
for (unsigned i = 0; i < processed.size(); ++i) {
if (d < distances[i]) {
distances[i] = d;
}
worths[ front[ processed[i] ].index] = rank;
//cout << "Boundary " << i << ' ' << front[processed[i]].index << ' ' << parents[ front[ processed[i] ].index]->fitness() << endl;
if (distances[i] >= distances[selected]) {
boundaries.push_back( front[ processed[i] ] );
}
rank--;
// clean up processed (make unique)
sort(processed.begin(), processed.end()); // swap out last first
for (unsigned i = 1; i < processed.size(); ++i) {
if (processed[i] == processed[i-1]) {
std::swap(processed[i], processed.back());
processed.pop_back();
--i;
}
}
// remove boundaries from front
while (processed.size()) {
unsigned idx = processed.back();
//std::cout << idx << ' ' ;
processed.pop_back();
std::swap(front.back(), front[idx]);
front.pop_back();
}
//std::cout << std::endl;
// calculate distances
std::vector<double> distances(front.size(), std::numeric_limits<double>::infinity());
unsigned selected = 0;
// select based on maximum distance to nearest processed point
for (unsigned i = 0; i < front.size(); ++i) {
for (unsigned k = 0; k < boundaries.size(); ++k) {
double d = calc_distance( front[i].fitness, boundaries[k].fitness);
if (d < distances[i]) {
distances[i] = d;
}
}
if (distances[i] > distances[selected]) {
selected = i;
}
}
while (!front.empty()) {
detail::FitnessInfo last = front[selected];
std::swap(front[selected], front.back());
front.pop_back();
std::swap(distances[selected], distances.back());
distances.pop_back();
// set worth
worths[last.index] = rank--;
if (front.empty()) break;
selected = 0;
for (unsigned i = 0; i < front.size(); ++i) {
double d = calc_distance(front[i].fitness, last.fitness);
if (d < distances[i]) {
distances[i] = d;
}
if (distances[i] >= distances[selected]) {
selected = i;
}
}
}
return rank;
}
return rank;
unsigned assign_worths2(std::vector<detail::FitnessInfo> front, unsigned rank, std::vector<double>& worths) {
unsigned nDim = front[0].fitness.size();
// find boundary points
std::vector<unsigned> processed(nDim);
for (unsigned i = 1; i < front.size(); ++i) {
for (unsigned dim = 0; dim < nDim; ++dim) {
if (front[i].fitness[dim] > front[processed[dim]].fitness[dim]) {
processed[dim] = i;
}
}
}
// assign fitness to processed and store in boundaries
std::vector<detail::FitnessInfo> boundaries;
for (unsigned i = 0; i < processed.size(); ++i) {
worths[ front[ processed[i] ].index] = rank;
//cout << "Boundary " << i << ' ' << front[processed[i]].index << ' ' << parents[ front[ processed[i] ].index]->fitness() << endl;
boundaries.push_back( front[ processed[i] ] );
}
rank--;
// clean up processed (make unique)
sort(processed.begin(), processed.end()); // swap out last first
for (unsigned i = 1; i < processed.size(); ++i) {
if (processed[i] == processed[i-1]) {
std::swap(processed[i], processed.back());
processed.pop_back();
--i;
}
}
// remove boundaries from front
while (processed.size()) {
unsigned idx = processed.back();
//std::cout << idx << ' ' ;
processed.pop_back();
std::swap(front.back(), front[idx]);
front.pop_back();
}
//std::cout << std::endl;
// calculate distances
std::vector< std::vector<double> > distances(front.size(), std::vector<double>(nDim, std::numeric_limits<double>::infinity()));
double bestsum = 0;
unsigned selected = 0;
// select based on maximum distance to nearest processed point
for (unsigned i = 0; i < front.size(); ++i) {
for (unsigned k = 0; k < boundaries.size(); ++k) {
for (unsigned dim = 0; dim < nDim; ++dim) {
double d = front[i].fitness[dim] - boundaries[k].fitness[dim];
if (d > 0 && d < distances[i][dim]) {
distances[i][dim] = d;
}
}
}
double sum = 0;
for (unsigned dim = 0; dim < nDim; ++dim) sum += distances[i][dim];
if (sum > bestsum) {
selected = i;
bestsum = sum;
}
}
while (!front.empty()) {
detail::FitnessInfo last = front[selected];
std::swap(front[selected], front.back());
front.pop_back();
std::swap(distances[selected], distances.back());
distances.pop_back();
// set worth
worths[last.index] = rank--;
if (front.empty()) break;
selected = 0;
bestsum = 0;
for (unsigned i = 0; i < front.size(); ++i) {
double sum = 0;
for (unsigned dim = 0; dim < nDim; ++dim) {
double d = front[i].fitness[dim] - last.fitness[dim];
if (d > 0 && d < distances[i][dim]) {
distances[i][dim] = d;
}
sum += distances[i][dim];
}
if (sum > bestsum) {
selected = i;
bestsum = sum;
}
}
}
return rank;
}
}