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;
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";
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;
292 xmin[i] = -
cities[f.Route(i)].longitude;
293 ymin[i] =
cities[f.Route(i)].lat;
331 unsigned int *
r = f.Route();
334 while (std::next_permutation(r+offset,r+
N_CITIES) ) {
336 double E = f.Energy();
365 const unsigned int *
r = f.Route();
387 const unsigned int *
r = f.Route();
409 std::vector<double>
dummy;
411 MySimAnFunc
f(dummy);
415 const unsigned int *
r = f.Route();
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: ");
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.
unsigned int r3[N_CITIES]
double k
parameters for the Boltzman distribution
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
unsigned long RndmInt(unsigned long 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.
unsigned int r1[N_CITIES]
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)
Set the marker style.
double city_distance(Stsp_city c1, Stsp_city c2)
static const double x1[5]
void Print(std::ostream &os, const OptionType &opt)
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.
double f2(const double *x)
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.
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
void print_distance_matrix()
unsigned int r2[N_CITIES]
structure holding the simulated annealing parameters