|
|
/*
Ant Colony System for solving a TSP problem
By Andres Perez-Uribe
LSA2 - STI - EPFL
January, 2002
Email: Andres.Perez-Uribe AT a3.epfl.ch
Based on:
Dorigo and Gambardella, Ant Colony System: A Cooperative Learning
Approach to the Traveling Salesman Problem,
IEEE, Trans Evol. Comp. 1, 1997, 53-66
WARNING: This is not a 100% implementation of Dorigo et al. algorithm,
I have only used the basic idea. In particular, I have used a simpler
exploration-exploitation strategy, and the value of the initial
pheromone level is based on what I've called the 'naive-distance'.
Please, refer to Dorigo et al. for the original algorithm. I think mine is
simpler, but may not work as well as theirs (especially because of
random exploration instead of their biased exploration.
Original distance calculation:
[1275068416] mincost=30.8785039 13 2 3 4 5 11 6 12 7 10 8 9 0 1
Elapsed: 19590719 ms = 5.5 hours
Updated distance calculation based on FAQ at TSPLIB
Elapsed 23234093ms = 6.5 hours
Best route: mincost=3323 0 1 13 2 3 4 5 11 6 12 7 10 8 9
The 3323 matches the best cost reported at TSPLIB for Burma14.
*/
#include <limits.h>
#include <time.h>
#include <iostream>
#include <iomanip>
using namespace std;
#include "Common.h"
#include "AntSwarm.h"
#include "Burma.h"
#include "Ulysses16.h"
#include "Ulysses22.h"
//---------------------------------------------------------
template < class facet, int dim, int numants >
void RunSuite(AntSwarmType < facet, dim, numants > & swarm,
int maxiter, int numruns, int best, int limit1, int limit2)
{
int belowbest = 0;
int numbest = 0;
int numlimit1 = 0;
int numlimit2 = 0;
int others = 0;
time_t start = time(0);
for (int run = 0; run < numruns; run++)
{
swarm.InitClass();
swarm.Initialize();
swarm.Run(maxiter);
cout << " Best : Len: " << swarm.BestSolutionCost << " " << swarm.BestSolution << endl;
if (swarm.BestSolutionCost.Cost() < best)
belowbest++;
else if (swarm.BestSolutionCost.Cost() == best)
numbest++;
else if (swarm.BestSolutionCost.Cost() < limit1)
numlimit1++;
else if (swarm.BestSolutionCost.Cost() < limit2)
numlimit2++;
else
others ++;
}
time_t end = time(0);
cout << endl;
cout << "Number of Runs : " << numruns << endl;
cout << "Iterations per Run: " << maxiter << endl;
cout << "Number of Ants : " << numants << endl;
cout << endl;
if (belowbest != 0)
cout << "# runs that better than the optimum: " << belowbest << "\n";
cout << "% runs == " << best << " (optimum) : " << (numbest * 100) / numruns << "\n";
cout << "% runs < " << limit1 << " : " << (numlimit1 * 100) / numruns << "\n";
cout << "% runs < " << limit2 << " : " << (numlimit2 * 100) / numruns << "\n";
cout << "% remaining : " << (others * 100) / numruns << "\n";
int totaltime = (end - start);
cout << "Took: " << totaltime << " seconds elapsed time\n";
cout << "Took: " << totaltime / numruns << " seconds per run\n";
}
//---------------------------------------------------------
int main (int argc, char** argv)
{
// CommandLine(argc, argv);
//
// if (!doExhaustive)
// {
// assert(NumAnts <= NumDimensions);
// }
//
// if (doExhaustive)
// {
// ExhaustiveSearch(NumDimensions);
// return 0;
// }
int t = time(0);
cout << "t=" << t << "\n";
srand(t);
// srand(time(0));
// srand(1066448641);
const int NumRuns = 100;
//AntSwarmType<Burma, Burma::NumDim, Burma::NumDim - 2> swarm;
//RunSuite(swarm, 160, NumRuns, 3323, 3400, 3800);
//AntSwarmType<Ulysses16, Ulysses16::NumDim, Ulysses16::NumDim - 2> swarm;
//RunSuite(swarm, 1100, NumRuns, 6859, 7000, 7200);
AntSwarmType < Ulysses22, Ulysses22::NumDim, Ulysses22::NumDim - 2 > swarm;
RunSuite(swarm, 32000, NumRuns, 7013, 7100, 7200);
return 0;
}
|
|
|
#ifndef Burma_h
#define Burma_h
#include "TspFacetBase.h"
class Burma : public TspFacetBase < 14, 3323 >
{
public:
SolutionCost InitializeFacet()
{
InitClass();
double Xpos[NumDim];
double Ypos[NumDim];
//These are not coordinates in a cartesion plane
//These are in degrees, ie. degrees.minutes
Xpos[0] = 16.47;
Ypos[0] = 96.10;
Xpos[1] = 16.47;
Ypos[1] = 94.44;
Xpos[2] = 20.09;
Ypos[2] = 92.54;
Xpos[3] = 22.39;
Ypos[3] = 93.37;
Xpos[4] = 25.23;
Ypos[4] = 97.24;
Xpos[5] = 22.00;
Ypos[5] = 96.05;
Xpos[6] = 20.47;
Ypos[6] = 97.02;
Xpos[7] = 17.20;
Ypos[7] = 96.29;
Xpos[8] = 16.30;
Ypos[8] = 97.38;
Xpos[9] = 14.05;
Ypos[9] = 98.12;
Xpos[10] = 16.53;
Ypos[10] = 97.38;
Xpos[11] = 21.52;
Ypos[11] = 95.59;
Xpos[12] = 19.41;
Ypos[12] = 97.13;
Xpos[13] = 20.09;
Ypos[13] = 94.55;
mDistances.Initialize(Xpos, Ypos);
return ComputeNaiveDistance();
}
} ;
#endif
|
|
|
#ifndef Ulysses16_h
#define Ulysses16_h
#include "TspFacetBase.h"
class Ulysses16 : public TspFacetBase < 16, 6859 >
{
public:
SolutionCost InitializeFacet()
{
InitClass();
double Xpos[NumDim];
double Ypos[NumDim];
//These are not coordinates in a cartesion plane
//These are in degrees, ie. degrees.minutes
Xpos[0] = 38.24; // 1 38.24 20.42
Ypos[0] = 20.42;
Xpos[1] = 39.57; // 2 39.57 26.15
Ypos[1] = 26.15;
Xpos[2] = 40.56; // 3 40.56 25.32
Ypos[2] = 25.32;
Xpos[3] = 36.26; // 4 36.26 23.12
Ypos[3] = 23.12;
Xpos[4] = 33.48; // 5 33.48 10.54
Ypos[4] = 10.54;
Xpos[5] = 37.56; // 6 37.56 12.19
Ypos[5] = 12.19;
Xpos[6] = 38.42; // 7 38.42 13.11
Ypos[6] = 13.11;
Xpos[7] = 37.52; // 8 37.52 20.44
Ypos[7] = 20.44;
Xpos[8] = 41.23; // 9 41.23 9.10
Ypos[8] = 9.10;
Xpos[9] = 41.17; // 10 41.17 13.05
Ypos[9] = 13.05;
Xpos[10] = 36.08; // 11 36.08 -5.21
Ypos[10] = -5.21;
Xpos[11] = 38.47; // 12 38.47 15.13
Ypos[11] = 15.13;
Xpos[12] = 38.15; // 13 38.15 15.35
Ypos[12] = 15.35;
Xpos[13] = 37.51; // 14 37.51 15.17
Ypos[13] = 15.17;
Xpos[14] = 35.49; // 15 35.49 14.32
Ypos[14] = 14.32;
Xpos[15] = 39.36; // 16 39.36 19.56
Ypos[15] = 19.56;
mDistances.Initialize(Xpos, Ypos);
return ComputeNaiveDistance();
}
} ;
#endif
|
|
|
#ifndef Ulysses22_h
#define Ulysses22_h
#include "TspFacetBase.h"
class Ulysses22 : public TspFacetBase < 22, 7013 >
{
public:
SolutionCost InitializeFacet()
{
InitClass();
double Xpos[NumDim];
double Ypos[NumDim];
//These are not coordinates in a cartesion plane
//These are in degrees, ie. degrees.minutes
Xpos[0] = 38.24; // 1 38.24 20.42
Ypos[0] = 20.42;
Xpos[1] = 39.57; // 2 39.57 26.15
Ypos[1] = 26.15;
Xpos[2] = 40.56; // 3 40.56 25.32
Ypos[2] = 25.32;
Xpos[3] = 36.26; // 4 36.26 23.12
Ypos[3] = 23.12;
Xpos[4] = 33.48; // 5 33.48 10.54
Ypos[4] = 10.54;
Xpos[5] = 37.56; // 6 37.56 12.19
Ypos[5] = 12.19;
Xpos[6] = 38.42; // 7 38.42 13.11
Ypos[6] = 13.11;
Xpos[7] = 37.52; // 8 37.52 20.44
Ypos[7] = 20.44;
Xpos[8] = 41.23; // 9 41.23 9.10
Ypos[8] = 9.10;
Xpos[9] = 41.17; // 10 41.17 13.05
Ypos[9] = 13.05;
Xpos[10] = 36.08; // 11 36.08 -5.21
Ypos[10] = -5.21;
Xpos[11] = 38.47; // 12 38.47 15.13
Ypos[11] = 15.13;
Xpos[12] = 38.15; // 13 38.15 15.35
Ypos[12] = 15.35;
Xpos[13] = 37.51; // 14 37.51 15.17
Ypos[13] = 15.17;
Xpos[14] = 35.49; // 15 35.49 14.32
Ypos[14] = 14.32;
Xpos[15] = 39.36; // 16 39.36 19.56
Ypos[15] = 19.56;
Xpos[16] = 38.09; // 17 38.09 24.36
Ypos[16] = 24.36;
Xpos[17] = 36.09; // 18 36.09 23.00
Ypos[17] = 23.00;
Xpos[18] = 40.44; // 19 40.44 13.57
Ypos[18] = 13.57;
Xpos[19] = 40.33; // 20 40.33 14.15
Ypos[19] = 14.15;
Xpos[20] = 40.37; // 21 40.37 14.23
Ypos[20] = 14.23;
Xpos[21] = 37.57; // 22 37.57 22.56
Ypos[21] = 22.56;
mDistances.Initialize(Xpos, Ypos);
return ComputeNaiveDistance();
}
} ;
#endif
|
|
|
#ifndef Ant_h
#define Ant_h
#include <string>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <assert.h>
using namespace std;
#include "Common.h"
#include "Solution.h"
#include "SolutionCost.h"
#include "Pheromone.h"
#include "DistanceCalculator.h"
//- NumUnvi
template < int NumDim, class PheromoneT, class SwarmFacetT > //, double q0 = 0.9>
class AntFacet
{
public:
const Solution& GetSolution()
{
return mSolution;
}
int NumUnvisitedCities()
{
int count = 0;
for (int i = 0; i < NumDim; ++i)
if (!mVisitedCities[i])
count++;
return count;
}
bool HasVisited(int i)
{
return mVisitedCities[i];
}
void InitializeSolution(void)
{
mSolution.Clear(NumDim);
InitializeVisitedCities();
}
void BuildSolutionAlongDimension(int d, PheromoneT& pheromone)
{
if (d == 0 || NumUnvisitedCities() == 0)
AddNextCity(d, mStartCity); /* ants start and end at the initial city */
else
AddNextCity(d, ChooseNextCity(mSolution.At(d - 1), pheromone));
}
protected:
void InitializeAntFacet(int k, SwarmFacetT& facet)
{
mFacet = &facet;
mStartCity = mFacet->AntIsInCity[k];
}
void AddNextCity(int d, int city)
{
mSolution.SetAt(d, city);
mVisitedCities[city] = true;
}
private:
int ChooseNextCity(int d, PheromoneT& tao)
{
/* exploration rate - the percentage of time we're going to explore vs exploit*/
const double q0 = 0.9;
/* exploration-exploitation */
double prob_explore = rand () / RAND_MAX;
return (prob_explore <= q0) ? CityWithMaxWeightedPheromone(d, tao) : RandomUnvisitedCity();
}
int CityWithMaxWeightedPheromone(int d, PheromoneT& tao)
{
double WeightedPheromone[NumDim];
int bestcity = 0;
for (int i = 0; i < NumDim; i++)
{
if (HasVisited(i))
WeightedPheromone[i] = 0.0;
else
{
double dist = mFacet->ComputeDistance(d, i);
WeightedPheromone[i] = tao.At(d, i) / (dist * dist);
if (WeightedPheromone[i] > WeightedPheromone[bestcity])
bestcity = i;
}
}
return bestcity;
}
int RandomUnvisitedCity()
{
//of the remaining unvisited cities, choose a random one
int RndCity;
int unvisited = NumUnvisitedCities();
if (unvisited == 1)
RndCity = 1;
else
RndCity = 1 + rand() % (unvisited - 1);
int j = 0;
int city;
for (city = 0; city < NumDim && j < RndCity; city++)
{
if (!HasVisited(city))
j++;
}
return city;
}
void InitializeVisitedCities()
{
for (int j = 0; j < NumDim; j++)
mVisitedCities[j] = false;
}
protected:
MySolution mSolution;
SwarmFacetT* mFacet;
int mStartCity;
bool mVisitedCities[NumDim];
} ;
template < int NumDim , class PheromoneT, class SwarmFacetT >
class Ant : public AntFacet < NumDim, PheromoneT, SwarmFacetT >
{
public:
void Initialize(int k, SwarmFacetT& facet)
{
char buf[20];
sprintf(buf, "%d", k);
mName = buf;
mIndex = k;
this->InitializeSolution();
this->InitializeAntFacet(k, facet);
}
SolutionCost GetCost()
{
return this->mSolutionCost;
}
bool IsBetterThan(const SolutionCost& cost)
{
return GetCost() < cost;
}
void ComputeSolution()
{
mSolutionCost = this->mFacet->ComputeCostOf(this->mSolution);
}
void PrintOn(ostream& os)
{
os << "Ant: [" << Name() << "] Len: " << this->mSolutionCost << " " << this->mSolution << endl;
}
const string& Name()
{
return mName;
}
private:
SolutionCost mSolutionCost;
string mName;
int mIndex;
} ;
#endif
|
|
|
#ifndef AntSwarm_h
#define AntSwarm_h
#include "Common.h"
#include "Ant.h"
#include "Pheromone.h"
#include "Solution.h"
#include "SolutionCost.h"
//
//SwarmFacet:
// the AntSwarm template contains all of the generic code for the Ant Swarm algorithm,
// SwarmFacet contains all the specific code.
// - must have a methods:
// - SolutionCost InitializeFacet()
// - void PrintIteration(int iter, const SolutionCost& cost, const Solution& solution)
// - void ComputeSolution(AntType& ant)
template < int NumAnts, class AntT, class PheromoneT, class SwarmFacetT >
class Ants
{
public:
void Initialize(SwarmFacetT& facet)
{
InitializeAnts(facet);
}
void InitializeIteration()
{
InitializeSolutions();
}
void UpdatePheromoneLevelLocal(PheromoneT& pheromone)
{
for (int k = 0; k < NumAnts; k++)
pheromone.UpdatePheromoneLevelLocal(mAnts[k].GetSolution(), mAnts[k].GetCost());
}
void BuildSolutionAlongDimension(int d, PheromoneT& pheromone)
{
for (int k = 0; k < NumAnts; k++)
mAnts[k].BuildSolutionAlongDimension(d, pheromone);
}
void FindBestSolution (SwarmFacetT& facet, SolutionCost& bestcost, Solution& bestsolution)
{
ComputeSolutions(facet);
FindBest(bestcost, bestsolution);
}
private:
void InitializeAnts(SwarmFacetT& facet)
{
for (int k = 0; k < NumAnts; k++)
mAnts[k].Initialize(k, facet);
}
void InitializeSolutions()
{
for (int k = 0; k < NumAnts; k++)
mAnts[k].InitializeSolution();
}
void ComputeSolutions(SwarmFacetT& facet)
{
for (int k = 0; k < NumAnts; k++)
mAnts[k].ComputeSolution();
}
void FindBest(SolutionCost& bestcost, Solution& bestsolution)
{
SolutionCost cost = bestcost;
for (int k = 0; k < NumAnts; k++)
{
if (mAnts[k].IsBetterThan(cost))
{
bestsolution = mAnts[k].GetSolution();
cost = mAnts[k].GetCost();
}
}
bestcost = cost;
}
private:
AntT mAnts[NumAnts];
} ;
template < int NumDim, class AntsT, class PheromoneT, class SwarmFacetT >
class AntSwarm : public SwarmFacetT
{
public:
Solution BestSolution; //The best solution found
SolutionCost BestSolutionCost; //The cost of the best solution
void Initialize()
{
mPheromone.Initialize(this->InitializeFacet());
mAnts.Initialize(*this);
BestSolutionCost = SolutionCost::WorstSolution();
}
void Run(int maxiterations)
{
SolutionCost bestoverallcost = SolutionCost::WorstSolution();
Solution bestoverallsolution;
mStagnantCount = 0;
for (int iter = 0; iter < maxiterations; iter++)
{
InitializeIteration();
BuildSolutions();
FindBestSolution();
this->PrintIteration(iter, BestSolutionCost, BestSolution);
UpdatePheromoneLevels();
if (mStagnantCount > 10)
{
//cout << "Stagnant!!\n";
if (BestSolutionCost < bestoverallcost)
{
bestoverallcost = BestSolutionCost;
bestoverallsolution = BestSolution;
}
BestSolutionCost = SolutionCost::WorstSolution();
mPheromone.Clear();
this->InitClass();
mAnts.Initialize(*this);
mStagnantCount = 0;
}
}
if (bestoverallcost < BestSolutionCost )
{
BestSolutionCost = bestoverallcost;
BestSolution = bestoverallsolution;
}
}
private:
void InitializeIteration()
{
mAnts.InitializeIteration();
}
void BuildSolutions()
{
for (int dimension = 0; dimension < NumDim + 1 ; ++dimension)
mAnts.BuildSolutionAlongDimension(dimension, mPheromone);
}
void FindBestSolution()
{
SolutionCost lastbest = BestSolutionCost;
mAnts.FindBestSolution(*this, BestSolutionCost, BestSolution);
if (lastbest == BestSolutionCost)
mStagnantCount++;
else
mStagnantCount = 0;
}
void UpdatePheromoneLevels()
{
mAnts.UpdatePheromoneLevelLocal(mPheromone);
mPheromone.Decay();
mPheromone.UpdatePheromoneLevelGlobal(BestSolution, BestSolutionCost);
}
int mStagnantCount;
AntsT mAnts;
PheromoneT mPheromone;
};
template < class facet, int dim, int numants >
struct AntSwarmType : public AntSwarm < dim,
Ants < numants,
Ant < dim, Pheromone < dim > , facet > ,
Pheromone < dim > ,
facet > ,
Pheromone < dim > ,
facet >
{}
;
#endif
|
|
|
#ifndef common_h
#define common_h
const bool gVerbose = false;
#endif
|
|
|
#ifndef DISTANCE_CALCULATOR_H
#define DISTANCE_CALCULATOR_H
#include <math.h>
class DistanceCalculator
{
public:
virtual double operator () (int i, int j) const = 0;
} ;
template <int numdim>
class DCalc : public DistanceCalculator
{
public:
double DCalc::operator () (int i, int j) const
{
return mDistanceVector[i * numdim + j];
}
void DCalc::Initialize(double* X, double* Y)
{
for (int i = 0; i < numdim; i++)
{
mDistanceVector[i * numdim + i] = 0;
for (int j = i + 1; j < numdim; j++)
{
int d = CalcDistanceBetween(X[i], Y[i], X[j], Y[j]);
mDistanceVector[i * numdim + j] = d;
mDistanceVector[j * numdim + i] = d;
}
}
}
private:
static const double PI = 3.141592;
inline double DecDegsToRadians(double x)
{
double deg = (int) x; //split off the degrees
double min = x - deg; //split off the minutes
// (equal to rad = PI * (deg + (100.0 * min) / 60.0) / 180.0; )
return (PI * (deg + 5.0 * min / 3.0)) / 180.0;
}
int DCalc::CalcDistanceBetween(double Xi, double Yi, double Xj, double Yj)
{
//Q: I get wrong distances for problems of type GEO.
//A: There has been some confusion of how to compute the distances. I use the following code.
//
//Convert coordinate input to longitude and latitude in radian:
double latitudei = DecDegsToRadians(Xi);
double longitudei = DecDegsToRadians(Yi);
double latitudej = DecDegsToRadians(Xj);
double longitudej = DecDegsToRadians(Yj);
//Compute the geographical distance:
const double RRR = 6378.388; //radius of the earth in kilometers
double q1 = cos(longitudei - longitudej);
double q2 = cos(latitudei - latitudej);
double q3 = cos(latitudei + latitudej);
return (int) ( RRR * acos( 0.5 * ((1.0 + q1)*q2 - (1.0 - q1)*q3) ) + 1.0);
}
int mDistanceVector[numdim * (numdim + 1)];
} ;
#endif
|
|
|
#ifndef Pheromone_h
#define Pheromone_h
#include <iostream>
#include <vector>
using std::cout;
using std::vector;
#include "Common.h"
#include "SolutionCost.h"
#include "Solution.h"
const double RHO = 0.1;
template < int NumDim >
class Pheromone
{
public:
void Initialize(const SolutionCost& guess)
{
Tao0 = GetInitial(guess);
if (gVerbose)
cout << "Tao0 = " << setprecision(9) << Tao0 << " Naivedist=" << guess << "\n";
Clear();
}
void Clear()
{
for (int i = 0; i < NumDim; i++)
{
for (int j = 0; j < NumDim; j++)
mTao[i][j] = Tao0;
}
}
double At(int i, int j)
{
return mTao[i][j];
}
void UpdatePheromoneLevelGlobal(const Solution& best, const SolutionCost& bestcost)
{
for (int d = 0; d < NumDim; d++)
{
int i = best.At(d);
int j = best.At(d + 1);
mTao[i][j] = (1 - RHO) * mTao[i][j] + RHO * (1.0 / bestcost.Cost());
}
}
void UpdatePheromoneLevelLocal(const Solution& soln, const SolutionCost& cost)
{
for (int d = 0; d < NumDim; d++)
{
int i = soln.At(d);
int j = soln.At(d + 1);
mTao[i][j] = (1 - RHO) * mTao[i][j] + RHO * UpdateAmount(mTao[i][j]);
}
}
void Decay()
{
for (int i = 0; i < NumDim; i++)
for (int j = 0; j < NumDim; j++)
mTao[i][j] -= DecayAmount(mTao[i][j]);
}
double DecayAmount(double current)
{
return current / 10.0;
}
double UpdateAmount(double current)
{
//return Tao0 / 10.0;
return current / 9.0;
}
double GetInitial(const SolutionCost& guess)
{
return 1.0 / (guess.Cost() * NumDim);
}
private:
double Tao0;
double mTao[NumDim][NumDim];
} ;
//#define RHO 0.8
#endif
|
|
|
#include <algorithm>
using namespace std;
#include "common.h"
#include "Solution.h"
#include "DistanceCalculator.h"
Solution::~Solution()
{}
void Solution::SetAt(int i, int val)
{
mSolution[i] = val;
}
int Solution::At(int i) const
{
return mSolution.at(i);
}
const Solution& Solution::operator = (const Solution& other)
{
mSolution = other.mSolution;
return *this;
}
void Solution::PrintOn(ostream& os) const
{
for (vector < int > ::const_iterator it = mSolution.begin(); it != mSolution.end(); ++it)
{
if (it == mSolution.begin())
os << setw(2) << (*it);
else
os << "-> " << setw(2) << (*it);
}
}
void MySolution::Clear(int numdimensions)
{
cNumDimensions = numdimensions;
mSolution.resize(cNumDimensions + 1);
mSolution.assign(cNumDimensions + 1, 0);
}
SolutionCost MySolution::Length(const DistanceCalculator& dc) const
{
double len = 0.0;
for (int d = 0; d < cNumDimensions; ++d)
len += dc(mSolution[d], mSolution[d + 1]);
len += dc(mSolution[cNumDimensions], mSolution[0]);
return SolutionCost(len);
}
void MySolution::Fill()
{
mSolution.resize(cNumDimensions + 1);
for (int i = 0; i < cNumDimensions; ++i)
mSolution[i] = i;
}
bool MySolution::Next()
{
bool res = next_permutation(mSolution.begin(), --mSolution.end());
mSolution[cNumDimensions + 1] = mSolution[0];
return res;
}
|
|
|
#ifndef Solution_h
#define Solution_h
#include <iomanip>
#include <ostream>
#include <vector>
using std::ostream;
using std::setw;
using std::vector;
#include "SolutionCost.h"
#include "DistanceCalculator.h"
class Solution
{
public:
int At(int i) const;
void SetAt(int i, int val);
const Solution& operator = (const Solution& other);
virtual void PrintOn(ostream& os) const;
virtual ~Solution();
protected:
vector < int > mSolution;
friend ostream& operator << (ostream& os, const Solution& solution);
} ;
inline ostream& operator << (ostream& os, const Solution& solution)
{
solution.PrintOn(os);
return os;
}
class MySolution : public Solution
{
public:
void Clear(int numdimensions);
SolutionCost Length(const DistanceCalculator& dc) const;
void Fill();
bool Next();
int cNumDimensions;
} ;
#endif
|
|
|
#ifndef SolutionCost_h
#define SolutionCost_h
#include "SolutionCostT.h"
#include <ostream>
#include <iomanip>
using std::ostream;
using std::setprecision;
typedef SolutionCostT < double > SolutionCost;
inline ostream& operator << (ostream& os, const SolutionCost& cost)
{
os << setprecision(9) << cost.Cost();
return os;
}
#endif
|
|
|
#ifndef SolutionCostT_h
#define SolutionCostT_h
template < typename CostT >
class SolutionCostT
{
public:
static const SolutionCostT WorstSolution()
{
return SolutionCostT(INT_MAX);
}
SolutionCostT()
: mCost(0.0)
{
//blah
}
explicit SolutionCostT(CostT c)
{
mCost = c;
}
SolutionCostT operator = (const SolutionCostT& other)
{
mCost = other.mCost;
return *this;
}
bool operator < (const SolutionCostT& other)
{
return mCost < other.mCost;
}
bool operator == (const SolutionCostT& other)
{
return mCost == other.mCost;
}
CostT Cost() const
{
return mCost;
}
private:
CostT mCost;
} ;
#endif
|
|
|
#ifndef TspFacetBase_h
#define TspFacetBase_h
#include <vector>
using std::vector;
#include "SolutionCost.h"
#include "Solution.h"
template < int numdim, int bestknown >
class TspFacetBase
{
public:
static const int NumDim = numdim;
static const int BestKnown = bestknown;
vector < int > AntIsInCity;
void InitClass()
{
static bool firsttime = true;
if (firsttime)
{
firsttime = false;
AntIsInCity.resize(NumDim);
for (int j = 0; j < NumDim; j++)
AntIsInCity[j] = j;
}
random_shuffle(AntIsInCity.begin(), AntIsInCity.end());
}
void PrintIteration(int iter, const SolutionCost& cost, const Solution& solution)
{
//if (gVerbose)
//cout << " Iter[" << iter << "] Len: " << cost << " " << solution << endl;
}
SolutionCost ComputeCostOf(Solution& solution)
{
return ((MySolution&)solution).Length(mDistances);
}
double ComputeDistance(int from, int to)
{
return mDistances(from, to);
}
SolutionCost ComputeNaiveDistance()
{
double SumNaiveDist = 0.0;
for (int i = 0; i < NumDim; i++)
SumNaiveDist += mDistances(i, (i + 1) % NumDim);
return SolutionCost(SumNaiveDist);
}
protected:
DCalc<numdim> mDistances;
} ;
#endif
|
|
|
#include <iostream>
using namespace std;
#include "Utilities.h"
#include "Solution.h"
#include "Burma.h"
void ExhaustiveSearch(int numdimensions)
{
cout << "Exhaustive Search on " << numdimensions << " dimensions...\n";
clock_t start = clock();
MySolution cities;
cities.Clear(numdimensions);
Burma burma;
burma.InitializeFacet();
cout << setprecision(9);
cities.Fill();
cities.SetAt(numdimensions, cities.At(0));
SolutionCost mincost = burma.ComputeCostOf(cities);
long j = 0;
cout << "[" << j << "] mincost=" << mincost << " ";
cities.PrintOn(cout);
cout << endl;
Solution bestcities = cities;
SolutionCost c;
while (cities.Next() ) // && j++ < 1000)
{
c = burma.ComputeCostOf(cities);
if (c < mincost)
{
mincost = c;
bestcities = cities;
cout << "[" << j << "] *** mincost=" << mincost << " ";
bestcities.PrintOn(cout);
cout << endl;
}
if ((j & 0x000FFFFF) == 0)
{
cout << "[" << j << "] mincost=" << mincost << " c=" << c << " ";
cities.PrintOn(cout);
cout << endl;
}
j++;
}
cout << "Best route: mincost=" << mincost << " ";
bestcities.PrintOn(cout);
//13 & 14 dim times are with a release build, the rest are debug builds
// time #dim best Route
// 5.5hrs 14 3323
// ? 13 3158
// 40.8m 12 3150 0-> 1-> 2-> 3-> 4-> 5->11-> 6-> 7->10-> 8-> 9-> 0
//191.9s 11 3136 0-> 1-> 2-> 3-> 4-> 5-> 6-> 7->10-> 8-> 9-> 0
// 16.4s 10 3114 0-> 1-> 2-> 3-> 4-> 5-> 6-> 7-> 8-> 9-> 0
// 1.5s 9 2626 0-> 1-> 2-> 3-> 4-> 5-> 6-> 7-> 8-> 0
// 0.2s 8 2382 0-> 1-> 2-> 3-> 4-> 5-> 6-> 7-> 0
// 15ms 7 2378 0-> 1-> 2-> 3-> 4-> 5-> 6-> 0
// 15ms 6 2336 0-> 1-> 2-> 3-> 4-> 5-> 0
// 15ms 5 2321 0-> 1-> 2-> 3-> 4-> 0
// 15ms 4 1570 0-> 1-> 2-> 3-> 0
// 15ms 3 1085 0-> 1-> 2-> 0
// 15ms 2 306 0-> 1-> 0
clock_t end = clock();
cout << "\nElapsed: " << end - start << endl;
}
//static bool doExhaustive = false;
//static int MAXITER = 160;//500;//20;
//static int NumDimensions = 14;
//static int NumAnts = 2 * NumDimensions / 3;
//
////---------------------------------------------------------
//void usage()
//{
// cout << "usage: pso [-x [numdim]] [[maxiter] [numdim] [numants]]\n"
// " numdim : number of dimensions (default=max=14)\n"
// " maxiter: number of iterations (default = 100)\n"
// " numants: number of ants to use (<= numdim)\n"
// " x : do exhaustive search\n";
//
// exit(0);
//}
//
////---------------------------------------------------------
//void CommandLine(int argc, char** argv)
//{
// //pso -x [maxiter] [numdim] [numants]
// bool doExhaustive = false;
// int posn = 0;
// for (int i = 1; i < argc; ++i)
// {
// if (strcmp(argv[i], "-x") == 0)
// doExhaustive = true;
// else if (posn == 0)
// {
// if (doExhaustive)
// NumDimensions = atoi(argv[i]);
// else
// MAXITER = atoi(argv[i]);
// posn++;
// }
// else if (posn == 1)
// {
// if (doExhaustive)
// usage();
// NumDimensions = atoi(argv[i]);
// NumAnts = 2 * NumDimensions / 3;
// posn++;
// }
// else if (posn == 2)
// {
// NumAnts = atoi(argv[i]);
// posn++;
// }
// else
// usage();
// }
//}
|
|
|
extern void ExhaustiveSearch(int numdimensions);
|