36 double lat, longitude;
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}};
53 #define N_CITIES (sizeof(cities)/sizeof(Stsp_city))
60 const double earth_radius = 6375.000;
67 double x1 = cla1*clo1;
68 double x2 = cla2*clo2;
70 double y1 = cla1*slo1;
71 double y2 = cla2*slo2;
76 double dot_product = x1*x2 + y1*y2 + z1*z2;
81 return angle*earth_radius;
92 printf(
"%15.8f ", distance_matrix[i][j]);
106 MySimAnFunc( std::vector<double> & allDist)
108 calculate_distance_matrix();
110 for (
unsigned int i = 0; i <
N_CITIES; ++i)
118 virtual ~MySimAnFunc() {}
120 unsigned int Route(
unsigned int i)
const {
return fRoute[i]; }
122 const unsigned int * Route()
const {
return fRoute; }
123 unsigned int * Route() {
return fRoute; }
125 virtual MySimAnFunc * Clone()
const {
return new MySimAnFunc(*
this); }
127 std::vector<double> & AllDist() {
return *fDist; }
129 virtual double Energy()
const {
134 for (
unsigned int i = 0; i <
N_CITIES; ++i) {
137 enrg += fDistanceMatrix( fRoute[i] , fRoute[ (i + 1) % N_CITIES ] );
145 const MySimAnFunc * f2 =
dynamic_cast<const MySimAnFunc *
> (&
f);
149 for (
unsigned int i = 0; i <
N_CITIES; ++i) {
150 d += (( fRoute[i] == f2->Route(i) ) ? 0 : 1 );
167 fRoute[
x1] = fRoute[
x2];
174 virtual void Print() {
176 for (
unsigned i = 0; i <
N_CITIES; ++i) {
177 printf(
" %d ", fRoute[i]);
180 fDist->push_back(Energy());
184 virtual MySimAnFunc & FastCopy(
const GSLSimAnFunc & f) {
185 const MySimAnFunc * rhs =
dynamic_cast<const MySimAnFunc *
>(&
f);
187 std::copy(rhs->fRoute, rhs->fRoute + N_CITIES, fRoute);
191 double Distance(
int i,
int j)
const {
return fDistanceMatrix(i,j); }
195 void SetRoute(
unsigned int * r) { std::copy(r,r+N_CITIES,fRoute); }
199 void calculate_distance_matrix();
204 Matrix fDistanceMatrix;
206 std::vector<double> * fDist;
212 void MySimAnFunc::calculate_distance_matrix()
218 for (j = 0; j <= i; ++j) {
224 fDistanceMatrix(i,j) =
dist;
229 void MySimAnFunc::PrintRoute() {
232 for (
unsigned int i = 0; i <
N_CITIES; ++i) {
233 std::cout << std::setw(20) << cities[Route(i)].name <<
" \t " << Route(i);
235 int k = Route( (i+ 1) % N_CITIES );
237 std::cout <<
"\tdistance [" << j <<
" , " << k <<
" ]\t= " <<
Distance(j,k) <<
"\tTotal Distance\t = " << dtot << std::endl;
239 std::cout <<
"Total Route energy is " << dtot << std::endl;
246 std::vector<double> allDist;
247 allDist.reserve(5000);
250 MySimAnFunc
f(allDist);
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;
272 std::cout <<
"Best Route is \n";
276 double x0[N_CITIES+1];
277 double y0[N_CITIES+1];
278 double xmin[N_CITIES+1];
279 double ymin[N_CITIES+1];
283 TH1 *
h1 =
new TH1D(
"h1",
"Distance",niter+1,0,niter+1);
284 for (
unsigned int i = 1; i <=niter; ++i) {
288 for (
unsigned int i = 0; i <
N_CITIES; ++i) {
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;
324 unsigned int nfiter = 0;
331 unsigned int * r = f.Route();
334 while (std::next_permutation(r+offset,r+N_CITIES) ) {
336 double E = f.Energy();
341 std::copy(r2,r2+N_CITIES,r3);
343 std::copy(r1,r1+N_CITIES,r2);
345 std::copy(r,r+N_CITIES,r1);
346 }
else if (E < second_E) {
348 std::copy(r2,r2+N_CITIES,r3);
350 std::copy(r,r+N_CITIES,r2);
351 }
else if (E < third_E) {
353 std::copy(r,r+N_CITIES,r3);
363 if (n == (N_CITIES-1)) {
365 const unsigned int * r = f.Route();
369 std::copy(r2,r2+N_CITIES,r3);
371 std::copy(r1,r1+N_CITIES,r2);
373 std::copy(r,r+N_CITIES,r1);
374 }
else if (
E < second_E) {
376 std::copy(r2,r2+N_CITIES,r3);
378 std::copy(r,r+N_CITIES,r2);
379 }
else if (
E < third_E) {
381 std::copy(r,r+N_CITIES,r3);
387 const unsigned int * r = f.Route();
388 std::copy(r,r+N_CITIES,newr);
409 std::vector<double>
dummy;
411 MySimAnFunc
f(dummy);
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);
426 std::cout <<
"number of calls " << nfiter << std::endl;
429 printf(
"\n# exhaustive best route: ");
430 std::cout << best_E << std::endl;
431 f.SetRoute(r1); f.PrintRoute();
434 printf(
"\n# exhaustive second_best route: ");
435 std::cout << second_E << std::endl;
436 f.SetRoute(r2); f.PrintRoute();
438 printf(
"\n# exhaustive third_best route: ");
439 std::cout << third_E << std::endl;
440 f.SetRoute(r3); f.PrintRoute();
448 int main(
int argc,
char **argv)
457 for (
Int_t i=1 ; i<argc ; i++) {
458 std::string arg = argv[i] ;
467 cerr <<
"Usage: " << argv[0] <<
" [-g] [-v]\n";
469 cerr <<
" -g : graphics mode\n";
470 cerr <<
" -v : verbose mode";
GSLSimAnnealing class for performing a simulated annealing search of a multidimensional function...
double dist(Rotation3D const &r1, Rotation3D const &r2)
GSLSimAnParams & Params()
GSLSimAnFunc class description.
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
void Step(const gsl_rng *r, void *xp, double step_size)
struct s_tsp_city Stsp_city
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
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.
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 ...
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)
void simanTSP(bool debug=true)
void do_all_perms(MySimAnFunc &f, int offset)
virtual void Draw(Option_t *option="")
Draw this histogram with options.
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...
virtual void SetMarkerStyle(Style_t mstyle=1)
double city_distance(Stsp_city c1, Stsp_city c2)
static const double x1[5]
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
double distance_matrix[N_CITIES][N_CITIES]
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.
A Graph is a graphics object made of two arrays X and Y with npoints each.
AxisAngle::Scalar Distance(const AxisAngle &r1, const R &r2)
Distance between two rotations.
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
void print_distance_matrix()
structure holding the simulated annealing parameters