PaGMO  1.1.5
tsp_cs.cpp
1 /*****************************************************************************
2  * Copyright (C) 2004-2014 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 #include "tsp_cs.h"
26 #include "../population.h"
27 
28 namespace pagmo { namespace problem {
29 
31 
36  tsp_cs::tsp_cs() : base_tsp(3, 0, 0 , base_tsp::RANDOMKEYS), m_weights(), m_values(), m_max_path_length(1.0)
37  {
38  std::vector<double> dumb(3,0);
39  m_weights = std::vector<std::vector<double> > (3,dumb);
40  m_weights[0][1] = 1;
41  m_weights[0][2] = 1;
42  m_weights[2][1] = 1;
43  m_weights[1][0] = 1;
44  m_weights[2][0] = 1;
45  m_weights[1][2] = 1;
46 
47  m_values = std::vector<double>(3,1.0);
48  m_max_edge_length = 1;
49  }
50 
52 
61  tsp_cs::tsp_cs(const std::vector<std::vector<double> >& weights, const std::vector<double>& values, const double max_path_length, const base_tsp::encoding_type & encoding):
62  base_tsp(weights.size(),
63  compute_dimensions(weights.size(), encoding)[0],
64  compute_dimensions(weights.size(), encoding)[1],
65  encoding
66  ), m_weights(weights), m_values(values), m_max_path_length(max_path_length)
67  {
68  check_weights(m_weights);
69  if (weights.size() != values.size())
70  {
71  pagmo_throw(value_error,"Size of weight matrix and values vector must be equal");
72  }
73 
74  m_max_edge_length = 0;
75  for (std::vector<std::vector<double> >::size_type i=0; i < weights.size(); ++i)
76  {
77  for (std::vector<std::vector<double> >::size_type j=0; j < weights.size(); ++j)
78  {
79  m_max_edge_length = (weights[i][j] > m_max_edge_length) ? (weights[i][j]) : (m_max_edge_length);
80  }
81  }
82  }
83 
86  {
87  return base_ptr(new tsp_cs(*this));
88  }
89 
91 
98  void tsp_cs::check_weights(const std::vector<std::vector<double> > &matrix) const
99  {
100  decision_vector::size_type n_cols = matrix.size();
101 
102  for (decision_vector::size_type i = 0; i < n_cols; ++i) {
103  decision_vector::size_type n_rows = matrix.at(i).size();
104  // check if the matrix is square
105  if (n_rows != n_cols)
106  pagmo_throw(value_error, "adjacency matrix is not square");
107 
108  for (size_t j = 0; j < n_rows; ++j) {
109  if (i == j && matrix.at(i).at(j) != 0)
110  pagmo_throw(value_error, "main diagonal elements must all be zeros.");
111  if (i != j && !matrix.at(i).at(j)) // fully connected
112  pagmo_throw(value_error, "adjacency matrix contains zero values.");
113  if (i != j && (!matrix.at(i).at(j)) == matrix.at(i).at(j)) // fully connected
114  pagmo_throw(value_error, "adjacency matrix contains NaN values.");
115  }
116  }
117  }
118 
119  boost::array<int, 2> tsp_cs::compute_dimensions(decision_vector::size_type n_cities, base_tsp::encoding_type encoding)
120  {
121  boost::array<int,2> retval;
122  switch( encoding ) {
123  case FULL:
124  retval[0] = n_cities*(n_cities-1)+2;
125  retval[1] = (n_cities-1)*(n_cities-2);
126  break;
127  case RANDOMKEYS:
128  retval[0] = 0;
129  retval[1] = 0;
130  break;
131  case CITIES:
132  retval[0] = 1;
133  retval[1] = 0;
134  break;
135  }
136  return retval;
137  }
138 
139  void tsp_cs::objfun_impl(fitness_vector &f, const decision_vector& x) const
140  {
141  f[0]=0;
142  decision_vector tour;
143  decision_vector::size_type n_cities = get_n_cities();
144  double cum_p, saved_length;
145  decision_vector::size_type dumb1, dumb2;
146 
147  switch( get_encoding() ) {
148  case FULL:
149  {
150  tour = full2cities(x);
151  break;
152  }
153  case RANDOMKEYS:
154  {
155  tour = randomkeys2cities(x);
156  break;
157  }
158  case CITIES:
159  {
160  tour = x;
161  break;
162  }
163  }
164  // We compute the maximal open walk
165  find_subsequence(tour, cum_p, saved_length, dumb1, dumb2);
166 
167  // We compute the length of the whole Hamiltonian path
168  double ham_path_len = 0;
169  for (decision_vector::size_type i=0; i<n_cities-1; ++i)
170  {
171  ham_path_len += m_weights[tour[i]][tour[i+1]];
172  }
173 
174  f[0] = -(cum_p) - (1 - ham_path_len / (n_cities * m_max_edge_length));
175  return;
176  }
177 
178  size_t tsp_cs::compute_idx(const size_t i, const size_t j, const size_t n) const
179  {
180  pagmo_assert( i!=j && i<n && j<n );
181  return i*(n-1) + j - (j>i? 1:0);
182  }
183 
184 
186 
199  void tsp_cs::find_subsequence(const decision_vector& tour, double& retval_p, double& retval_l, decision_vector::size_type& retval_it_l, decision_vector::size_type& retval_it_r) const
200  {
201  if (tour.size() != get_n_cities())
202  {
203  pagmo_throw(value_error, "tour dimension must be equal to the city number");
204  }
205 
206  // We declare the necessary variable
207  decision_vector::size_type n_cities = get_n_cities();
208  decision_vector::size_type it_l = 0, it_r = 0;
209  bool cond_r = true, cond_l = true;
210  double cum_p = m_values[tour[0]];
211  double saved_length = m_max_path_length;
212 
213  // We initialize the starting values
214  retval_p = cum_p;
215  retval_l = saved_length;
216  retval_it_l = it_l;
217  retval_it_r = it_r;
218 
219  // Main body of the double loop
220 
221  while(cond_l)
222  {
223  while(cond_r)
224  {
225  // We increment the right "pointer" updating the value and length of the path
226  saved_length -= m_weights[tour[it_r % n_cities]][tour[(it_r + 1) % n_cities]];
227  cum_p += m_values[tour[(it_r + 1) % n_cities]];
228  it_r += 1;
229 
230  // We update the various retvals only if the new subpath is valid
231  if (saved_length < 0 || (it_l % n_cities == it_r % n_cities))
232  {
233  cond_r = false;
234  }
235  else if (cum_p > retval_p)
236  {
237  retval_p = cum_p;
238  retval_l = saved_length;
239  retval_it_l = it_l % n_cities;
240  retval_it_r = it_r % n_cities;
241  }
242  else if (cum_p == retval_p)
243  {
244  if (saved_length > retval_l)
245  {
246  retval_p = cum_p;
247  retval_l = saved_length;
248  retval_it_l = it_l % n_cities;
249  retval_it_r = it_r % n_cities;
250  }
251  }
252  }
253  // We get out if all cities are included in the current path
254  if (it_l % n_cities == it_r % n_cities)
255  {
256  cond_l = false;
257  }
258  else
259  {
260  // We increment the left "pointer" updating the value and length of the path
261  saved_length += m_weights[tour[it_l % n_cities]][tour[(it_l + 1) % n_cities]];
262  cum_p -= m_values[tour[it_l]];
263  it_l += 1;
264  // We update the various retvals only if the new subpath is valid
265  if (saved_length > 0)
266  {
267  cond_r = true;
268  if (cum_p > retval_p)
269  {
270  retval_p = cum_p;
271  retval_l = saved_length;
272  retval_it_l = it_l % n_cities;
273  retval_it_r = it_r % n_cities;
274  }
275  else if (cum_p == retval_p)
276  {
277  if (saved_length > retval_l)
278  {
279  retval_p = cum_p;
280  retval_l = saved_length;
281  retval_it_l = it_l % n_cities;
282  retval_it_r = it_r % n_cities;
283  }
284  }
285  }
286  if (it_l == n_cities)
287  {
288  cond_l = false;
289  }
290  }
291  }
292  }
293 
294  void tsp_cs::compute_constraints_impl(constraint_vector &c, const decision_vector& x) const
295  {
296  decision_vector::size_type n_cities = get_n_cities();
297 
298  switch( get_encoding() )
299  {
300  case FULL:
301  {
302  // 1 - We set the equality constraints
303  for (size_t i = 0; i < n_cities; i++) {
304  c[i] = 0;
305  c[i+n_cities] = 0;
306  for (size_t j = 0; j < n_cities; j++) {
307  if(i==j) continue; // ignoring main diagonal
308  decision_vector::size_type rows = compute_idx(i, j, n_cities);
309  decision_vector::size_type cols = compute_idx(j, i, n_cities);
310  c[i] += x[rows];
311  c[i+n_cities] += x[cols];
312  }
313  c[i] = c[i]-1;
314  c[i+n_cities] = c[i+n_cities]-1;
315  }
316 
317  //2 - We set the inequality constraints
318  //2.1 - First we compute the uj (see http://en.wikipedia.org/wiki/Travelling_salesman_problem#Integer_linear_programming_formulation)
319  // we start always out tour from the first city, without loosing generality
320  size_t next_city = 0,current_city = 0;
321  std::vector<int> u(n_cities);
322  for (size_t i = 0; i < n_cities; i++) {
323  u[current_city] = i+1;
324  for (size_t j = 0; j < n_cities; j++)
325  {
326  if (current_city==j) continue;
327  if (x[compute_idx(current_city, j, n_cities)] == 1)
328  {
329  next_city = j;
330  break;
331  }
332  }
333  current_city = next_city;
334  }
335  int count=0;
336  for (size_t i = 1; i < n_cities; i++) {
337  for (size_t j = 1; j < n_cities; j++)
338  {
339  if (i==j) continue;
340  c[2*n_cities+count] = u[i]-u[j] + (n_cities+1) * x[compute_idx(i, j, n_cities)] - n_cities;
341  count++;
342  }
343  }
344  break;
345  }
346  case RANDOMKEYS:
347  break;
348  case CITIES:
349  {
350  std::vector<population::size_type> range(n_cities);
351  for (std::vector<population::size_type>::size_type i=0; i<range.size(); ++i)
352  {
353  range[i]=i;
354  }
355  c[0] = !std::is_permutation(x.begin(),x.end(),range.begin());
356  break;
357  }
358  }
359  return;
360  }
361 
363  double tsp_cs::distance(decision_vector::size_type i, decision_vector::size_type j) const
364  {
365  return m_weights[i][j];
366  }
367 
369 
372  const std::vector<std::vector<double> >& tsp_cs::get_weights() const
373  {
374  return m_weights;
375  }
376 
378 
381  const std::vector<double>& tsp_cs::get_values() const
382  {
383  return m_values;
384  }
385 
387 
391  {
392  return m_max_path_length;
393  }
394 
396  std::string tsp_cs::get_name() const
397  {
398  return "City-selection Travelling Salesman Problem (TSP-CS)";
399  }
400 
402 
405  std::string tsp_cs::human_readable_extra() const
406  {
407  std::ostringstream oss;
408  oss << "\n\tNumber of cities: " << get_n_cities() << '\n';
409  oss << "\tEncoding: ";
410  switch( get_encoding() ) {
411  case FULL:
412  oss << "FULL" << '\n';
413  break;
414  case RANDOMKEYS:
415  oss << "RANDOMKEYS" << '\n';
416  break;
417  case CITIES:
418  oss << "CITIES" << '\n';
419  break;
420  }
421  oss << "\tCities Values: " << m_values << std::endl;
422  oss << "\tMax path length: " << m_max_path_length << '\n';
423  oss << "\tWeight Matrix: \n";
424  for (decision_vector::size_type i=0; i<get_n_cities() ; ++i)
425  {
426  oss << "\t\t" << m_weights.at(i) << '\n';
427  if (i>5)
428  {
429  oss << "\t\t..." << '\n';
430  break;
431  }
432  }
433  return oss.str();
434  }
435 
436 
437 }} //namespaces
438 
439 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::problem::tsp_cs)
void find_subsequence(const decision_vector &, double &, double &, decision_vector::size_type &, decision_vector::size_type &) const
Computes the best subpath of an hamilonian path satisfying the max_path_length constraint.
Definition: tsp_cs.cpp:199
Root PaGMO namespace.
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base problem.
Definition: problem/base.h:62
std::vector< double > decision_vector
Decision vector type.
Definition: types.h:40
double distance(decision_vector::size_type, decision_vector::size_type) const
Definition of distance function.
Definition: tsp_cs.cpp:363
const std::vector< std::vector< double > > & get_weights() const
Getter for m_weights.
Definition: tsp_cs.cpp:372
const std::vector< double > & get_values() const
Getter for m_values.
Definition: tsp_cs.cpp:381
double get_max_path_length() const
Getter for m_max_path_value.
Definition: tsp_cs.cpp:390
The City-Selection Travelling Salesman Problem.
Definition: tsp_cs.h:57
tsp_cs()
Constructors.
Definition: tsp_cs.cpp:36
std::string human_readable_extra() const
Extra human readable info for the problem.
Definition: tsp_cs.cpp:405
base_ptr clone() const
Copy constructor for polymorphic objects.
Definition: tsp_cs.cpp:85
pagmo::decision_vector full2cities(const pagmo::decision_vector &) const
From FULL to CITIES encoding.
Definition: base_tsp.cpp:74
As a vector of doubles in [0,1].
Definition: base_tsp.h:69
decision_vector::size_type get_n_cities() const
Getter for m_n_cities.
Definition: base_tsp.cpp:193
encoding_type
Mechanism used to encode the sequence of vertices to be visited.
Definition: base_tsp.h:68
std::string get_name() const
Returns the problem name.
Definition: tsp_cs.cpp:396
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
std::vector< double > constraint_vector
Constraint vector type.
Definition: types.h:44
As a matrix with ones and zeros.
Definition: base_tsp.h:70
Base TSP (Travelling Salesman Problem).
Definition: base_tsp.h:64
As a sequence of cities ids.
Definition: base_tsp.h:71
pagmo::decision_vector randomkeys2cities(const pagmo::decision_vector &) const
From RANDOMKEYS to CITIES encoding.
Definition: base_tsp.cpp:127
decision_vector::size_type size_type
Problem's size type: the same as pagmo::decision_vector's size type.
Definition: problem/base.h:160
encoding_type get_encoding() const
Getter for m_encoding.
Definition: base_tsp.cpp:184