PaGMO  1.1.5
wfg.cpp
1 /*****************************************************************************
2  * Copyright (C) 2004-2015 The PaGMO development team, *
3  * Advanced Concepts Team (ACT), European Space Agency (ESA) *
4  * *
5  * https://github.com/esa/pagmo *
6  * *
7  * act@esa.int *
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  * This program is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17  * GNU General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this program; if not, write to the *
21  * Free Software Foundation, Inc., *
22  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
23  *****************************************************************************/
24 
25 
26 #include "hv2d.h"
27 #include "wfg.h"
28 #include "base.h"
29 #include <algorithm>
30 #include <boost/bind.hpp>
31 
32 namespace pagmo { namespace util { namespace hv_algorithm {
33 
35 
38 bool wfg::cmp_points(double* a, double* b) const
39 {
40  for(int i = m_current_slice - 1; i >= 0 ; --i){
41  if (a[i] > b[i]) {
42  return true;
43  } else if(a[i] < b[i]) {
44  return false;
45  }
46  }
47  return false;
48 }
49 
51 wfg::wfg(const unsigned int stop_dimension) : m_current_slice(0), m_stop_dimension(stop_dimension)
52 {
53  if (stop_dimension < 2 ) {
54  pagmo_throw(value_error, "Stop dimension for WFG must be greater than or equal to 2");
55  }
56 }
57 
59 
67 double wfg::compute(std::vector<fitness_vector> &points, const fitness_vector &r_point) const
68 {
69  allocate_wfg_members(points, r_point);
70  double hv = compute_hv(1);
71  free_wfg_members();
72  return hv;
73 }
74 
76 
88 std::vector<double> wfg::contributions(std::vector<fitness_vector> &points, const fitness_vector &r_point) const
89 {
90  std::vector<double> c;
91  c.reserve(points.size());
92 
93  // Allocate the same members as for 'compute' method
94  allocate_wfg_members(points, r_point);
95 
96  // Prepare the memory for first front
97  double** fr = new double*[m_max_points];
98  for(unsigned int i = 0 ; i < m_max_points ; ++i) {
99  fr[i]=new double[m_current_slice];
100  }
101  m_frames[m_n_frames] = fr;
102  m_frames_size[m_n_frames] = 0;
103  ++m_n_frames;
104 
105  for(unsigned int p_idx = 0 ; p_idx < m_max_points ; ++p_idx) {
106  limitset(0, p_idx, 1);
107  c.push_back(exclusive_hv(p_idx, 1));
108  }
109 
110  // Free the contributions and the remaining WFG members
111  free_wfg_members();
112 
113  return c;
114 }
115 
117 void wfg::allocate_wfg_members(std::vector<fitness_vector> &points, const fitness_vector &r_point) const
118 {
119  m_max_points = points.size();
120  m_max_dim = r_point.size();
121 
122  m_refpoint = new double[m_max_dim];
123  for(unsigned int d_idx = 0 ; d_idx < m_max_dim ; ++d_idx) {
124  m_refpoint[d_idx] = r_point[d_idx];
125  }
126 
127  // Reserve the space beforehand for each level or recursion.
128  // WFG with slicing feature will not go recursively deeper than the dimension size.
129  m_frames = new double**[m_max_dim];
130  m_frames_size = new unsigned int[m_max_dim];
131 
132  // Copy the initial set into the frame at index 0.
133  double** fr = new double*[m_max_points];
134  for(unsigned int p_idx = 0 ; p_idx < m_max_points ; ++p_idx) {
135  fr[p_idx] = new double[m_max_dim];
136  for(unsigned int d_idx = 0 ; d_idx < m_max_dim ; ++d_idx) {
137  fr[p_idx][d_idx] = points[p_idx][d_idx];
138  }
139  }
140  m_frames[0] = fr;
141  m_frames_size[0] = m_max_points;
142  m_n_frames = 1;
143 
144  // Variable holding the current "depth" of dimension slicing. We progress by slicing dimensions from the end.
145  m_current_slice = m_max_dim;
146 }
147 
149 void wfg::free_wfg_members() const
150 {
151  // Free the memory.
152  delete[] m_refpoint;
153 
154  for(unsigned int fr_idx = 0 ; fr_idx < m_n_frames ; ++fr_idx) {
155  for(unsigned int p_idx = 0; p_idx < m_max_points ; ++p_idx) {
156  delete[] m_frames[fr_idx][p_idx];
157  }
158  delete[] m_frames[fr_idx];
159  }
160  delete[] m_frames;
161  delete[] m_frames_size;
162 }
163 
165 void wfg::limitset(const unsigned int begin_idx, const unsigned int p_idx, const unsigned int rec_level) const
166 {
167  double **points = m_frames[rec_level - 1];
168  unsigned int n_points = m_frames_size[rec_level - 1];
169 
170  int no_points = 0;
171 
172  double* p = points[p_idx];
173  double** frame = m_frames[rec_level];
174 
175  for(unsigned int idx = begin_idx; idx < n_points; ++idx) {
176  if (idx == p_idx) {
177  continue;
178  }
179 
180  for(fitness_vector::size_type f_idx = 0; f_idx < m_current_slice; ++f_idx) {
181  frame[no_points][f_idx] = std::max(points[idx][f_idx], p[f_idx]);
182  }
183 
184  std::vector<int> cmp_results;
185  cmp_results.resize(no_points);
186  double* s = frame[no_points];
187 
188  bool keep_s = true;
189 
190  // Check whether any point is dominating the point 's'.
191  for(int q_idx = 0; q_idx < no_points; ++q_idx) {
192  cmp_results[q_idx] = base::dom_cmp(s, frame[q_idx], m_current_slice);
193  if (cmp_results[q_idx] == base::DOM_CMP_B_DOMINATES_A) {
194  keep_s = false;
195  break;
196  }
197  }
198  // If neither is, remove points dominated by 's' (we store that during the first loop).
199  if( keep_s ) {
200  int prev = 0;
201  int next = 0;
202  while(next < no_points) {
203  if( cmp_results[next] != base::DOM_CMP_A_DOMINATES_B && cmp_results[next] != base::DOM_CMP_A_B_EQUAL) {
204  if(prev < next) {
205  for(unsigned int d_idx = 0; d_idx < m_current_slice ; ++d_idx) {
206  frame[prev][d_idx] = frame[next][d_idx];
207  }
208  }
209  ++prev;
210  }
211  ++next;
212  }
213  // Append 's' at the end, if prev==next it's not necessary as it's already there.
214  if(prev < next) {
215  for(unsigned int d_idx = 0; d_idx < m_current_slice ; ++d_idx) {
216  frame[prev][d_idx] = s[d_idx];
217  }
218  }
219  no_points = prev + 1;
220  }
221  }
222 
223  m_frames_size[rec_level] = no_points;
224 }
225 
227 double wfg::compute_hv(const unsigned int rec_level) const
228 {
229  double **points = m_frames[rec_level - 1];
230  unsigned int n_points = m_frames_size[rec_level - 1];
231 
232  // Simple inclusion-exclusion for one and two points
233  if (n_points == 1) {
234  return base::volume_between(points[0], m_refpoint, m_current_slice);
235  }
236  else if (n_points == 2) {
237  double hv = base::volume_between(points[0], m_refpoint, m_current_slice)
238  + base::volume_between(points[1], m_refpoint, m_current_slice);
239  double isect = 1.0;
240  for(unsigned int i=0;i<m_current_slice;++i) {
241  isect *= (m_refpoint[i] - std::max(points[0][i], points[1][i]));
242  }
243  return hv - isect;
244  }
245 
246  // If already sliced to dimension at which we use another algorithm.
247  if (m_current_slice == m_stop_dimension) {
248 
249  if (m_stop_dimension == 2) {
250  // Use a very efficient version of hv2d
251  return hv2d().compute(points, n_points, m_refpoint);
252  } else {
253  // Let hypervolume object pick the best method otherwise.
254  std::vector<fitness_vector> points_cpy;
255  points_cpy.reserve(n_points);
256  for(unsigned int i = 0 ; i < n_points ; ++i) {
257  points_cpy.push_back(fitness_vector(points[i], points[i] + m_current_slice));
258  }
259  fitness_vector r_cpy(m_refpoint, m_refpoint + m_current_slice);
260 
261  hypervolume hv = hypervolume(points_cpy, false);
262  hv.set_copy_points(false);
263  return hv.compute(r_cpy);
264  }
265  } else {
266  // Otherwise, sort the points in preparation for the next recursive step
267  // Bind the object under "this" pointer to the cmp_points method so it can be used as a valid comparator function for std::sort
268  // We need that in order for the cmp_points to have acces to the m_current_slice member variable.
269  std::sort(points, points + n_points, boost::bind(&wfg::cmp_points, this, _1, _2));
270  }
271 
272  double H = 0.0;
273  --m_current_slice;
274 
275  if(rec_level >= m_n_frames) {
276  double** fr = new double*[m_max_points];
277  for(unsigned int i = 0 ; i < m_max_points ; ++i) {
278  fr[i]=new double[m_current_slice];
279  }
280  m_frames[m_n_frames] = fr;
281  m_frames_size[m_n_frames] = 0;
282  ++m_n_frames;
283  }
284 
285  for(unsigned int p_idx = 0 ; p_idx < n_points ; ++p_idx) {
286  limitset(p_idx + 1, p_idx, rec_level);
287 
288  H += fabs((points[p_idx][m_current_slice] - m_refpoint[m_current_slice]) * exclusive_hv(p_idx, rec_level));
289  }
290  ++m_current_slice;
291  return H;
292 }
293 
295 double wfg::exclusive_hv(const unsigned int p_idx, const unsigned int rec_level) const
296 {
297  //double H = base::volume_between(points[p_idx], m_refpoint, m_current_slice);
298  double H = base::volume_between(m_frames[rec_level - 1][p_idx], m_refpoint, m_current_slice);
299 
300  if (m_frames_size[rec_level] == 1) {
301  H -= base::volume_between(m_frames[rec_level][0], m_refpoint, m_current_slice);
302  } else if (m_frames_size[rec_level] > 1) {
303  H -= compute_hv(rec_level + 1);
304  }
305 
306  return H;
307 }
308 
310 
318 void wfg::verify_before_compute(const std::vector<fitness_vector> &points, const fitness_vector &r_point) const
319 {
320  base::assert_minimisation(points, r_point);
321 }
322 
325 {
326  return base_ptr(new wfg(*this));
327 }
328 
330 std::string wfg::get_name() const
331 {
332  return "WFG algorithm";
333 }
334 
335 } } }
336 
337 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::util::hv_algorithm::wfg)
Root PaGMO namespace.
base_ptr clone() const
Clone method.
Definition: wfg.cpp:324
static double volume_between(const fitness_vector &, const fitness_vector &, unsigned int=0)
Compute volume between two points.
void assert_minimisation(const std::vector< fitness_vector > &, const fitness_vector &) const
Assert that reference point dominates every other point from the set.
static int dom_cmp(double *, double *, unsigned int)
Dominance comparison method.
void verify_before_compute(const std::vector< fitness_vector > &, const fitness_vector &) const
Verify before compute method.
Definition: wfg.cpp:318
wfg(const unsigned int stop_dimension=2)
Constructor.
Definition: wfg.cpp:51
WFG hypervolume algorithm.
Definition: wfg.h:48
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
double compute(std::vector< fitness_vector > &, const fitness_vector &) const
Compute hypervolume.
Definition: wfg.cpp:67
boost::shared_ptr< base > base_ptr
Base hypervolume algorithm class.
std::vector< double > contributions(std::vector< fitness_vector > &, const fitness_vector &) const
Contributions method.
Definition: wfg.cpp:88
std::string get_name() const
Algorithm name.
Definition: wfg.cpp:330