ROOT  6.06/09
Reference Guide
simanTSP.cxx
Go to the documentation of this file.
1 // example of simulated annealing
2 // traveling salesman problem
3 // example from GSL siman_tsp, see also
4 // http://www.gnu.org/software/gsl/manual/html_node/Traveling-Salesman-Problem.html
5 //
6 //
7 // minimize total distance when visiting a set of cities
8 //
9 // you can run in ROOT by doing
10 //
11 //.x simanTSP.cxx+
12 // running FulLSearch() you can check that the result is correct
13 //
14 #include <cmath>
15 #include <vector>
16 #include <algorithm>
17 
18 #include "Math/GSLSimAnnealing.h"
19 #include "Math/GSLRndmEngines.h"
20 #include "Math/SMatrix.h"
21 #include "Math/Math.h"
22 #include "TH1.h"
23 #include "TGraph.h"
24 #include "TCanvas.h"
25 #include "TApplication.h"
26 
27 bool showGraphics = false;
28 
29 using namespace ROOT::Math;
30 
31 /* in this table, latitude and longitude are obtained from the US
32  Census Bureau, at http://www.census.gov/cgi-bin/gazetteer */
33 
34 struct s_tsp_city {
35  const char * name;
36  double lat, longitude; /* coordinates */
37 };
38 typedef struct s_tsp_city Stsp_city;
39 
40 Stsp_city cities[] = {{"Santa Fe", 35.68, 105.95},
41  {"Phoenix", 33.54, 112.07},
42  {"Albuquerque", 35.12, 106.62},
43  {"Clovis", 34.41, 103.20},
44  {"Durango", 37.29, 107.87},
45  {"Dallas", 32.79, 96.77},
46  {"Tesuque", 35.77, 105.92},
47  {"Grants", 35.15, 107.84},
48  {"Los Alamos", 35.89, 106.28},
49  {"Las Cruces", 32.34, 106.76},
50  {"Cortez", 37.35, 108.58},
51  {"Gallup", 35.52, 108.74}};
52 
53 #define N_CITIES (sizeof(cities)/sizeof(Stsp_city))
54 
56 
57 /* distance between two cities */
59 {
60  const double earth_radius = 6375.000; /* 6000KM approximately */
61  /* sin and std::cos of lat and long; must convert to radians */
62  double sla1 = std::sin(c1.lat*M_PI/180), cla1 = std::cos(c1.lat*M_PI/180),
63  slo1 = std::sin(c1.longitude*M_PI/180), clo1 = std::cos(c1.longitude*M_PI/180);
64  double sla2 = std::sin(c2.lat*M_PI/180), cla2 = std::cos(c2.lat*M_PI/180),
65  slo2 = std::sin(c2.longitude*M_PI/180), clo2 = std::cos(c2.longitude*M_PI/180);
66 
67  double x1 = cla1*clo1;
68  double x2 = cla2*clo2;
69 
70  double y1 = cla1*slo1;
71  double y2 = cla2*slo2;
72 
73  double z1 = sla1;
74  double z2 = sla2;
75 
76  double dot_product = x1*x2 + y1*y2 + z1*z2;
77 
78  double angle = std::acos(dot_product);
79 
80  /* distance is the angle (in radians) times the earth radius */
81  return angle*earth_radius;
82 }
83 
84 
86 {
87  unsigned int i, j;
88 
89  for (i = 0; i < N_CITIES; ++i) {
90  printf("# ");
91  for (j = 0; j < N_CITIES; ++j) {
92  printf("%15.8f ", distance_matrix[i][j]);
93  }
94  printf("\n");
95  }
96 }
97 
98 // re-implement GSLSimAnFunc for re-defining the metric and function
99 
100 
101 
102 class MySimAnFunc : public GSLSimAnFunc {
103 
104 public:
105 
106  MySimAnFunc( std::vector<double> & allDist)
107  {
108  calculate_distance_matrix();
109  // initial route is just the sequantial order
110  for (unsigned int i = 0; i < N_CITIES; ++i)
111  fRoute[i] = i;
112 
113  // keep track of all the found distances
114  fDist = &allDist;
115  }
116 
117 
118  virtual ~MySimAnFunc() {}
119 
120  unsigned int Route(unsigned int i) const { return fRoute[i]; }
121 
122  const unsigned int * Route() const { return fRoute; }
123  unsigned int * Route() { return fRoute; }
124 
125  virtual MySimAnFunc * Clone() const { return new MySimAnFunc(*this); }
126 
127  std::vector<double> & AllDist() { return *fDist; }
128 
129  virtual double Energy() const {
130  // calculate the energy
131 
132 
133  double enrg = 0;
134  for (unsigned int i = 0; i < N_CITIES; ++i) {
135  /* use the distance_matrix to optimize this calculation; it had
136  better be allocated!! */
137  enrg += fDistanceMatrix( fRoute[i] , fRoute[ (i + 1) % N_CITIES ] );
138  }
139 
140  //std::cout << "energy is " << enrg << std::endl;
141  return enrg;
142  }
143 
144  virtual double Distance(const GSLSimAnFunc & f) const {
145  const MySimAnFunc * f2 = dynamic_cast<const MySimAnFunc *> (&f);
146  assert (f2 != 0);
147  double d = 0;
148  // use change in permutations
149  for (unsigned int i = 0; i < N_CITIES; ++i) {
150  d += (( fRoute[i] == f2->Route(i) ) ? 0 : 1 );
151  }
152  return d;
153  }
154  virtual void Step(const GSLRandomEngine & r, double ) {
155  // swap to city in the matrix
156  int x1, x2, dummy;
157 
158  /* pick the two cities to swap in the matrix; we leave the first
159  city fixed */
160  x1 = r.RndmInt(N_CITIES);
161  do {
162  x2 = r.RndmInt(N_CITIES);
163  } while (x2 == x1);
164 
165  // swap x1 and x2
166  dummy = fRoute[x1];
167  fRoute[x1] = fRoute[x2];
168  fRoute[x2] = dummy;
169 
170  //std::cout << "make step -swap " << x1 << " " << x2 << std::endl;
171 
172  }
173 
174  virtual void Print() {
175  printf(" [");
176  for (unsigned i = 0; i < N_CITIES; ++i) {
177  printf(" %d ", fRoute[i]);
178  }
179  printf("] ");
180  fDist->push_back(Energy()); // store all found distances
181  }
182 
183  // fast copy (need to keep base class type for using virtuality
184  virtual MySimAnFunc & FastCopy(const GSLSimAnFunc & f) {
185  const MySimAnFunc * rhs = dynamic_cast<const MySimAnFunc *>(&f);
186  assert (rhs != 0);
187  std::copy(rhs->fRoute, rhs->fRoute + N_CITIES, fRoute);
188  return *this;
189  }
190 
191  double Distance(int i, int j) const { return fDistanceMatrix(i,j); }
192 
193  void PrintRoute(); // used for debugging at the end
194 
195  void SetRoute(unsigned int * r) { std::copy(r,r+N_CITIES,fRoute); } // set a new route (used by exh. search)
196 
197 private:
198 
199  void calculate_distance_matrix();
200 
201  // variable of the system - order how cities are visited
202  unsigned int fRoute[N_CITIES];
204  Matrix fDistanceMatrix;
205 
206  std::vector<double> * fDist; // pointer to all distance vector
207 
208 };
209 
210 
211 // calculate distance between the cities
212 void MySimAnFunc::calculate_distance_matrix()
213 {
214  unsigned int i, j;
215  double dist;
216 
217  for (i = 0; i < N_CITIES; ++i) {
218  for (j = 0; j <= i; ++j) {
219  if (i == j) {
220  dist = 0;
221  } else {
222  dist = city_distance(cities[i], cities[j]);
223  }
224  fDistanceMatrix(i,j) = dist;
225  }
226  }
227 }
228 
229 void MySimAnFunc::PrintRoute() {
230  // print the route and distance
231  double dtot = 0;
232  for (unsigned int i = 0; i < N_CITIES; ++i) {
233  std::cout << std::setw(20) << cities[Route(i)].name << " \t " << Route(i);
234  int j = Route(i);
235  int k = Route( (i+ 1) % N_CITIES );
236  dtot += Distance(j,k);
237  std::cout << "\tdistance [" << j << " , " << k << " ]\t= " << Distance(j,k) << "\tTotal Distance\t = " << dtot << std::endl;
238  }
239  std::cout << "Total Route energy is " << dtot << std::endl;
240 }
241 
242 
243 // minimize using simulated annealing
244 void simanTSP(bool debug = true) {
245 
246  std::vector<double> allDist;
247  allDist.reserve(5000); // keep track of all distance for histogramming later on
248 
249  // create class
250  MySimAnFunc f(allDist);
251 
252  GSLSimAnnealing siman;
253 
254  GSLSimAnParams & p = siman.Params();
255  p.n_tries = 200;
256  p.iters_fixed_T = 10;
257  p.step_size = 1; // not used
258  p.k = 1;
259  p.t_initial = 5000;
260  p.mu = 1.01;
261  p.t_min = 0.5;
262 
263  // set the parameters
264 
265  // solve
266  siman.Solve(f, debug);
267 
268  unsigned int niter = allDist.size();
269  std::cout << "minimum found is for distance " << f.Energy() << std::endl;
270  if (debug) std::cout << "number of iterations is " << niter << std::endl;
271 
272  std::cout << "Best Route is \n";
273  f.PrintRoute();
274 
275  // plot configurations
276  double x0[N_CITIES+1];
277  double y0[N_CITIES+1];
278  double xmin[N_CITIES+1];
279  double ymin[N_CITIES+1];
280 
281 
282  // plot histograms with distances
283  TH1 * h1 = new TH1D("h1","Distance",niter+1,0,niter+1);
284  for (unsigned int i = 1; i <=niter; ++i) {
285  h1->SetBinContent(i,allDist[i]);
286  }
287 
288  for (unsigned int i = 0; i < N_CITIES; ++i) {
289  // invert long to have west on left side
290  x0[i] = -cities[i].longitude;
291  y0[i] = cities[i].lat;
292  xmin[i] = -cities[f.Route(i)].longitude;
293  ymin[i] = cities[f.Route(i)].lat;
294  }
295  // to close the curve
296  x0[N_CITIES] = x0[0];
297  y0[N_CITIES] = y0[0];
298  xmin[N_CITIES] = xmin[0];
299  ymin[N_CITIES] = ymin[0];
300 
301  if ( showGraphics )
302  {
303 
304  TGraph *g0 = new TGraph(N_CITIES+1,x0,y0);
305  TGraph *gmin = new TGraph(N_CITIES+1,xmin,ymin);
306 
307  TCanvas * c1 = new TCanvas("c","TSP",10,10,1000,800);
308  c1->Divide(2,2);
309  c1->cd(1);
310  g0->Draw("alp");
311  g0->SetMarkerStyle(20);
312  c1->cd(2);
313  gmin->SetMarkerStyle(20);
314  gmin->Draw("alp");
315  c1->cd(3);
316  h1->Draw();
317  }
318 
319 }
320 
321 unsigned int r1[N_CITIES];
322 unsigned int r2[N_CITIES];
323 unsigned int r3[N_CITIES];
324 unsigned int nfiter = 0;
326 
327 // use STL algorithms for permutations
328 // algorithm does also the inverse
329 void do_all_perms(MySimAnFunc & f, int offset )
330 {
331  unsigned int * r = f.Route();
332 
333  // do all permutation of N_CITIES -1 (keep first one at same place)
334  while (std::next_permutation(r+offset,r+N_CITIES) ) {
335 
336  double E = f.Energy();
337  nfiter++;
338  /* now save the best 3 energies and routes */
339  if (E < best_E) {
340  third_E = second_E;
341  std::copy(r2,r2+N_CITIES,r3);
342  second_E = best_E;
343  std::copy(r1,r1+N_CITIES,r2);
344  best_E = E;
345  std::copy(r,r+N_CITIES,r1);
346  } else if (E < second_E) {
347  third_E = second_E;
348  std::copy(r2,r2+N_CITIES,r3);
349  second_E = E;
350  std::copy(r,r+N_CITIES,r2);
351  } else if (E < third_E) {
352  third_E = E;
353  std::copy(r,r+N_CITIES,r3);
354  }
355  }
356 }
357 
358 #ifdef O
359 /* James Theiler's recursive algorithm for generating all routes */
360 // version used in GSL
361 void do_all_perms(MySimAnFunc & f, int n) {
362 
363  if (n == (N_CITIES-1)) {
364  /* do it! calculate the energy/cost for that route */
365  const unsigned int * r = f.Route();
366  /* now save the best 3 energies and routes */
367  if (E < best_E) {
368  third_E = second_E;
369  std::copy(r2,r2+N_CITIES,r3);
370  second_E = best_E;
371  std::copy(r1,r1+N_CITIES,r2);
372  best_E = E;
373  std::copy(r,r+N_CITIES,r1);
374  } else if (E < second_E) {
375  third_E = second_E;
376  std::copy(r2,r2+N_CITIES,r3);
377  second_E = E;
378  std::copy(r,r+N_CITIES,r2);
379  } else if (E < third_E) {
380  third_E = E;
381  std::copy(r,r+N_CITIES,r3);
382  }
383  } else {
384  unsigned int newr[N_CITIES];
385  unsigned int j;
386  int swap_tmp;
387  const unsigned int * r = f.Route();
388  std::copy(r,r+N_CITIES,newr);
389  for (j = n; j < N_CITIES; ++j) {
390  // swap j and n
391  swap_tmp = newr[j];
392  newr[j] = newr[n];
393  newr[n] = swap_tmp;
394  f.SetRoute(newr);
395  do_all_perms(f, n+1);
396  }
397  }
398 }
399 #endif
400 
401 
402 // search by brute force the best route
403 // the full permutations will be Factorial (N_CITIES-1)
404 // which is approx 4 E+7 for 12 cities
405 
407 {
408  // use still MySimAnFunc for initial routes and distance
409  std::vector<double> dummy;
410 
411  MySimAnFunc f(dummy);
412 
413  // intitial config
414 
415  const unsigned int * r = f.Route();
416  std::copy(r,r+N_CITIES,r1);
417  std::copy(r,r+N_CITIES,r2);
418  std::copy(r,r+N_CITIES,r3);
419  best_E = f.Energy();
420  second_E = 1.E100;
421  third_E = 1.E100;
422 
423 
424  do_all_perms(f, 1 );
425 
426  std::cout << "number of calls " << nfiter << std::endl;
427 
428 
429  printf("\n# exhaustive best route: ");
430  std::cout << best_E << std::endl;
431  f.SetRoute(r1); f.PrintRoute();
432 
433 
434  printf("\n# exhaustive second_best route: ");
435  std::cout << second_E << std::endl;
436  f.SetRoute(r2); f.PrintRoute();
437 
438  printf("\n# exhaustive third_best route: ");
439  std::cout << third_E << std::endl;
440  f.SetRoute(r3); f.PrintRoute();
441 
442 }
443 
444 
445 
446 
447 #ifndef __CINT__
448 int main(int argc, char **argv)
449 {
450  using std::cout;
451  using std::endl;
452  using std::cerr;
453 
454  bool verbose = false;
455 
456  // Parse command line arguments
457  for (Int_t i=1 ; i<argc ; i++) {
458  std::string arg = argv[i] ;
459  if (arg == "-g") {
460  showGraphics = true;
461  }
462  if (arg == "-v") {
463  showGraphics = true;
464  verbose = true;
465  }
466  if (arg == "-h") {
467  cerr << "Usage: " << argv[0] << " [-g] [-v]\n";
468  cerr << " where:\n";
469  cerr << " -g : graphics mode\n";
470  cerr << " -v : verbose mode";
471  cerr << endl;
472  return -1;
473  }
474  }
475 
476  if ( showGraphics )
477  {
478  TApplication* theApp = 0;
479  theApp = new TApplication("App",&argc,argv);
480  theApp->Run();
481  simanTSP(verbose);
482  delete theApp;
483  theApp = 0;
484  }
485  else
486  {
487  simanTSP(verbose);
488  }
489 
490  // to check that the result is correct
491  // FullSearch();
492  return 0;
493 }
494 #endif
495 
GSLSimAnnealing class for performing a simulated annealing search of a multidimensional function...
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
float xmin
Definition: THbookFile.cxx:93
#define N_CITIES
Definition: simanTSP.cxx:53
TCanvas * c1
Definition: legend1.C:2
float ymin
Definition: THbookFile.cxx:93
#define assert(cond)
Definition: unittest.h:542
GSLSimAnFunc class description.
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:659
void Step(const gsl_rng *r, void *xp, double step_size)
int Int_t
Definition: RtypesCore.h:41
struct s_tsp_city Stsp_city
Definition: simanTSP.cxx:38
Stsp_city cities[]
Definition: simanTSP.cxx:40
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:740
THist< 1, double > TH1D
Definition: THist.h:314
double cos(double)
int Solve(const ROOT::Math::IMultiGenFunction &func, const double *x0, const double *scale, double *xmin, bool debug=false)
solve the simulated annealing given a multi-dim function, the initial vector parameters and a vector ...
static const double x2[5]
SMatrix: a generic fixed size D1 x D2 Matrix class.
double acos(double)
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
double k
parameters for the Boltzman distribution
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
double third_E
Definition: simanTSP.cxx:325
TH1F * h1
Definition: legend1.C:5
double sin(double)
void FullSearch()
Definition: simanTSP.cxx:406
unsigned int RndmInt(unsigned int max) const
Generate an integer number between [0,max-1] (including 0 and max-1) if max is larger than available ...
int main(int argc, char **argv)
Definition: simanTSP.cxx:448
void simanTSP(bool debug=true)
Definition: simanTSP.cxx:244
void do_all_perms(MySimAnFunc &f, int offset)
Definition: simanTSP.cxx:329
#define M_PI
Definition: Rotated.cxx:105
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2878
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8543
bool verbose
Double_t E()
Definition: TMath.h:54
virtual void SetMarkerStyle(Style_t mstyle=1)
Definition: TAttMarker.h:53
The Canvas class.
Definition: TCanvas.h:48
double city_distance(Stsp_city c1, Stsp_city c2)
Definition: simanTSP.cxx:58
return c2
Definition: legend2.C:14
bool showGraphics
Definition: simanTSP.cxx:27
static const double x1[5]
double f(double x)
void Print(std::ostream &os, const OptionType &opt)
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
static RooMathCoreReg dummy
The TH1 histogram class.
Definition: TH1.h:80
double second_E
Definition: simanTSP.cxx:325
double distance_matrix[N_CITIES][N_CITIES]
Definition: simanTSP.cxx:55
#define name(a, b)
Definition: linkTestLib0.cpp:5
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1077
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
AxisAngle::Scalar Distance(const AxisAngle &r1, const R &r2)
Distance between two rotations.
Definition: AxisAngle.h:326
bool debug
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:45
const Int_t n
Definition: legend1.C:16
void print_distance_matrix()
Definition: simanTSP.cxx:85
structure holding the simulated annealing parameters
double best_E
Definition: simanTSP.cxx:325