TSP.cpp
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief A 10-city traveling salesman problem.
5 
6  * The implementation follows:
7  * <p>
8  * D. E. Goldberg and R. Lingle, Alleles, loci, and traveling
9  * salesman problem. In <em> Proc. of the International Conference on
10  * Genetic Algorithms and Their Applications</em>, pages 154-159,
11  * Pittsburg, PA, 1985 </p>
12 
13  * The traveling salsman problem is a combinatorial optimization task. A
14  * salesman is supposed to visit $n$ cities. Each travelling connection
15  * is associated with a cost (i.e. the time fot the trip). The problem is
16  * to find the cheapest round-route that visits each city exactly once
17  * and returns to the starting point.
18 
19  *
20  *
21  * \author T. Voss
22  * \date -
23  *
24  *
25  * \par Copyright 1995-2015 Shark Development Team
26  *
27  * <BR><HR>
28  * This file is part of Shark.
29  * <http://image.diku.dk/shark/>
30  *
31  * Shark is free software: you can redistribute it and/or modify
32  * it under the terms of the GNU Lesser General Public License as published
33  * by the Free Software Foundation, either version 3 of the License, or
34  * (at your option) any later version.
35  *
36  * Shark is distributed in the hope that it will be useful,
37  * but WITHOUT ANY WARRANTY; without even the implied warranty of
38  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
39  * GNU Lesser General Public License for more details.
40  *
41  * You should have received a copy of the GNU Lesser General Public License
42  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
43  *
44  */
45 
47 
51 
52 #include <shark/LinAlg/Base.h>
53 
54 #include <boost/format.hpp>
55 #include <boost/graph/adjacency_matrix.hpp>
56 #include <boost/graph/adjacency_list.hpp>
57 #include <boost/graph/graphml.hpp>
58 #include <boost/graph/graphviz.hpp>
59 #include <boost/graph/metric_tsp_approx.hpp>
60 #include <boost/property_map/dynamic_property_map.hpp>
61 #include <boost/range/algorithm_ext/iota.hpp>
62 
63 typedef shark::IntVector Tour;
64 
65 /**
66  * \brief Defines the problem, i.e., its cost matrix.
67  */
68 const double cities[10][10] = {
69  { 0, 28, 57, 72, 81, 85, 80, 113, 89, 80 },
70  { 28, 0, 28, 45, 54, 57, 63, 85, 63, 63 },
71  { 57, 28, 0, 20, 30, 28, 57, 57, 40, 57 },
72  { 72, 45, 20, 0, 10, 20, 72, 45, 20, 45 },
73  { 81, 54, 30, 10, 0, 22, 81, 41, 10, 41 },
74  { 85, 57, 28, 20, 22, 0, 63, 28, 28, 63 },
75  { 80, 63, 57, 72, 81, 63, 0, 80, 89, 113 },
76  { 113, 85, 57, 45, 41, 28, 80, 0, 40, 80 },
77  { 89, 63, 40, 20, 10, 28, 89, 40, 0, 40 },
78  { 80, 63, 57, 45, 41, 63, 113, 80, 40, 0 }
79 };
80 
81 namespace shark {
82 
83  typedef boost::adjacency_matrix< boost::undirectedS,
84  boost::property< boost::vertex_color_t, std::string >,
85  boost::property< boost::edge_weight_t, double,
86  boost::property< boost::edge_color_t, std::string >
87  >,
88  boost::no_property
89  > Graph;
90  typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
91  typedef boost::graph_traits<Graph>::edge_descriptor Edge;
92  typedef boost::property_map<Graph, boost::edge_weight_t>::type WeightMap;
93  typedef boost::property_map<Graph, boost::edge_color_t>::type ColorMap;
94 
96  typedef std::vector< IndividualType > Population;
97 
98  bool compare_fitness( const IndividualType & a, const IndividualType & b ) {
99  return( a.unpenalizedFitness() < b.unpenalizedFitness() );
100  }
101 
102  /**
103  * \brief Calculates the cost of a tour w.r.t. to a cost matrix.
104  */
105  struct TspTourLength : public shark::SingleObjectiveFunction{
106  typedef shark::SingleObjectiveFunction base_type;
107 
108  /**
109  * \brief Default c'tor, initializes the cost matrix.
110  */
111  TspTourLength( const shark::RealMatrix & costMatrix = shark::RealMatrix() ) : m_costMatrix( costMatrix )
112  { }
113 
114  std::string name() const
115  { return "TspTourLength"; }
116 
117  std::size_t numberOfVariables() const
118  { return m_costMatrix.size1(); }
119 
120  /**
121  * \brief Calculates the costs of the supplied tour.
122  */
123  base_type::ResultType eval( const base_type::SearchPointType & input ) const {
124  SIZE_CHECK( input.size() == m_costMatrix.size1() );
125  m_evaluationCounter++;
126  base_type::ResultType result( 0 );
127  for( std::size_t i = 0; i < input.size() - 1; i++ ) {
128  result += m_costMatrix( input( i ), input( i+1 ) );
129  }
130  return( result );
131  }
132 
133  shark::RealMatrix m_costMatrix;
134  };
135 
136  /**
137  * @brief Implements partially mapped crossover
138  *
139  * Makes sure that only correct tours on the graph
140  * are generated.
141  */
142  struct PartiallyMappedCrossover {
143 
144  /**
145  * @brief Performs the crossover.
146  *
147  * @param [in] parent1 First parent individual.
148  * @param [in] parent2 Second parent individual.
149  * @returns A pair of offspring individuals.
150  */
151  template<typename IndividualType>
152  std::pair<IndividualType, IndividualType> operator()(
153  const IndividualType & parent1,
154  const IndividualType & parent2 ) const
155  {
156  std::pair< IndividualType, IndividualType > offspring( parent1, parent2 );
157 
158  const Tour & t1 = parent1.searchPoint();
159  const Tour & t2 = parent2.searchPoint();
160 
161  std::size_t cuttingPoint1 = shark::Rng::discrete( 0, t1.size() - 1 );
162  std::size_t cuttingPoint2 = shark::Rng::discrete( 0, t2.size() - 1 );
163 
164  while( cuttingPoint1 == cuttingPoint2 )
165  cuttingPoint2 = shark::Rng::discrete( 0, t2.size() - 1 );
166 
167  if( cuttingPoint1 > cuttingPoint2 )
168  std::swap( cuttingPoint1, cuttingPoint2 );
169 
170  Tour r1( t1.size(), -1 ), r2( t2.size(), -1 );
171 
172  for( std::size_t i = cuttingPoint1; i <= cuttingPoint2; i++ ) {
173  offspring.first.searchPoint()( i ) = t2( i );
174  offspring.second.searchPoint()( i ) = t1( i );
175 
176  r1[ t2( i ) ] = t1( i );
177  r2[ t1( i ) ] = t2( i );
178  }
179 
180  for( std::size_t i = 0; i < t1.size(); i++) {
181  if ((i >= cuttingPoint1) && (i <= cuttingPoint2)) continue;
182 
183  std::size_t n1 = t1[i] ;
184  std::size_t m1 = r1[n1] ;
185 
186  std::size_t n2 = t2[i] ;
187  std::size_t m2 = r2[n2] ;
188 
189  while (m1 != std::numeric_limits<std::size_t>::max()) {
190  n1 = m1 ;
191  m1 = r1[m1] ;
192  }
193  while (m2 != std::numeric_limits<std::size_t>::max()) {
194  n2 = m2 ;
195  m2 = r2[m2] ;
196  }
197  offspring.first.searchPoint()[i] = n1 ;
198  offspring.second.searchPoint()[i] = n2 ;
199  }
200 
201  return offspring;
202  }
203  };
204 }
205 
206 
207 int main( int argc, char ** argv ) {
208 
209  // Define the problem instance
210  shark::Graph g( 10 );
211  // Iterate the graph and assign a label to every vertex
212  boost::graph_traits<shark::Graph>::vertex_iterator v, v_end;
213  for( boost::tie(v,v_end) = boost::vertices(g); v != v_end; ++v )
214  boost::put(boost::vertex_color_t(), g, *v, ( boost::format( "City_%1%" ) % *v ).str() );
215  // Get hold of the weight map and the color map for calculation and visualization purposes.
216  shark::WeightMap weightMap = boost::get( boost::edge_weight, g );
217  shark::ColorMap colorMap = boost::get( boost::edge_color, g );
218 
219  // Iterate the graph and insert distances.
220  shark::Edge e;
221  bool inserted = false;
222 
223  shark::RealMatrix costMatrix( 10, 10 );
224  for( std::size_t i = 0; i < costMatrix.size1(); i++ ) {
225  for( std::size_t j = 0; j < costMatrix.size1(); j++ ) {
226 
227  if( i == j ) continue;
228 
229  costMatrix(i,j) = cities[i][j];
230 
231  // Add the edge
232  boost::tie( e, inserted ) = boost::add_edge( i, j, g );
233  if( inserted ) {
234  // Remember distance between cities i and j.
235  weightMap[ e ] = cities[i][j];
236  // Mark the edge as blue.
237  colorMap[ e ] = "blue";
238  }
239  }
240  }
241 
242  // Mating selection operator
244  // Variation (crossover) operator
245  shark::PartiallyMappedCrossover pmc;
246  // Fitness function instance.
247  shark::TspTourLength ttl( costMatrix );
248 
249  // Size of the parent population
250  const std::size_t mu = 100;
251  // Size of the offspring population
252  const std::size_t lambda = 100;
253 
254  // Parent population
255  shark::Population parents( mu );
256  // Offspring population
257  shark::Population offspring( lambda );
258 
259  // Default tour: 0,...,9
260  Tour t( 10 ); boost::iota( t, 0 );
261 
262  // Initialize parents
263  for( std::size_t i = 0; i < parents.size(); i++ ) {
264  parents[i].searchPoint() = t;
265  // Generate a permutation.
266  std::random_shuffle( parents[ i ].searchPoint().begin(), parents[ i ].searchPoint().end() );
267  // Evaluate the individual.
268  parents[i].penalizedFitness() = parents[i].unpenalizedFitness() = ttl( parents[i].searchPoint() );
269  }
270 
271  // Loop until maximum number of fitness function evaluations is reached.
272  while( ttl.evaluationCounter() < 10000 ) {
273  shark::RealVector selectionProbabilities(parents.size());
274  for(std::size_t i = 0; i != parents.size(); ++i){
275  selectionProbabilities(i) = parents[i].unpenalizedFitness();
276  }
277  selectionProbabilities/= sum(selectionProbabilities);
278  // Crossover candidates.
279  for( std::size_t i = 0; i < offspring.size() - 1; i+=2 ) {
280  // Carry out fitness proportional fitness selection and
281  // perform partially mapped crossover on parent individuals.
282  std::pair<
284  shark::IndividualType
285  > result = pmc(
286  *rws( parents.begin(), parents.end(), selectionProbabilities ),
287  *rws( parents.begin(), parents.end(), selectionProbabilities )
288  );
289  offspring[ i ] = result.first;
290  offspring[ i + 1 ] = result.second;
291 
292  // Evaluate offspring individuals.
293  offspring[ i ].penalizedFitness() =
294  offspring[ i ].unpenalizedFitness() = ttl( offspring[ i ].searchPoint() );
295 
296  offspring[ i+1 ].penalizedFitness() =
297  offspring[ i+1 ].unpenalizedFitness() = ttl( offspring[ i+1 ].searchPoint() );
298  }
299 
300  // Swap parent and offspring population.
301  std::swap( parents, offspring );
302  }
303 
304  // Sort in ascending order to find best individual.
305  std::sort( parents.begin(), parents.end(), shark::compare_fitness);
306 
307  // Extract the best tour currently known.
308  Tour final = parents.front().searchPoint();
309 
310  // Mark the final tour green in the graph.
311  bool extracted = false;
312  for( std::size_t i = 0; i < final.size() - 1; i++ ) {
313  boost::tie( e, extracted ) = boost::edge( shark::Vertex( final( i ) ), shark::Vertex( final( i + 1 ) ), g );
314 
315  if( extracted )
316  colorMap[ e ] = "green";
317  }
318 
319  // Calculate approximate solution
320  double len = 0.0;
321  std::vector< shark::Vertex > approxTour;
322  boost::metric_tsp_approx( g, boost::make_tsp_tour_len_visitor( g, std::back_inserter( approxTour ) , len, weightMap ) );
323  // Mark the approximation red in the graph.
324  for( std::size_t i = 0; i < approxTour.size() - 1; i++ ) {
325  boost::tie( e, extracted ) = boost::edge( approxTour[ i ], approxTour[ i+1 ], g );
326 
327  if( extracted )
328  colorMap[ e ] = "red";
329  }
330 
331  // Output the graph and the final path
332  std::ofstream outGraphviz( "graph.dot" );
333  boost::dynamic_properties dp;
334  dp.property( "node_id", boost::get( boost::vertex_color, g ) );
335  dp.property( "weight", boost::get( boost::edge_weight, g ) );
336  dp.property( "label", boost::get( boost::edge_weight, g ) );
337  dp.property( "color", boost::get( boost::edge_color, g ) );
338  boost::write_graphviz_dp( outGraphviz, g, dp );
339 
340  // Output best solution and corresponding fitness.
341  std::cout << parents.front().searchPoint() << " -> " << parents.front().unpenalizedFitness() << std::endl;
342  // Output approximate solution and corresponding fitness.
343  std::cout << "Approx: " << len << " vs. GA: " << parents.front().unpenalizedFitness() << std::endl;
344 
345  return( EXIT_SUCCESS );
346 }