Commit 42e6658e authored by snuverink_j's avatar snuverink_j
Browse files

fix 236: add option to add reference point for hypervolume calculation

parent dec39823
...@@ -61,11 +61,13 @@ public: ...@@ -61,11 +61,13 @@ public:
* @param[in] dim number of objectives * @param[in] dim number of objectives
* @param[in] comms available to the optimizer * @param[in] comms available to the optimizer
* @param[in] args the user passed on the command line * @param[in] args the user passed on the command line
* @param[in] hypervolRef hypervolume reference point
*/ */
FixedPisaNsga2(Expressions::Named_t objectives, FixedPisaNsga2(Expressions::Named_t objectives,
Expressions::Named_t constraints, Expressions::Named_t constraints,
DVarContainer_t dvars, size_t dim, Comm::Bundle_t comms, DVarContainer_t dvars, size_t dim, Comm::Bundle_t comms,
CmdArguments_t args); CmdArguments_t args,
std::vector<double> hypervolRef);
~FixedPisaNsga2(); ~FixedPisaNsga2();
...@@ -175,7 +177,7 @@ private: ...@@ -175,7 +177,7 @@ private:
DVarContainer_t dvars_m; DVarContainer_t dvars_m;
/// command line arguments specified by the user /// command line arguments specified by the user
CmdArguments_t args_; CmdArguments_t args_m;
/// size of initial population /// size of initial population
size_t alpha_m; size_t alpha_m;
...@@ -205,6 +207,9 @@ private: ...@@ -205,6 +207,9 @@ private:
double conv_hvol_progress_; double conv_hvol_progress_;
double hvol_progress_; double hvol_progress_;
/// hypervolume reference point
std::vector<double> hvol_ref_m;
/// file header for result files contains this parameter description /// file header for result files contains this parameter description
std::string file_param_descr_; std::string file_param_descr_;
......
...@@ -32,15 +32,17 @@ FixedPisaNsga2<CO, MO>::FixedPisaNsga2( ...@@ -32,15 +32,17 @@ FixedPisaNsga2<CO, MO>::FixedPisaNsga2(
Expressions::Named_t constraints, Expressions::Named_t constraints,
DVarContainer_t dvars, DVarContainer_t dvars,
size_t dim, Comm::Bundle_t comms, size_t dim, Comm::Bundle_t comms,
CmdArguments_t args) CmdArguments_t args,
std::vector<double> hypervolRef)
: Optimizer(comms.opt) : Optimizer(comms.opt)
, statistics_(new Statistics<size_t>("individuals")) , statistics_(new Statistics<size_t>("individuals"))
, comms_(comms) , comms_(comms)
, objectives_m(objectives) , objectives_m(objectives)
, constraints_m(constraints) , constraints_m(constraints)
, dvars_m(dvars) , dvars_m(dvars)
, args_(args) , args_m(args)
, dim_m(dim) , dim_m(dim)
, hvol_ref_m(hypervolRef)
{ {
my_local_pid_ = 0; my_local_pid_ = 0;
MPI_Comm_rank(comms_.opt, &my_local_pid_); MPI_Comm_rank(comms_.opt, &my_local_pid_);
...@@ -150,7 +152,7 @@ void FixedPisaNsga2<CO, MO>::initialize() { ...@@ -150,7 +152,7 @@ void FixedPisaNsga2<CO, MO>::initialize() {
bool compHyvol = (objectives_m.size() > (hyper_opt / 2 + 1)); bool compHyvol = (objectives_m.size() > (hyper_opt / 2 + 1));
if (compHyvol) if (compHyvol)
current_hvol_ = current_hvol_ =
variator_m->population()->computeHypervolume(comms_.island_id); variator_m->population()->computeHypervolume(comms_.island_id, hvol_ref_m);
boost::chrono::duration<double> total = boost::chrono::duration<double> total =
boost::chrono::system_clock::now() - run_clock_start_; boost::chrono::system_clock::now() - run_clock_start_;
...@@ -262,9 +264,8 @@ bool FixedPisaNsga2<CO, MO>::onMessage(MPI_Status status, size_t length) { ...@@ -262,9 +264,8 @@ bool FixedPisaNsga2<CO, MO>::onMessage(MPI_Status status, size_t length) {
ind->objectives.clear(); ind->objectives.clear();
//XXX: check order of genes //XXX: check order of genes
reqVarContainer_t::iterator itr; reqVarContainer_t::iterator itr = res.begin();
std::map<std::string, double> vars; for(; itr != res.end(); ++itr) {
for(itr = res.begin(); itr != res.end(); itr++) {
// mark invalid if expression could not be evaluated or constraint does not hold // mark invalid if expression could not be evaluated or constraint does not hold
if(!itr->second.is_valid || (itr->second.value.size() > 1 && !itr->second.value[0])) { if(!itr->second.is_valid || (itr->second.value.size() > 1 && !itr->second.value[0])) {
std::ostringstream dump; std::ostringstream dump;
...@@ -321,7 +322,7 @@ void FixedPisaNsga2<CO, MO>::postPoll() { ...@@ -321,7 +322,7 @@ void FixedPisaNsga2<CO, MO>::postPoll() {
bool compHyvol = (objectives_m.size() > (hyper_opt / 2 + 1)); bool compHyvol = (objectives_m.size() > (hyper_opt / 2 + 1));
if (compHyvol) { if (compHyvol) {
double hvol = double hvol =
variator_m->population()->computeHypervolume(comms_.island_id); variator_m->population()->computeHypervolume(comms_.island_id, hvol_ref_m);
hvol_progress_ = fabs(current_hvol_ - hvol) / current_hvol_; hvol_progress_ = fabs(current_hvol_ - hvol) / current_hvol_;
current_hvol_ = hvol; current_hvol_ = hvol;
} }
...@@ -370,7 +371,7 @@ void FixedPisaNsga2<CO, MO>::exchangeSolutionStates() { ...@@ -370,7 +371,7 @@ void FixedPisaNsga2<CO, MO>::exchangeSolutionStates() {
typedef typename FixedPisaNsga2::Individual_t individual; typedef typename FixedPisaNsga2::Individual_t individual;
size_t num_masters = args_->getArg<size_t>("num-masters", 1, false); size_t num_masters = args_m->getArg<size_t>("num-masters", 1, false);
if(num_masters <= 1 || if(num_masters <= 1 ||
exchangeSolStateFreq_m == 0 || exchangeSolStateFreq_m == 0 ||
...@@ -660,7 +661,7 @@ template< template <class> class CO, template <class> class MO > ...@@ -660,7 +661,7 @@ template< template <class> class CO, template <class> class MO >
void FixedPisaNsga2<CO, MO>::dumpPopulationToFile() { void FixedPisaNsga2<CO, MO>::dumpPopulationToFile() {
// only dump old data format if the user requests it // only dump old data format if the user requests it
if(! args_->getArg<bool>("dump-dat", false, false)) return; if(! args_m->getArg<bool>("dump-dat", false, false)) return;
typedef typename FixedPisaNsga2::Individual_t individual; typedef typename FixedPisaNsga2::Individual_t individual;
boost::shared_ptr<individual> temp; boost::shared_ptr<individual> temp;
......
...@@ -162,7 +162,9 @@ public: ...@@ -162,7 +162,9 @@ public:
} }
double computeHypervolume(size_t island_id) { double computeHypervolume(size_t island_id, const std::vector<double>& referencePoint) {
// protection check
if (individuals.empty() == true) return -1;
std::ofstream file; std::ofstream file;
std::ostringstream filename; std::ostringstream filename;
...@@ -186,7 +188,7 @@ public: ...@@ -186,7 +188,7 @@ public:
file.flush(); file.flush();
file.close(); file.close();
hypervolume_ = Hypervolume::FromFile(filename.str()); hypervolume_ = Hypervolume::FromFile(filename.str(), referencePoint);
return hypervolume_; return hypervolume_;
} }
......
...@@ -113,13 +113,15 @@ public: ...@@ -113,13 +113,15 @@ public:
functionDictionary_t known_expr_funcs, functionDictionary_t known_expr_funcs,
const DVarContainer_t &dvar, const DVarContainer_t &dvar,
const Expressions::Named_t &obj, const Expressions::Named_t &obj,
const Expressions::Named_t &cons) const Expressions::Named_t &cons,
std::vector<double> hypervolRef = {})
: Poller(comm->mpiComm()) : Poller(comm->mpiComm())
, comm_(comm) , comm_(comm)
, cmd_args_(args) , cmd_args_(args)
, objectives_(obj) , objectives_(obj)
, constraints_(cons) , constraints_(cons)
, dvars_(dvar) , dvars_(dvar)
, hypervolRef_(hypervolRef)
{ {
setup(known_expr_funcs); setup(known_expr_funcs);
} }
...@@ -157,9 +159,11 @@ protected: ...@@ -157,9 +159,11 @@ protected:
bool has_opt_converged_; bool has_opt_converged_;
bool continue_polling_; bool continue_polling_;
Expressions::Named_t objectives_; Expressions::Named_t objectives_; ///< objectives
Expressions::Named_t constraints_; Expressions::Named_t constraints_; ///< constraints
DVarContainer_t dvars_; DVarContainer_t dvars_; ///< design variables
std::vector<double> hypervolRef_; ///< hypervolume reference point
// keep track of state of all workers // keep track of state of all workers
std::vector<bool> is_worker_idle_; std::vector<bool> is_worker_idle_;
...@@ -248,7 +252,7 @@ protected: ...@@ -248,7 +252,7 @@ protected:
boost::scoped_ptr<Opt_t> opt( boost::scoped_ptr<Opt_t> opt(
new Opt_t(objectives_, constraints_, dvars_, objectives_.size(), new Opt_t(objectives_, constraints_, dvars_, objectives_.size(),
comm_->getBundle(), cmd_args_)); comm_->getBundle(), cmd_args_, hypervolRef_));
opt->initialize(); opt->initialize();
std::cout << "Stop Opt.." << std::endl; std::cout << "Stop Opt.." << std::endl;
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
/** /**
* \class CmdArguments * \class CmdArguments
* \brief Parsing commando line arguments * \brief Parsing command line arguments
* *
* In order to have a flexible framework, each component implementation gets * In order to have a flexible framework, each component implementation gets
* access to all command line arguments. * access to all command line arguments.
......
...@@ -35,80 +35,72 @@ ...@@ -35,80 +35,72 @@
#include <sys/time.h> #include <sys/time.h>
#include <sys/resource.h> #include <sys/resource.h>
#include "wfg.h" #include "wfg.h"
#include "avl.h" #include "avl.h"
#include <string> #include <string>
#include "hypervolume.h" #include "hypervolume.h"
#define MAXIMISING true #define MAXIMISING true
#if MAXIMISING #if MAXIMISING
#define BEATS(x,y) (x > y) #define BEATS(x,y) (x > y)
//#define BEATSEQ(x,y) (x >= y)
#else #else
#define BEATS(x,y) (x < y) #define BEATS(x,y) (x < y)
//#define BEATSEQ(x,y) (x <= y)
#endif #endif
#define WORSE(x,y) (BEATS(y,x) ? (x) : (y)) #define WORSE(x,y) (BEATS(y,x) ? (x) : (y))
//#define BETTER(x,y) (BEATS(y,x) ? (y) : (x))
namespace Hypervolume {
int n; // the number of objectives
POINT ref; // the reference point
int n; // the number of objectives FRONT *fs; // memory management stuff
POINT ref; // the reference point int fr = 0; // current depth
int frmax = -1; // max depth malloced so far (for hyper_opt = 0)
int maxm = 0; // identify the biggest fronts in the file
int maxn = 0;
FRONT *fs; // memory management stuff static avl_tree_t *tree;
int fr = 0; // current depth double hv(FRONT);
int frmax = -1; // max depth malloced so far (for hyper_opt = 0)
int maxm = 0; // identify the biggest fronts in the file
int maxn = 0;
static int compare_tree_asc( const void *p1, const void *p2)
static avl_tree_t *tree; {
double hv(FRONT);
using namespace Hypervolume;
static int compare_tree_asc( const void *p1, const void *p2)
{
const double x1= *((const double *)p1+1); const double x1= *((const double *)p1+1);
const double x2= *((const double *)p2+1); const double x2= *((const double *)p2+1);
if (x1 != x2) return (x1 > x2) ? -1 : 1; if (x1 != x2) return (x1 > x2) ? -1 : 1;
else return 0; else return 0;
} }
int greater(const void *v1, const void *v2) int greater(const void *v1, const void *v2)
// this sorts points improving in the last objective // this sorts points improving in the last objective
{ {
POINT p = *(POINT*)v1; POINT p = *(POINT*)v1;
POINT q = *(POINT*)v2; POINT q = *(POINT*)v2;
#if hyper_opt == 1 #if hyper_opt == 1
for (int i = n - fr - 1; i >= 0; i--) for (int i = n - fr - 1; i >= 0; i--)
#else #else
for (int i = n - 1; i >= 0; i--) for (int i = n - 1; i >= 0; i--)
#endif #endif
if BEATS(p.objectives[i],q.objectives[i]) return 1; if BEATS(p.objectives[i],q.objectives[i]) return 1;
else else
if BEATS(q.objectives[i],p.objectives[i]) return -1; if BEATS(q.objectives[i],p.objectives[i]) return -1;
return 0; return 0;
} }
int dominates2way(POINT p, POINT q) int dominates2way(POINT p, POINT q)
// returns -1 if p dominates q, 1 if q dominates p, 2 if p == q, 0 otherwise // returns -1 if p dominates q, 1 if q dominates p, 2 if p == q, 0 otherwise
{ {
// domination could be checked in either order // domination could be checked in either order
#if hyper_opt == 1 #if hyper_opt == 1
for (int i = n - fr - 1; i >= 0; i--) for (int i = n - fr - 1; i >= 0; i--)
#else #else
for (int i = n - 1; i >= 0; i--) for (int i = n - 1; i >= 0; i--)
#endif #endif
if BEATS(p.objectives[i],q.objectives[i]) if BEATS(p.objectives[i],q.objectives[i])
{for (int j = i - 1; j >= 0; j--) {for (int j = i - 1; j >= 0; j--)
if BEATS(q.objectives[j],p.objectives[j]) return 0; if BEATS(q.objectives[j],p.objectives[j]) return 0;
...@@ -119,14 +111,14 @@ int dominates2way(POINT p, POINT q) ...@@ -119,14 +111,14 @@ int dominates2way(POINT p, POINT q)
if BEATS(p.objectives[j],q.objectives[j]) return 0; if BEATS(p.objectives[j],q.objectives[j]) return 0;
return 1;} return 1;}
return 2; return 2;
} }
void makeDominatedBit(FRONT ps, int p) void makeDominatedBit(FRONT ps, int p)
// creates the front ps[p+1 ..] in fs[fr], with each point bounded by ps[p] and dominated points removed // creates the front ps[p+1 ..] in fs[fr], with each point bounded by ps[p] and dominated points removed
{ {
// when hyper_opt = 0 each new frame is allocated as needed, because the worst-case needs #frames = #points // when hyper_opt = 0 each new frame is allocated as needed, because the worst-case needs #frames = #points
#if hyper_opt == 0 #if hyper_opt == 0
if (fr > frmax) if (fr > frmax)
{frmax = fr; {frmax = fr;
fs[fr].points = (POINT*) malloc( maxm * sizeof(POINT)); fs[fr].points = (POINT*) malloc( maxm * sizeof(POINT));
...@@ -135,7 +127,7 @@ void makeDominatedBit(FRONT ps, int p) ...@@ -135,7 +127,7 @@ void makeDominatedBit(FRONT ps, int p)
fs[fr].points[j].objectives = (OBJECTIVE*) malloc(maxn * sizeof(OBJECTIVE)); fs[fr].points[j].objectives = (OBJECTIVE*) malloc(maxn * sizeof(OBJECTIVE));
} }
} }
#endif #endif
int z = ps.nPoints - 1 - p; int z = ps.nPoints - 1 - p;
for (int i = 0; i < z; i++) for (int i = 0; i < z; i++)
...@@ -163,39 +155,39 @@ void makeDominatedBit(FRONT ps, int p) ...@@ -163,39 +155,39 @@ void makeDominatedBit(FRONT ps, int p)
fs[fr].nPoints++;} fs[fr].nPoints++;}
} }
fr++; fr++;
} }
double hv2(FRONT ps) double hv2(FRONT ps)
// returns the hypervolume of ps[0 ..] in 2D // returns the hypervolume of ps[0 ..] in 2D
// assumes that ps is sorted improving // assumes that ps is sorted improving
{ {
double volume = fabs((ps.points[0].objectives[0] - ref.objectives[0]) * double volume = fabs((ps.points[0].objectives[0] - ref.objectives[0]) *
(ps.points[0].objectives[1] - ref.objectives[1])); (ps.points[0].objectives[1] - ref.objectives[1]));
for (int i = 1; i < ps.nPoints; i++) for (int i = 1; i < ps.nPoints; i++)
volume += fabs((ps.points[i].objectives[0] - ref.objectives[0]) * volume += fabs((ps.points[i].objectives[0] - ref.objectives[0]) *
(ps.points[i].objectives[1] - ps.points[i - 1].objectives[1])); (ps.points[i].objectives[1] - ps.points[i - 1].objectives[1]));
return volume; return volume;
} }
double hv3_AVL(FRONT ps) double hv3_AVL(FRONT ps)
/* hv3_AVL: 3D algorithm code taken from version hv-1.2 available at /* hv3_AVL: 3D algorithm code taken from version hv-1.2 available at
http://iridia.ulb.ac.be/~manuel/hypervolume and proposed by: http://iridia.ulb.ac.be/~manuel/hypervolume and proposed by:
Carlos M. Fonseca, Lus Paquete, and Manuel Lpez-Ibez. An improved Carlos M. Fonseca, Lus Paquete, and Manuel Lpez-Ibez. An improved
dimension-sweep algorithm for the hypervolume indicator. In IEEE dimension-sweep algorithm for the hypervolume indicator. In IEEE
Congress on Evolutionary Computation, pages 1157-1163, Vancouver, Congress on Evolutionary Computation, pages 1157-1163, Vancouver,
Canada, July 2006. Canada, July 2006.
Copyright (c) 2009 Copyright (c) 2009
Carlos M. Fonseca <cmfonsec@ualg.pt> Carlos M. Fonseca <cmfonsec@ualg.pt>
Manuel Lopez-Ibanez <manuel.lopez-ibanez@ulb.ac.be> Manuel Lopez-Ibanez <manuel.lopez-ibanez@ulb.ac.be>
Luis Paquete <paquete@dei.uc.pt> Luis Paquete <paquete@dei.uc.pt>
*/ */
// returns the hypervolume of ps[0 ..] in 3D // returns the hypervolume of ps[0 ..] in 3D
// assumes that ps is sorted improving // assumes that ps is sorted improving
{ {
avl_init_node(ps.points[ps.nPoints-1].tnode,ps.points[ps.nPoints-1].objectives); avl_init_node(ps.points[ps.nPoints-1].tnode,ps.points[ps.nPoints-1].objectives);
avl_insert_top(tree,ps.points[ps.nPoints-1].tnode); avl_insert_top(tree,ps.points[ps.nPoints-1].tnode);
...@@ -284,22 +276,22 @@ Canada, July 2006. ...@@ -284,22 +276,22 @@ Canada, July 2006.
} }
avl_clear_tree(tree); avl_clear_tree(tree);
return hyperv; return hyperv;
} }
double inclhv(POINT p) double inclhv(POINT p)
// returns the inclusive hypervolume of p // returns the inclusive hypervolume of p
{ {
double volume = 1; double volume = 1;
for (int i = 0; i < n; i++) for (int i = 0; i < n; i++)
volume *= fabs(p.objectives[i] - ref.objectives[i]); volume *= fabs(p.objectives[i] - ref.objectives[i]);
return volume; return volume;
} }
double exclhv(FRONT ps, int p) double exclhv(FRONT ps, int p)
// returns the exclusive hypervolume of ps[p] relative to ps[p+1 ..] // returns the exclusive hypervolume of ps[p] relative to ps[p+1 ..]
{ {
double volume = inclhv(ps.points[p]); double volume = inclhv(ps.points[p]);
if (ps.nPoints > p + 1) if (ps.nPoints > p + 1)
{ {
...@@ -308,27 +300,27 @@ double exclhv(FRONT ps, int p) ...@@ -308,27 +300,27 @@ double exclhv(FRONT ps, int p)
fr--; fr--;
} }
return volume; return volume;
} }
double hv(FRONT ps) double hv(FRONT ps)
// returns the hypervolume of ps[0 ..] // returns the hypervolume of ps[0 ..]
{ {
#if hyper_opt > 0 #if hyper_opt > 0
qsort(ps.points, ps.nPoints, sizeof(POINT), greater); qsort(ps.points, ps.nPoints, sizeof(POINT), greater);
#endif #endif
#if hyper_opt == 2 #if hyper_opt == 2
if (n == 2) return hv2(ps); if (n == 2) return hv2(ps);
#elif hyper_opt == 3 #elif hyper_opt == 3
if (n == 3) return hv3_AVL(ps); if (n == 3) return hv3_AVL(ps);
#endif #endif
double volume = 0; double volume = 0;
#if hyper_opt <= 1 #if hyper_opt <= 1
for (int i = 0; i < ps.nPoints; i++) volume += exclhv(ps, i); for (int i = 0; i < ps.nPoints; i++) volume += exclhv(ps, i);
#else #else
n--; n--;
for (int i = ps.nPoints - 1; i >= 0; i--) for (int i = ps.nPoints - 1; i >= 0; i--)
// we can ditch dominated points here, // we can ditch dominated points here,
...@@ -336,26 +328,27 @@ double hv(FRONT ps) ...@@ -336,26 +328,27 @@ double hv(FRONT ps)
volume += fabs(ps.points[i].objectives[n] - ref.objectives[n]) * exclhv(ps, i); volume += fabs(ps.points[i].objectives[n] - ref.objectives[n]) * exclhv(ps, i);
n++; n++;
#endif #endif
return volume; return volume;
} }