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 #if BOOST_VERSION == 106000
56 #include <boost/type_traits/ice.hpp>//Required because of boost 1.60 bug
57 #endif
58 #include <boost/graph/adjacency_matrix.hpp>
59 #include <boost/graph/adjacency_list.hpp>
60 #include <boost/graph/graphml.hpp>
61 #include <boost/graph/graphviz.hpp>
62 #include <boost/graph/metric_tsp_approx.hpp>
63 #include <boost/property_map/dynamic_property_map.hpp>
64 #include <boost/range/algorithm_ext/iota.hpp>
65 
66 typedef shark::IntVector Tour;
67 
68 /**
69  * \brief Defines the problem, i.e., its cost matrix.
70  */
71 const double cities[10][10] = {
72  { 0, 28, 57, 72, 81, 85, 80, 113, 89, 80 },
73  { 28, 0, 28, 45, 54, 57, 63, 85, 63, 63 },
74  { 57, 28, 0, 20, 30, 28, 57, 57, 40, 57 },
75  { 72, 45, 20, 0, 10, 20, 72, 45, 20, 45 },
76  { 81, 54, 30, 10, 0, 22, 81, 41, 10, 41 },
77  { 85, 57, 28, 20, 22, 0, 63, 28, 28, 63 },
78  { 80, 63, 57, 72, 81, 63, 0, 80, 89, 113 },
79  { 113, 85, 57, 45, 41, 28, 80, 0, 40, 80 },
80  { 89, 63, 40, 20, 10, 28, 89, 40, 0, 40 },
81  { 80, 63, 57, 45, 41, 63, 113, 80, 40, 0 }
82 };
83 
84 namespace shark {
85 
86  typedef boost::adjacency_matrix< boost::undirectedS,
87  boost::property< boost::vertex_color_t, std::string >,
88  boost::property< boost::edge_weight_t, double,
89  boost::property< boost::edge_color_t, std::string >
90  >,
91  boost::no_property
92  > Graph;
93  typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
94  typedef boost::graph_traits<Graph>::edge_descriptor Edge;
95  typedef boost::property_map<Graph, boost::edge_weight_t>::type WeightMap;
96  typedef boost::property_map<Graph, boost::edge_color_t>::type ColorMap;
97 
99  typedef std::vector< IndividualType > Population;
100 
101  bool compare_fitness( const IndividualType & a, const IndividualType & b ) {
102  return( a.unpenalizedFitness() < b.unpenalizedFitness() );
103  }
104 
105  /**
106  * \brief Calculates the cost of a tour w.r.t. to a cost matrix.
107  */
108  struct TspTourLength : public shark::SingleObjectiveFunction{
109  typedef shark::SingleObjectiveFunction base_type;
110 
111  /**
112  * \brief Default c'tor, initializes the cost matrix.
113  */
114  TspTourLength( const shark::RealMatrix & costMatrix = shark::RealMatrix() ) : m_costMatrix( costMatrix )
115  { }
116 
117  std::string name() const
118  { return "TspTourLength"; }
119 
120  std::size_t numberOfVariables() const
121  { return m_costMatrix.size1(); }
122 
123  /**
124  * \brief Calculates the costs of the supplied tour.
125  */
126  base_type::ResultType eval( const base_type::SearchPointType & input ) const {
127  SIZE_CHECK( input.size() == m_costMatrix.size1() );
128  m_evaluationCounter++;
129  base_type::ResultType result( 0 );
130  for( std::size_t i = 0; i < input.size() - 1; i++ ) {
131  result += m_costMatrix( input( i ), input( i+1 ) );
132  }
133  return( result );
134  }
135 
136  shark::RealMatrix m_costMatrix;
137  };
138 
139  /**
140  * @brief Implements partially mapped crossover
141  *
142  * Makes sure that only correct tours on the graph
143  * are generated.
144  */
145  struct PartiallyMappedCrossover {
146 
147  /**
148  * @brief Performs the crossover.
149  *
150  * @param [in] parent1 First parent individual.
151  * @param [in] parent2 Second parent individual.
152  * @returns A pair of offspring individuals.
153  */
154  template<typename IndividualType>
155  std::pair<IndividualType, IndividualType> operator()(
156  const IndividualType & parent1,
157  const IndividualType & parent2 ) const
158  {
159  std::pair< IndividualType, IndividualType > offspring( parent1, parent2 );
160 
161  const Tour & t1 = parent1.searchPoint();
162  const Tour & t2 = parent2.searchPoint();
163 
164  std::size_t cuttingPoint1 = shark::Rng::discrete( 0, t1.size() - 1 );
165  std::size_t cuttingPoint2 = shark::Rng::discrete( 0, t2.size() - 1 );
166 
167  while( cuttingPoint1 == cuttingPoint2 )
168  cuttingPoint2 = shark::Rng::discrete( 0, t2.size() - 1 );
169 
170  if( cuttingPoint1 > cuttingPoint2 )
171  std::swap( cuttingPoint1, cuttingPoint2 );
172 
173  Tour r1( t1.size(), -1 ), r2( t2.size(), -1 );
174 
175  for( std::size_t i = cuttingPoint1; i <= cuttingPoint2; i++ ) {
176  offspring.first.searchPoint()( i ) = t2( i );
177  offspring.second.searchPoint()( i ) = t1( i );
178 
179  r1[ t2( i ) ] = t1( i );
180  r2[ t1( i ) ] = t2( i );
181  }
182 
183  for( std::size_t i = 0; i < t1.size(); i++) {
184  if ((i >= cuttingPoint1) && (i <= cuttingPoint2)) continue;
185 
186  std::size_t n1 = t1[i] ;
187  std::size_t m1 = r1[n1] ;
188 
189  std::size_t n2 = t2[i] ;
190  std::size_t m2 = r2[n2] ;
191 
192  while (m1 != std::numeric_limits<std::size_t>::max()) {
193  n1 = m1 ;
194  m1 = r1[m1] ;
195  }
196  while (m2 != std::numeric_limits<std::size_t>::max()) {
197  n2 = m2 ;
198  m2 = r2[m2] ;
199  }
200  offspring.first.searchPoint()[i] = n1 ;
201  offspring.second.searchPoint()[i] = n2 ;
202  }
203 
204  return offspring;
205  }
206  };
207 }
208 
209 
210 int main( int argc, char ** argv ) {
211 
212  // Define the problem instance
213  shark::Graph g( 10 );
214  // Iterate the graph and assign a label to every vertex
215  boost::graph_traits<shark::Graph>::vertex_iterator v, v_end;
216  for( boost::tie(v,v_end) = boost::vertices(g); v != v_end; ++v )
217  boost::put(boost::vertex_color_t(), g, *v, ( boost::format( "City_%1%" ) % *v ).str() );
218  // Get hold of the weight map and the color map for calculation and visualization purposes.
219  shark::WeightMap weightMap = boost::get( boost::edge_weight, g );
220  shark::ColorMap colorMap = boost::get( boost::edge_color, g );
221 
222  // Iterate the graph and insert distances.
223  shark::Edge e;
224  bool inserted = false;
225 
226  shark::RealMatrix costMatrix( 10, 10 );
227  for( std::size_t i = 0; i < costMatrix.size1(); i++ ) {
228  for( std::size_t j = 0; j < costMatrix.size1(); j++ ) {
229 
230  if( i == j ) continue;
231 
232  costMatrix(i,j) = cities[i][j];
233 
234  // Add the edge
235  boost::tie( e, inserted ) = boost::add_edge( i, j, g );
236  if( inserted ) {
237  // Remember distance between cities i and j.
238  weightMap[ e ] = cities[i][j];
239  // Mark the edge as blue.
240  colorMap[ e ] = "blue";
241  }
242  }
243  }
244 
245  // Mating selection operator
247  // Variation (crossover) operator
248  shark::PartiallyMappedCrossover pmc;
249  // Fitness function instance.
250  shark::TspTourLength ttl( costMatrix );
251 
252  // Size of the parent population
253  const std::size_t mu = 100;
254  // Size of the offspring population
255  const std::size_t lambda = 100;
256 
257  // Parent population
258  shark::Population parents( mu );
259  // Offspring population
260  shark::Population offspring( lambda );
261 
262  // Default tour: 0,...,9
263  Tour t( 10 ); boost::iota( t, 0 );
264 
265  // Initialize parents
266  for( std::size_t i = 0; i < parents.size(); i++ ) {
267  parents[i].searchPoint() = t;
268  // Generate a permutation.
269  std::random_shuffle( parents[ i ].searchPoint().begin(), parents[ i ].searchPoint().end() );
270  // Evaluate the individual.
271  parents[i].penalizedFitness() = parents[i].unpenalizedFitness() = ttl( parents[i].searchPoint() );
272  }
273 
274  // Loop until maximum number of fitness function evaluations is reached.
275  while( ttl.evaluationCounter() < 10000 ) {
276  shark::RealVector selectionProbabilities(parents.size());
277  for(std::size_t i = 0; i != parents.size(); ++i){
278  selectionProbabilities(i) = parents[i].unpenalizedFitness();
279  }
280  selectionProbabilities/= sum(selectionProbabilities);
281  // Crossover candidates.
282  for( std::size_t i = 0; i < offspring.size() - 1; i+=2 ) {
283  // Carry out fitness proportional fitness selection and
284  // perform partially mapped crossover on parent individuals.
285  std::pair<
287  shark::IndividualType
288  > result = pmc(
289  *rws( shark::Rng::globalRng, parents.begin(), parents.end(), selectionProbabilities ),
290  *rws( shark::Rng::globalRng, parents.begin(), parents.end(), selectionProbabilities )
291  );
292  offspring[ i ] = result.first;
293  offspring[ i + 1 ] = result.second;
294 
295  // Evaluate offspring individuals.
296  offspring[ i ].penalizedFitness() =
297  offspring[ i ].unpenalizedFitness() = ttl( offspring[ i ].searchPoint() );
298 
299  offspring[ i+1 ].penalizedFitness() =
300  offspring[ i+1 ].unpenalizedFitness() = ttl( offspring[ i+1 ].searchPoint() );
301  }
302 
303  // Swap parent and offspring population.
304  std::swap( parents, offspring );
305  }
306 
307  // Sort in ascending order to find best individual.
308  std::sort( parents.begin(), parents.end(), shark::compare_fitness);
309 
310  // Extract the best tour currently known.
311  Tour final = parents.front().searchPoint();
312 
313  // Mark the final tour green in the graph.
314  bool extracted = false;
315  for( std::size_t i = 0; i < final.size() - 1; i++ ) {
316  boost::tie( e, extracted ) = boost::edge( shark::Vertex( final( i ) ), shark::Vertex( final( i + 1 ) ), g );
317 
318  if( extracted )
319  colorMap[ e ] = "green";
320  }
321 
322  // Calculate approximate solution
323  double len = 0.0;
324  std::vector< shark::Vertex > approxTour;
325  boost::metric_tsp_approx( g, boost::make_tsp_tour_len_visitor( g, std::back_inserter( approxTour ) , len, weightMap ) );
326  // Mark the approximation red in the graph.
327  for( std::size_t i = 0; i < approxTour.size() - 1; i++ ) {
328  boost::tie( e, extracted ) = boost::edge( approxTour[ i ], approxTour[ i+1 ], g );
329 
330  if( extracted )
331  colorMap[ e ] = "red";
332  }
333 
334  // Output the graph and the final path
335  std::ofstream outGraphviz( "graph.dot" );
336  boost::dynamic_properties dp;
337  dp.property( "node_id", boost::get( boost::vertex_color, g ) );
338  dp.property( "weight", boost::get( boost::edge_weight, g ) );
339  dp.property( "label", boost::get( boost::edge_weight, g ) );
340  dp.property( "color", boost::get( boost::edge_color, g ) );
341  boost::write_graphviz_dp( outGraphviz, g, dp );
342 
343  // Output best solution and corresponding fitness.
344  std::cout << parents.front().searchPoint() << " -> " << parents.front().unpenalizedFitness() << std::endl;
345  // Output approximate solution and corresponding fitness.
346  std::cout << "Approx: " << len << " vs. GA: " << parents.front().unpenalizedFitness() << std::endl;
347 
348  return( EXIT_SUCCESS );
349 }