Dość dawno to robiłem komuś, chyba na zaliczenie... no to jakoś tak wygląda: :)
Tam są chyba ze dwie operacje tylko, które wykonujemy losowo... no i ostatecznie mamy trasę, która jest optymalna, albo prawie.
Tak jest w przypadku żywego komiwojażera - on na czuja znajduje niemal optymalne rozwiązania... taka działa natura.
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "ag.h"
#pragma inline
TFPoint *pts; // tablica punktow -> miast
int *ti; // tablica pomocnicza
int *ruletka; // fit/sumfit = k[i]/N, k - krotnosc i-tego w tablicy N-elem.
int n_pts = 30;
int pop_siz = 200;
int pcross = 700, pmutation = 330;
float am=1, bm=0; // f = am*d + bm, wsp. do przeksztalcenia odleglosci
// na funkcję przystosowania
float max_tr;
float anilingt, bestnt;
inline void swap(int &a, int &b)
{ int x = a; a = b; b = x; }
int flip(int prob)
{ return random(1000) < prob; }
double dist(TFPoint *a, TFPoint *b)
{
return sqrt(sqrf(b->x-a->x)+sqrf(b->y-a->y));
}
//#define dist(a, b) hypot(b->x-a->x, b->y-a->y)
double droga(Node t[])
{
TFPoint *prv = pts + t[n_pts-1], // n -> 1 -> 2 -> 3 ... -> n
*nxt;
Node *e = t + n_pts;
double d = 0;
do{
nxt = pts + *t++;
d += dist(prv, nxt);
prv = nxt;
}while( t != e );
return d;
}
inline double xfit(float x)
{
// return 1/x; // bez skalowania = marnie
float f;
return (f = x*am + bm) >= 0 ? f : 0;
}
float Trasa::decode()
{
return fit = xfit(x = droga(trasa));
}
TKomi::TKomi()
{
pts = new TFPoint[MAX_PTS+8];
ti = new int[2*MAX_PTS+2 + P_P+2];
ruletka = ti + 2*MAX_PTS+2;
srand(GetTickCount());
reset();
}
TKomi::~TKomi()
{
delete []ti;
delete []pts;
}
// ===aniling===========================
float revcst(int iorder[], int ncity, int n[])
{
float de;
TFPoint *p[7];
n[3] = 1 + ((n[1]+ncity-2) % ncity); //Find the city before n[1] ..
n[4] = 1 + (n[2] % ncity); //.. and the city after n[2].
for(int j = 1; j <= 4; j++)
p[j] = pts+iorder[n[j]]; // Find coordinates for the four cities involved.
// Calculate cost of disconnecting the segment
// at both ends and reconnecting in the opposite order.
de = -dist(p[1], p[3]);
de -= dist(p[2], p[4]);
de += dist(p[1], p[4]);
de += dist(p[2], p[3]);
return de;
}
void reverse(int iorder[], int ncity, int n[])
{
int nn = (1+((n[2]-n[1]+ncity) % ncity))/2; // This many cities must be swapped to
for(int j = 1; j <= nn; j++) // effect the reversal.
{
int k = 1 + ((n[1]+j-2) % ncity); // Start at the ends of the segment and
int l = 1 + ((n[2]-j+ncity) % ncity); // swap pairs of cities, moving toward the center.
int itmp = iorder[k];
iorder[k] = iorder[l];
iorder[l] = itmp;
}
}
float trncst(int iorder[], int ncity, int n[])
{
float de;
TFPoint *p[7];
n[4] = 1 + (n[3] % ncity); //Find the city following n[3]..
n[5] = 1 + ((n[1]+ncity-2) % ncity); //..and the one preceding n[1]..
n[6] = 1 + (n[2] % ncity); //..and the one following n[2].
for(int j = 1; j <= 6; j++)
p[j] = pts + iorder[n[j]]; // Determine coordinates for the six cities involved.
de = -dist(p[2],p[6]); // Calculate the cost of disconnecting the
de -= dist(p[1],p[5]); // path segment from n[1] to n[2],
de -= dist(p[3],p[4]); // opening a space between n[3] and
de += dist(p[1],p[3]); // n[4], connecting the segment in the
de += dist(p[2],p[4]); // space, and connecting n[5] to n[6].
de += dist(p[5],p[6]);
return de;
}
void trnspt(int iorder[], int ncity, int n[])
{
int m1,m2,m3,nn,j,jj,*jorder;
jorder = ti;
m1 = 1 + ((n[2]-n[1]+ncity) % ncity); // Find number of cities from n[1] to n[2]
m2 = 1 + ((n[5]-n[4]+ncity) % ncity); // ...and the number from n[4] to n[5]
m3 = 1 + ((n[3]-n[6]+ncity) % ncity); // ...and the number from n[6] to n[3].
nn = 1;
for(j = 1; j <= m1; j++) //Copy the chosen segment.
{
jj = 1 + ((j+n[1]-2) % ncity);
jorder[nn++]=iorder[jj];
}
for(j = 1; j <= m2; j++) // Then copy the segment from n[4] to n[5].
{
jj = 1+((j+n[4]-2) % ncity);
jorder[nn++]=iorder[jj];
}
for(j = 1; j <= m3; j++) // Finally, the segment from n[6] to n[3].
{
jj = 1 + ((j+n[6]-2) % ncity);
jorder[nn++]=iorder[jj];
}
for(j = 1; j <= ncity; j++) // Copy jorder back into iorder.
iorder[j]=jorder[j];
}
int metrop(float de, float t)
{
return de < 0.0 || randf() < exp(-de/t);
}
inline int flip() { return rand() & 1; }
#define TFACTR 0.9
void aniling(int *iorder, int ncity)
{
int nover,nlimit,i1,i2;
static int n[7];
float path,de,t;
nover = 100*ncity; //Maximum number of paths tried at any temperature.
nlimit = 10*ncity; //Maximum number of successful path changes before continuing.
t = 0.5*MAX_XY;
for(int j = 1; j <= 100; j++) // Try up to 100 temperature steps.
{
int nsucc = 0, nn;
for(int k = 1; k <= nover; k++)
{
do{
n[1] = 1+(int)(ncity*randf()); // Choose beginning of segment
n[2] = 1+(int)((ncity-1)*randf()); // ..and.. end of segment.
if( n[2] >= n[1] ) ++n[2];
nn = 1+((n[1]-n[2]+ncity-1) % ncity); //nn is the number of cities not on the segment.
}while( nn < 3);
// Decide whether to do a segment reversal or transport.
if( flip() ) // Do a transport.
{
n[3] = n[2]+(int)(abs(nn-2)*randf())+1;
n[3] = 1+((n[3]-1) % ncity); // Transport to a location not on the path.
de = trncst(iorder,ncity,n); // Calculate cost.
if( metrop(de,t) ) // Consult the oracle.
{
++nsucc;
trnspt(iorder,ncity,n); // Carry out the transport.
}
}
else // Do a path reversal.
{
de = revcst(iorder,ncity,n); // Calculate cost.
if( metrop(de,t) )
{
++nsucc;
reverse(iorder,ncity,n); // Carry out the reversal.
}
}
if( nsucc >= nlimit ) break; // Finish early if we have enough successful changes.
}
// printf("\n %s %10.6f %s %12.6f \n","T =",t," Path Length =",path);
// printf("Successful Moves: %6d\n",nsucc);
t *= TFACTR; // Annealing schedule.
if( nsucc == 0 ) return; // If no success, we are done.
}
}
void doaniling(int *p)
{
int *ix = ti+n_pts+1;
for(int i = 1; i <= n_pts; i++) ix[i] = *p++;
// aproxt = droga(ix+1);
aniling(ix, n_pts);
anilingt = droga(ix+1);
}
//==============================
double bestn(int k)
{
for(int i = 0; i < n_pts; i++) ti[i] = 0;
ti[k] = 1;
double tr = 0;
TFPoint *p = pts+k;
for(int i = 1; i < n_pts; i++)
{
int nxt;
double mx = 1e6, d;
for(int j = 0; j < n_pts; j++)
if( !ti[j] && (d = dist(p, pts+j)) < mx )
{
mx = d;
nxt = j;
}
tr += mx;
p = pts+nxt;
ti[nxt] = i+1;
}
return tr+dist(pts+k, p);
}
void dobestn()
{
bestnt = 1e6;
float d;
for(int i = 0; i < n_pts; i++)
if( (d = bestn(i)) < bestnt ) bestnt = d;
}
void kola(int v)
{
double m = 2*M_PI/n_pts;
for(int i = 0; i < n_pts; i++)
{
float x = m*i, y = x;
if( v == 8 ) x *= 1.11, y *= 2;
x = cos(x); y = sin(y);
pts[i].set(MAX_XY*(0.5*x+0.5), MAX_XY*(0.5*y+0.5));
if( v != 8 )
pts[++i].set(MAX_XY*(0.36*y+0.5), MAX_XY*(0.36*x+0.5));
}
}
void siatka(int v)
{
int n = sqrt(n_pts)+0.5;
int i = 0, j = 0;
for(int k = 0; k < n_pts; k++)
{
pts[k].set(MAX_XY*i/n, MAX_XY*j/n);
if( ++i >= n )
{
i = 0; j++;
}
}
}
void mkdane(int typ=0)
{
int v = random(9);
if( v > 6 ) kola(v);
else if( v > 4 ) siatka(v);
else for(int i = 0; i < n_pts; i++)
{
float x = randf(), y = randf();
// pts[i].x = random(MAX_XY);
// pts[i].y = random(MAX_XY);
pts[i].set(x *= MAX_XY, y *= MAX_XY);
if( v & 3 == 3 ) // symetrycznie
{
pts[++i].set(MAX_XY-x, y);
pts[++i].set(MAX_XY-x, MAX_XY-y);
pts[++i].set(x, MAX_XY-y);
}
}
}
void TKomi::reset(int restartag)
{
if( !restartag )
mkdane(); // przygotowanie danych
curpop = pops;
prvpop = curpop+MAX_POP;
initpop();
if( !restartag )
{
dobestn();
doaniling(curpop->trasa);
}
}
void TKomi::skaluj()
{
double m;
if( fabs(m=minx-avex) < 1e-6 ) return;
am = m = 0.5/m; // dla najlepszego fit = 1, a dla sredniaka fit = 0.5
bm = 0.5-m*avex;
Trasa *o = curpop;
sumfit = 0;
for(int i = 0; i < pop_siz; i++, o++)
sumfit += o->fit = xfit(o->x);
}
void TKomi::initpop()
{
am = 1; bm = 0;
double s2 = 0;
sumfit = 0; avex = 0;
ZeroMemory(curpop, sizeof(Trasa)*MAX_POP);
// tworzymy losową populację
Trasa *o = curpop, *b = o;
for(int i = 0; i < pop_siz; i++, o++)
{
for(int k = 0; k < n_pts; k++) o->trasa[k] = k;
// losowa permutacja liczb
for(int k = 0; k < n_pts-1; k++)
swap(o->trasa[k], o->trasa[k+random(n_pts-k)]);
sumfit += o->decode();
avex += o->x; s2 += sqrf(o->x); // suma x i suma kwadratow do obl. sredniej i wariancji
if( *o > *b ) b = o;
}
avex /= pop_siz; minx = b->x;
// s2 = sqrt(s2/pop_siz - sqrf(avex)); // odchylenie standardowe
max_tr = minx*1.1; // maksimum y na wykresie
skaluj();
trasabest = *b;
nprob = 0;
mint[0] = minx;
popbest = kpop = 1;
}
// wypelnimy tablicę indeksow do ruletki -> przyspieszy dzialanie
void liczrul(Trasa *p, float sumfit)
{
float s = 0;
double m = P_P/sumfit;
for(int i = 0, j = 0; i < pop_siz; i++, p++)
{
int k = m*(s += p->fit);
while( j < k ) ruletka[j++] = i;
}
}
// wybor met. ruletki
Trasa *selectt(Trasa *p)
{
return p + ruletka[rand() & (P_P-1)];
}
// iteracyjna ruletka -
Trasa *select(Trasa *p, float sumfit)
{
float v = randf() * sumfit, partsum = 0;
for(int k = pop_siz; --k > 0; p++)
if( (partsum += p->fit) >= v ) break;
return p;
}
//=============== mutacje =================
typedef void (*FMut)(Node *);
FMut fmutation;
void invert(Node *t)
{
// zamiana pary losowo wybranych
swap(t[random(n_pts)], t[random(n_pts)]);
}
void revers(Node *t)
{
int i = random(n_pts), j = random(n_pts);
if( i == j && --i < 0 ) i += n_pts;
if( i > j ) swap(i,j);
TFPoint *p[6];
p[1] = pts+t[i]; p[4] = pts+t[j];
p[2] = pts+t[i+1]; p[3] = pts+t[j-1];
// (p,i,...j,n) => (p,j,...i,n)
p[0] = pts+t[(i ? i : n_pts)-1];
p[5] = pts+t[j+1 < n_pts ? j+1 : 0];
int v = 0;
if( j-i+1 != n_pts )
{
float r = dist(p[0],p[4]);
r += dist(p[1],p[5]);
r -= dist(p[0],p[1]);
r -= dist(p[4],p[5]);
if( r < 0 ) v = 1;
}
/* // i/lub rewersja dopelnienia
if( j-i-1 != 0 )
if( (dist(p[4],p[2]) + dist(p[3],p[1])) <
(dist(p[1],p[2]) + dist(p[3],p[4])) ) v |= 2;
*/
if( v )
if( v != 3 )
{
if( v == 2 ) { i++; j--; } // ...->i->|i+1-> ... ->j-1|->j->...
do swap(t[i],t[j]);
while( ++i < --j );
}
else swap(t[i],t[j]);
}
inline int rdn(int d, int n=n_pts) { return (d+n) % n; }
inline int rd(int d, int n=n_pts) { return d % n; }
void trans(Node *t)
{
int i = random(n_pts), n = 3+random(n_pts-3); // 3...N-1 -> poza i..j
int j = (i-1 + n_pts-n) % n_pts;
int k = (j+1 + random(n-1)) % n_pts; // i..j wstawimy za k
// (p,i)+(i..j)+(j,n) + (k,l) => (p,n)+(k,i)+(i..j)+(j,l)
// r = (p,n)+(k,i)+(j,l) -(p,i)-(j,n)-(k,l)
TFPoint *p[6];
p[1] = pts+t[i]; p[2] = pts+t[j]; // i, j
p[0] = pts+t[rdn(i-1)]; p[3] = pts+t[rd(j+1)]; // p, n
p[4] = pts+t[k]; p[5] = pts+t[k=rd(k+1)]; // k, l; teraz k = k+1 = l
float r = dist(p[0], p[3]); // p-n
r += dist(p[4], p[1]); // k-i
r += dist(p[2], p[5]); // j-l
r -= dist(p[0], p[1]); // p-i
r -= dist(p[2], p[3]); // j-n
r -= dist(p[4], p[5]); // k-l
if( r < 0 ) // gdy skraca
{
for(int i = 0; i < n_pts; i++) ti[i] = t[i];
// n..k = bz.
// l... = i..j
// ...j = l..p
int l = k;
while( 1 ) // i..j
{
t[l] = ti[i];
if( ++l == n_pts ) l = 0;
if( i == j ) break;
if( ++i == n_pts ) i = 0;
}
while( 1 ) // l..p
{
t[l] = ti[k];
if( l == j ) break;
if( ++l == n_pts ) l = 0;
if( ++k == n_pts ) k = 0;
}
}
}
void mutmix(Node *t)
{
revers(t);
trans(t);
}
// krzyżowanie PMX - wymiana losowego podciągu
void crosspmx(Node *c1, Node *c2, Node *p1, Node *p2)
{
int *i1 = ti, *i2 = ti+n_pts;
for(int i = 0; i < n_pts; i++)
{
i1[c1[i] = p1[i]] = i;
i2[c2[i] = p2[i]] = i; // i[k] = poz k -> c[poz] = k
}
if( flip(pcross) )
{
int i = random(n_pts);
int n = random(n_pts);
if( i > n ) swap(i,n);
n -= i;
do
{
int k = i1[c2[i]];
Node c = c1[i];
c1[i] = c1[k]; c1[k] = c;
k = i2[c];
c = c2[i];
c2[i] = c2[k]; c2[k] = c;
i++;
}while( --n > 0 );
}
if( flip(pmutation) ) fmutation(c1);
if( flip(pmutation) ) fmutation(c2);
}
// krzyżowanie CX - cykliczne
void crosscx(Node *c1, Node *c2, Node *p1, Node *p2)
{
int *&ix = ti;
for(int i = 0; i < n_pts; i++)
{
c1[i] = p2[i];
ix[c2[i] = p1[i]] = i; // ix[k] = poz k -> c2[poz] = k
}
if( flip(pcross) )
{
int k = 0;
do // idziemy az do konca cyklu
{
int c = c1[k];
c1[k] = c2[k];
c2[k] = c;
k = ix[c];
}while( k != 0 );
}
if( flip(pmutation) ) fmutation(c1);
if( flip(pmutation) ) fmutation(c2);
}
void TKomi::generacja()
{
Trasa *c = prvpop, *pb = c;
prvpop = curpop; curpop = c;
liczrul(prvpop, sumfit);
double sf = 0, sx = 0;
int n = pop_siz >> 1;
do{
Trasa *a = selectt(prvpop),
*b = selectt(prvpop);
crosscx(c->trasa,(c+1)->trasa, a->trasa,b->trasa);
sf += c->decode(); sx += c->x;
if( *c > *pb ) pb = c;
c++;
sf += c->decode(); sx += c->x;
if( *c > *pb ) pb = c;
c++;
}while( --n > 0 );
sumfit = sf;
kpop++;
if( *pb > trasabest )
{
trasabest = *pb;
popbest = kpop;
}
minx = pb->x;
avex = sx/pop_siz;
int k = kpop >> Sh_P; // srednia z 8
if( k < N_P )
{
mint[nprob] += minx;
if( (kpop & (M_Pr-1)) == 0 ) mint[++nprob] = 0;
}
if( !(kpop & 31) ) skaluj();
}
void TKomi::setMut(int im)
{
if( im == 0 ) fmutation = invert;
else if( im == 1 ) fmutation = revers;
else fmutation = trans;
}
int TKomi::oblicz(int iter)
{
while( --iter >= 0 )
{
generacja();
}
return nprob > 250 && kpop-popbest > 500;
}
//---------------------------------------------------------------------------
using namespace std;
void czytaj(char* fname)
{
fstream f(fname, ios_base::in);
f.unsetf(ios_base::skipws);
int n = 0, v = 0;
float *p = (float*)pts;
f.peek();
*p = 0;
for(int l = 0; l < 2*MAX_PTS && !f.eof(); )
{
char c;
f >> c;
if( isdigit(c) )
{
v = 10*v + c - '0';
n++;
}
else if( n )
{
*p++ = v;
v = n = 0;
l++;
}
}
if( n ) p++;
n_pts = (p - (float*)pts)/2;
if( n_pts < MIN_PTS ) mkdane();
}
void zapisz(char* fname)
{
fstream f(fname, ios_base::out|ios_base::trunc);
// f.precision(0);
int n = 0;
for(int l = 0; l < n_pts ; l++)
{
f << (int)pts[l].x << ',' << (int)pts[l].y << '\t';
if( ++n == 8 ) { n = 0; f << endl; }
}
}
#pragma package(smart_init)
Oryginalnie był to algorytm genetyczny, no ale dość słabo się spisywał... chociaż coś rzeczywiście tam robił;
no ale wyparzanie lepiej.
Chyba w oparciu o ten opis to robiłem, więc możesz sobie poczytać:
http://bookstaber.com/david/opinion/SATSP.pdf
Było to robione w Buildrze C++, ale teraz nie mam go zainstalowanego...
program pięknie rysował te trasy na żywo w czasie obliczania.