2009-07-01

Optimizacija mravlje kolonije

Sad prelazimo sa analogije metalurških procesa na biologiju :) E ta veštačka inteligencija, uvek traže inspiraciju u nečemu postojećem. Ideja je nastala na osnovu proučavanja konunikacije i upravljanja koji se dešava među insektima. Mravi kad traže neki izvor ispuštaju feromone. Feromoni vremenom isparavaju a količina feromona na nekoj stazi zavisi i od dužine staze. Mrav koji nađe kraći put na tom putu će ostaviti veću količinu feromona koji će drugim mravima poslužiti za navođenje na bolju putanju.

Pomoću ove ideje moguće je napisati heuristiku za rešavanje problema trgovačkog putnika (nalaženje što kraće Hamiltonove konture). 

Dovoljno je imati 3 proste formule, toliko proste da su matematičari morali da ih zakomplikuju lupanjem nekih čudnih oznaka, samo im rusko slovo Я fali. To me podseti na epizodu kad sam polagao statistiku kod jedne nemaštovite asistentkinje kojoj je bilo čudno što formulu uslovne verovatnoće nisam izveo identično kao u knjizi, tj. što nisam bubao već pravio delikt mišljenja. Pa sam joj rekao da dok god je formula ispravna mogu da koristim i Э kao oznaku :D Matematičar su uvek najveći problem matematike.

Zato ću ja formule prikazati uprošteno (a i ne znam kako bih napisao sve one indekse, stepene, klince, palce). Prva formula služi da se izračuna količina feromona za neku ivicu (direktnu vezu dva grada). Formula je rekurentna (programerski rečeno rekurzivna) jer se nova količina feromona izračunava preko stare prostim dodavanjem:

kolicina_feromona = kolicina_feromona + Q/duzina_Hamiltonove_konture *RO

Ili uprošćeno:

kolicina_feromona +=  Q/duzina_Hamiltonove_konture *RO

Q je neka konstanta koja označava količinu feromona kojom svaki mrav raspolaže. Pa mrav koji nađe kraću putanju ima ivicama da doda više feromona jer je Q/mala_vrednost veći broj nego Q/velika_vrednost. Konstanta RO je između 0 i 1 i služi da odredi koliko novog feromona treba dodati. Pošto je formula rekurentna mora postojati početna vrednost pa se uzima da je kolicina_feromona = 1 / broj_gradova

Druga formula simulira isparavanje feromona:

kolicina_feromona = kolicina_feromona * (1 - RO)

ili:

kolicina_feromona *= (1 - RO)

RO je ista konstanta kao iz prethodne formule. Npr. ako je RO = 0,7 to znači da će se na putanju kojom je mrav prošao staviti 70% vrednosti Q/duzina_Hamiltonove_konture, a od stare količine feromona 70% će ispariti.

I treća formula izračunava verovatnoću da se krene nekom ivicom tj. da se obiđe neki neposećeni grad iz tekućeg grada. Za svaku moguću ivicu (tj. za svaku direktnu vezu tekućeg i neposećenih gradova) se izračunava poželjnost. Poželjnost je proporcionalna količini feromona na ivici i njenoj kratkoći. Kratkoću definišemo kao 1 / duzina_ivice a količinu feromona čitamo direktno. Udeo pozeljnosti i-te moguće ivice u zbiru pozeljnosti svih mogućih ivica je upravo verovatnoća da se mrav zaputi tom ivicom. Količinu feromona matematičari označavaju ubercool slovom Tau, a kratkoću slovom Eta. Formula je (** je operator stepenovanja):

P = (Tau**ALFA * Eta**BETA) / SUM(Tau**ALFA * Eta**BETA)

Da bi se kontrolisao značaj kratkoće i feromona uvode se ALFA i BETA parametri kojim se stepenuju količine feromona (Tau) i kratkoće (Eta). Pošto su i Tau i Eta manji od 0 znači da što je stepen veći to se umanjuje značaj te karakteristike u poželjnosti. 


I konačno kod:



import random, copy, math, time

class Grad(object):
def __init__(self, x, y, ime):
self.x = x
self.y = y
self.ime = ime

def rastojanje(self, grad):
"""
rastojanje ovog grada od grada primljenog kao prarametar
preko euklidske metrike (odnosno funkcije hypot).
"""
return math.hypot(self.x - grad.x, self.y - grad.y)

def __str__(self):
return "%s(%d, %d)" % (self.ime, self.x, self.y)

class Ivica(object):
def __init__(self, feromoni, duzina):
self.feromoni = feromoni
self.duzina = duzina

class Mrav(object):
# atribut klase u kojem se cuva broj svih stvorenih mrava
brojac = -1
gradovi = None
ivice = None

def inicijalizuj_atribute(self):
self.tek_grad = self.pocetni_grad
self.putanja = [self.tek_grad]
self.duzina_putanje = 0.0

def __init__(self, gradovi, alfa, beta, ro, q):
Mrav.ALFA = alfa
Mrav.BETA = beta
Mrav.RO = ro
Mrav.Q = q

if not Mrav.gradovi:
Mrav.gradovi = gradovi

if not Mrav.ivice:
Mrav.ivice = {}
init_feromoni = 1.0 / len(Mrav.gradovi)
for i in Mrav.gradovi:
for j in Mrav.gradovi:
if i != j :
rastojanje = i.rastojanje(j)
Mrav.ivice[i, j] = Ivica(init_feromoni, rastojanje)

Mrav.brojac += 1
self.pocetni_grad = Mrav.gradovi[Mrav.brojac % len(Mrav.gradovi)]

self.inicijalizuj_atribute()

def _tau_x_eta(self, do_grad):
# izracunava vrednost pozeljnosti ivice od tekuceg grada do grada
# primljenog kao parametar. Pozeljnost ivice je direktno proporcionalna
# kolicini feromona na njoj (sto se oznacava sa tau), a obratni proporcionalna
# njenoj duzini (eta). Stepenovanjem tau i eta sa parametrima ALFA i BETA
# omogucava nam da dajemo prednost jednom ili drugom kriterijumu pozeljnosti
# i to tako sto BETA oznacava koliko umanjujemo znacaj blizine a ALFA
# koliko umanjujemo znacaj feromona. Npr. za ALFA: 1 i BETA: 5 znaci da
# nam je vaznija kolicina feromona od blizine
ivica = Mrav.ivice[self.tek_grad, do_grad]
tau = ivica.feromoni
eta = 1.0 / ivica.duzina

return tau**Mrav.ALFA * eta**Mrav.BETA

def _odaberi_sledeci_grad(self):
imenilac = 0.0

# svi gradovi razlika poseceni = neposeceni, skupovna operacija
neposeceni_gradovi = set(Mrav.gradovi).difference(self.putanja)

# ako nema neposecenih vrati None
if not neposeceni_gradovi:
return None

# za svaki neposeceni grad izracunaj _tau_x_eta i saberi
# tj. racunamo zbir svih pozeljnosti posete neposecenim gradovima
for grad in neposeceni_gradovi:
imenilac += self._tau_x_eta(grad)

# za svaki neposecen grad izracunaj verovatnocu posete p
# i pokusaj da je ostvaris, ako je ne ostvaris probaj
# ponovo sve iz pocetka dok konacno ne uspes u tome.
# Verovatnoca posete se izravunava tako sto podelimo pozeljnost posete
# neposecenog grada sa zbirom svih pozeljnosti posete neposecenim gradovima
while True:
for grad in neposeceni_gradovi:
p = self._tau_x_eta(grad) / imenilac

if random.random() < p:
return grad

def pusti_mrava(self):
sled_grad = self._odaberi_sledeci_grad()

# dok postoji sledeci grad za obilazak stavi ga u putanju obishenih
# azuriraj duzinu putanje mrava i tekuci grad a zatim nadji sledeci grad
while sled_grad:
self.putanja.append(sled_grad)
self.duzina_putanje += Mrav.ivice[self.tek_grad, sled_grad].duzina
self.tek_grad = sled_grad
sled_grad = self._odaberi_sledeci_grad()

# sad je sled_grad None znaci da smo prosli sve gradove
# a Hamiltonovu konturu zatvaramo vezujuci 1. i poslednji
# poseceni grad koji se cuvaju u self.putanja (0. i -1. poziciju)
self.duzina_putanje += Mrav.ivice[self.putanja[0], self.putanja[-1]].duzina

def pospi_feromone(self):
br_gradova = len(self.putanja)

# svakoj ivici putanje mrava (Hamiltonove konture) u oba smera dodaj
# odgovarajucu kolicinu feromona po formuli Q / duzina_putanje * RO
for i in xrange(br_gradova):
od_grada = self.putanja[i]
do_grada = self.putanja[(i + 1) % br_gradova]

Mrav.ivice[od_grada, do_grada].feromoni += Mrav.Q / self.duzina_putanje * Mrav.RO
Mrav.ivice[do_grada, od_grada].feromoni = Mrav.ivice[od_grada, do_grada].feromoni

@staticmethod
def ispari_feromone():
# smanji kolicinu feromona na svakoj ivici po rekurentnoj formuli feromoni *= 1 / RO
for ivica in Mrav.ivice.itervalues():
ivica.feromoni *= (1.0 - Mrav.RO)

def koordinate(self):
x, y = [], []

for grad in self.putanja:
x.append(grad.x)
y.append(grad.y)

return x, y

# toString u javi
def __str__(self):
return "duzina puta: %f" % self.duzina_putanje

# compare u javi
def __cmp__(self, drugi):
return cmp(self.duzina_putanje, drugi.duzina_putanje)


class MravljiSimulator(object):
def __init__(self, gradovi, ALFA = 1.0, BETA = 5.0, RO = 0.5, Q = 100.0):
self.mravi = []
for i in xrange(len(gradovi)):
self.mravi.append(Mrav(gradovi, ALFA, BETA, RO, Q))
self.najbolji_mrav = self.mravi[0] # proizvoljnog mrava proglasimo najboljim

def start(self, ponavljanja=20):
# za svako ponavljanje
for i in xrange(ponavljanja):
# pusti svakog mrava da napravi Hamiltonovu konturu
for mrav in self.mravi:
mrav.pusti_mrava()

Mrav.ispari_feromone()

# svaki mrav nek pospe feromone na svojoj konturi
for mrav in self.mravi:
mrav.pospi_feromone()

najuspesniji_mrav_generacije = min(self.mravi)

# pamtimo sveukupno najuspesnijeg mrava
self.najbolji_mrav = copy.deepcopy(
min(self.najbolji_mrav, najuspesniji_mrav_generacije))

# pripremamo mrave za ponavljanje tako sto ih resetujemo
for mrav in self.mravi:
mrav.inicijalizuj_atribute()

def stampaj_resenje(self):
if self.najbolji_mrav:
try:
# ako je instalirana matplot biblioteka koristi je za graficki prikaz
import matplotlib.pyplot as plt

x, y = self.najbolji_mrav.koordinate()

# 1. tacku dodajemo ponovo da bi zatvorili konturu od poslednje ka 1. tacki
x.append(x[0])
y.append(y[0])

plt.title('prikaz nadjenog resenja')
plt.plot(x, y, 'g', marker='.', markerfacecolor='r', markersize=20, label='najbolje r.')
plt.xlabel('x koordinata')
plt.ylabel('y koordinata')
plt.text(2, 95, 'najbolja duzina: %f.2' % self.najbolji_mrav.duzina_putanje)
plt.legend()

for grad in self.najbolji_mrav.putanja:
plt.text(grad.x, grad.y - 2.5, grad.ime)

plt.axis([min(x) -1, max(x) + 1, min(y) -1, max(y) +1])

plt.show()
except ImportError:
# nema matplot biblioteke stampaj resenje rucno
print "Resenje"
print self.najbolji_mrav.putanja
print "duzina puta je: ", self.najbolji_mrav.duzina_putanje


def main():
try:
import psyco # ako je instalirana biblioteka koristi je (ubrzava izvsenje koda)
psyco.full()
except ImportError:
pass

gradovi = [Grad(1, 2, 'a'), Grad(29, 21, 'b'),Grad(100, 60, 'c'), Grad(12, 86, 'd'),
Grad(92, 46, 'e'), Grad(83, 38, 'f'), Grad(55, 36, 'g'), Grad(71, 99, 'h'),
Grad(12, 41, 'i'), Grad(34, 48, 'j'), Grad(69, 33, 'k'), Grad(78, 10, 'l'),
Grad(86, 68, 'm'), Grad(79, 27, 'n'), Grad(22, 69, 'o'), Grad(75, 55, 'p'),
Grad(51, 68, 'r'), Grad(91, 23, 's'), Grad(22, 42, 't'), Grad(47, 80, 'u'),
Grad(60, 10, 'w'), Grad(91, 79, 'v'), Grad(5, 66, 'x'), Grad(42, 90, 'y'),
Grad(23, 59, '1'), Grad(46, 83, '2'), Grad(93, 63, '3'), Grad(47, 17, '4'),
Grad(53, 79, '5'), Grad(76, 23, '6'), Grad(91, 62, '7'), Grad(44, 97, '8'),
]

simulator = MravljiSimulator(gradovi, ALFA=1.0, BETA = 2.0, RO=0.5)

t0 = time.time()
simulator.start(ponavljanja=50)
print "Izvrsenje trajalo", time.time() - t0, "sekundi"

simulator.stampaj_resenje()

if __name__ == "__main__":
main()





Gradovi koje sam ovde koristio su isti oni koje sam korisitio sa Simuliranim kaljenjem. Rastojanja gradova se unapred izračunaju da bi se program brže izvršavao. Brzo se dobijaju dobra rešenja. Otprilike da je putanja koju pronađe u proseku dugačka oko 530. A svaki 15-20 put ubode se i rešenje čiju sam sliku postavio. Za razliku od Simuliranog kaljena ovde ima mnoštvo parametara sa kojima se treba igrati da se nađe njihova dobra vrednost.

2009-06-23

Kako se kale kraljice

Još jedan problem s kojim se ljudi koji se bave umetnošću programiraju vaćaju u koštac je problem N kraljica. Na tabli N x N treba postaviti N kraljca a da jedna drugu ne jedu.

I ovaj problem se može rešiti Simuliranim kaljenjem. Funkcija cilja je u ovom slučaju broj sukoba, samo za razliku od prošlog primera sa trgovcem ovde se zna da je optimalno rešenje kad je f. cilja = 0.

Počentno rešenje se lako nalazi. Prvo ne sme se desti da više od jedne kraljce bude u bilo kom redu i bilo kojoj koloni. Dakle najprostije je da se poređaju po dijagonali. U mom rešenju se takvo rešenje tvikuje par puta i dobija se početno rešenje.

Tvikovanje se sastoji u nasumičnoj zameni 2 kolone.

Implementacija je takva da se u klasi Resenje u nizu koji se zove kraljice čuva pozicija kraljice u redu za svaku kolonu. Npr. niz [1, 0, 2] znači da se u 1. koloni kraljica nalazi u 2. redu, 2. koloni u redu 1, a u 3. koloni u redu 3 (i u Pythonu su prvi elementi niza na poziciji 0). Dakle u početnom rešenju ćemo pre tvikovanja imati ovakav niz: [0, 1, 2, 3, ...] jer kraljice postavljamo po dijagonali.

Dalje formiramo tabelu N x N i koja se puni None objektima (u pythonu je sve objekat pa i Null vrednost ) osim što se na dijagonalama postavi string K kao oznaka da je tu kraljica.

Samo simulirano kaljenje je slično kao u prošlom primeru. Jedino se manje vrši resetovanje tekućeg rešenja jer je takav zadatak da se tvikovanje može obaviti na već tvikovano rešenje bez većeg problema.

A pošto znamo koje je optimalno rešenje, kaljenje se ponavlja iznova sve dok ne dođemo do tog rešenja kad se štampa optimalni raspored na tabli.

import random, copy, math, time

class Resenje(object):
def __init__(self, br_kraljica):
self.br_kraljica = br_kraljica
# resenje: pozicija kraljice u koloni 0 do len - 1
self.kraljice = range(br_kraljica)
# energija: treba da je 0
self.__broj_sukoba = None

## postavljamo pocetno resenje
# generisemo tabelu br_kraljica X br_kraljica popunjenu None elementima
self.tabla = [[None for i in range(br_kraljica)] for j in range(br_kraljica)]

# postavimo kraljice po dijagonali
for i in xrange(br_kraljica):
self.tabla[i][i] = 'K'

# promesamo kolone
for i in xrange(br_kraljica):
self.tvikuj_resenje()

# f-ja cilja, optimalno resenje ima 0 sukoba
def get_broj_sukoba(self):
# Ovo je optimizacija da se broj sukoba racuna samo ako vec nije izrazcunat
if not self.__broj_sukoba:
self.__broj_sukoba = self.__racunaj_broj_sukoba()
return self.__broj_sukoba

# racuna energiju resenja, prakticno funkciju cilja koja se minimizuje
# a je ovde broj sukoba svih kraljica na tabli
def __racunaj_broj_sukoba(self):
# povratna vr.
# zato sto se sukobljava sam sa sobom i to pri proveri svake od 4
# dijagonala a to se ne racuna
br_sukoba = -4 * self.br_kraljica

dx = (-1, 1, -1, 1) # pomeraj koordinate x za kretanje po dijagonalama
dy = (-1, 1, 1, -1) # pomeraj koordinate y -//-

# provera dijagonala
for x, y in enumerate(self.kraljice): # (x,y) koord. kraljica
for i in xrange(4): # za sva 4 kraka dijagonala
pomx, pomy = x, y # pretraga krece od koordinata kraljice
# dok se pomeranjem po dijagonali nalazimo na tabli ispitujemo
# da li smo nabasali na kraljicu. Svaki od 4 ciklusa za 4 dijagonale
# nailazi na sopstvenu kraljicu pa je zato br_sukoba inicijalizovan
# na -4 * self.br_kraljica
while (0 <= pomx < self.br_kraljica) and (0 <= pomy < self.br_kraljica):
if self.tabla[pomx][pomy] == 'K':
br_sukoba += 1
pomx += dx[i]
pomy += dy[i]

return br_sukoba

broj_sukoba = property(get_broj_sukoba)

# zamenjuje 2 kolone
def tvikuj_resenje(self):
r1 = random.randrange(self.br_kraljica)
r2 = random.randrange(self.br_kraljica)

# kolone moraju biti razlicite
while r2 == r1:
r2 = random.randrange(self.br_kraljica)

# zameni redove r1 i r2
self.tabla[r1], self.tabla[r2] = self.tabla[r2], self.tabla[r1]

# zameni elemente
self.kraljice[r1], self.kraljice[r2] = self.kraljice[r2], self.kraljice[r1]

# prethodno izracunata duzina puta vise ne vazi pa se setuje na None
self.__broj_sukoba = None

def __str__(self):
rez = ''
for red in self.tabla:
s = ''
for element in red:
if element: s += element
else: s += '.'
rez += s
rez += '\n'
return rez


class SimulatorKaljenja(object):
def __init__(self, br_kraljica):
self.br_kraljica = br_kraljica
self.resenje = None

def start(self, pocetna_temperatura=90.0, alpha=0.997, ponavljanja=None):
najbolje_resenje = Resenje(self.br_kraljica)
tekuca_temp = pocetna_temperatura
if not ponavljanja: ponavljanja = self.br_kraljica / 2


while tekuca_temp > 0.10 and najbolje_resenje.broj_sukoba > 0:
privremeno_resenje = copy.deepcopy(najbolje_resenje)

# Monte Carlo na tekucoj temperaturi
for i in xrange(ponavljanja):
privremeno_resenje.tvikuj_resenje()

delta_e = privremeno_resenje.broj_sukoba - najbolje_resenje.broj_sukoba

# ako je resenje bolje uzimamo ga a ako nije tada uz
# odredjenu verovatnocu odabiramo losije resenje kao tekuce
if delta_e < 0.0 or (math.exp(-delta_e / tekuca_temp) > random.random()):
if privremeno_resenje.broj_sukoba < najbolje_resenje.broj_sukoba:
najbolje_resenje = copy.deepcopy(privremeno_resenje)
else:
privremeno_resenje = copy.deepcopy(najbolje_resenje)

tekuca_temp *= alpha

self.resenje = najbolje_resenje

def main():
try:
import psyco
psyco.full()
except ImportError:
pass

simulator = SimulatorKaljenja(32)

t0 = time.time()

i = 1
while True:
print 'ciklus', i
i += 1
simulator.start()
print 'najbolje resenje', simulator.resenje.broj_sukoba
if simulator.resenje.broj_sukoba == 0:
break

print simulator.resenje
print "Izvrsenje trajalo", time.time() - t0, "sekundi"


if __name__ == "__main__":
main()

Primer rezultat ovoga koda je:

ciklus 1 
najbolje resenje 2 
ciklus 2 
najbolje resenje 0 
.K.............................. 
....................K........... 
..................K............. 
............K................... 
..K............................. 
...........................K.... 
.................K.............. 
K............................... 
..........K..................... 
.....K.......................... 
...............................K 
...................K............ 
..........................K..... 
.............K.................. 
.....................K.......... 
............................K... 
...............K................ 
.............................K.. 
......................K......... 
......K......................... 
.........................K...... 
.......K........................ 
..............K................. 
...........K.................... 
..............................K. 
....K........................... 
................K............... 
........................K....... 
.........K...................... 
.......................K........ 
...K............................ 
........K....................... 

Izvrsenje trajalo 35.1647040844 sekundi

2009-06-22

Kako se kalio čelik

Nešto mi je ovijeh dana dosadno bilo, a dokon um đavolje dvorište. Te se dovatim jedne knjige o AI algoritmima. U pa tu ima svakakve 'euristike, ovi iz Operacionih bi zinuli. Ovde ću kako budem stizao da postavljam neke algoritme. Knjiga ima neke C kodove, sad malo ko parla C, a i ti kodovi su tako odvratni i u knjizi dati iz delova i nepotpuno. Pa sam ih ja sredio u mom omiljenom prog. jeziku.

Prvi algoritam je Simulirano kaljenje. Ima nešto o tome u knjizi iz Operacionih, ali to su pisali matematičari (čitaj: dosadni tipovi). Ideja je da se podražava proces kaljenja. To smo valjda učili iz Teknologije, eto neke koristi i od toga Pa kako se kali? Podigne se temperatura pa se čuka metalni predmet i izvaja u neki oblik, pa se polako 'ladi a i dalje čuka.
U ovom algoritmu se prvo nađe neko najobičnije rešenje, zatim se to rešenje tvikuje (malo se promeni) i gleda kakvo je u odnosu na prethodno. Ako je bolje super, usvajamo ga, ako nije onda po nekoj verovatnoći usvaja. Jer lošije rešenje može u nekoj sledećoj iteraciji voditi u bolje. Usvajanjem ponekog lošijeg rešenja se izbegava nalaženje nekog lokalnog minimuma koji je daleko od globalnog.

E sad da se ne bi mlogo lutalo ta verovatnoća da se usvoji lošije rešenje se iz ciklusa u ciklus smanjuje. To podražava smanjenje se temperatura metala pri kaljenju, što je metal ladniji teže ga je oblikovati. A verovatnoća usvajanja lošijeg rešenja zavisi i od toga koliko je to rešenje lošije. Ako je mlogo loše onda je mala verovatnoća da ga usvojimo.
Znam da vam ništa nije jasno Evo primer primene algoritma. Setite se problema trgovačkog putnika. Postoje algoritmi ali su oni NP-kompletni (u prevodu, teški u p.m za više od 30 čvorova). E sad Simuliranim kaljenjem se dobija rešenje blisko optimalnom ako ne i optimalno (ako vas baš ukenja).

I kako to radi? Pa uzme se neka putanja. Može nasumična može sistemom najbližih čvorova (ovo je bolje početno rešenje). Pa se onda u tom rešenju zamene mesta nekih čvorova i izračuna nova dužina i poredi sa prethodnom. Ako je bolja usvaja se rešenje, ako nije ide po onoj verovatnoći. I tako se iz ciklusa u ciklus smanjuje temperatura sve dok ne padne na neku zaustavnu a pamti se najbolje rešenje. Prosto ko pasulj!

Nemojte da vas uplaši kod, to je prost jezik, a najveći deo koda odlazi na štampanje, što vam nije potrebno da razumete jer je to pythonizovan matlab. Prvo pravimo klasu Grad koja će predstavljati jedan čvor u grafu (__init__ metoda je konstruktor u Pythonu) i definišemo funkciju koja za dati grad iz liste gradova nalazi najbliži:

class Grad(object):
def __init__(self, x, y, label):
self.x = x
self.y = y
self.label = label

def rastojanje(self, grad):
"""
rastojanje ovog grada od grada primljenog kao prarametar
preko euklidske metrike (odnosno funkcije hypot).
"""
return math.hypot(self.x - grad.x, self.y - grad.y)

def __str__(self):
return "%s(%d, %d)" % (self.label, self.x, self.y)


def _najblizi(pivot, gradovi):
rez = None
min_rastojanje = sys.maxint
for grad in gradovi:
rastojanje = grad.rastojanje(pivot)
if rastojanje < min_rastojanje:
min_rastojanje = rastojanje
rez = grad
return rez

Potom prvimo klasu za rešenje koje se sastoji od redosleda gradova koje treba obići. Ima metodu za dužinu rešenja (vrednost funkcije cija što bi rekli ćelavi matematičari), metodu za tvikovanja i jednu pomoćnu za grafike. Objekat u konstruktoru automatski nalazi početno rešenje obilaska po metodu idem od čvora do prvog najbližeg dok ne prođem kroz sve:

class Resenje(object):
def __init__(self, gradovi=[]):
self.gradovi = [] # gradovi
self.__duzina_puta = None # energija resenja

tek = gradovi.pop(0) # ukloni se prvi od primenljenih
self.gradovi.append(tek) # i stavi kao pocetni

while gradovi:
tek = _najblizi(tek, gradovi) # nadje se najblizi tekucem
gradovi.remove(tek) # on se ukloni iz primeljenih
self.gradovi.append(tek) # i stavi kao sledeci u putanji

def get_duzina_puta(self):
# ako je duzina puta izracunata (prethodno je pozivana grad.duzina_puta)
# i ako u medjuvremenu nije bilo zamene tada je ta duzina i dalje postoji
# i ispravna je. Ovo je jedna optimizacija.
if not self.__duzina_puta:
self.__duzina_puta = self.__racunaj_put()
return self.__duzina_puta

def __racunaj_put(self):
duzina = 0.0

for i in range(len(self.gradovi)):
tek_grad = self.gradovi[i]
sledeci_grad = self.gradovi[(i+1) % len(self.gradovi)]
duzina += tek_grad.rastojanje(sledeci_grad)

return duzina

# duzina_puta je property kao u C#, get metoda je get_duzina_puta
# set metoda ne postoji (read-only property)
duzina_puta = property(get_duzina_puta)

def zameni_put(self):
g1 = random.randrange(len(self.gradovi))
g2 = random.randrange(len(self.gradovi))

while g2 == g1: # obezbedjivanje da g2 != g1
g2 = random.randrange(len(self.gradovi))

# zamena u jednoj liniji bez pomocne promenjive
self.gradovi[g1], self.gradovi[g2] = self.gradovi[g2], self.gradovi[g1]

# prethodno izracunata duzina puta vise ne vazi pa se setuje na None
self.__duzina_puta = None

# pomocna metoda za stampanje resenja
def koordinate(self):
x, y = [], []

for grad in self.gradovi:
x.append(grad.x)
y.append(grad.y)

return x, y
I konačno klasa u kojoj se simulirano kaljenje obavlja. U konstruktoru se prihvata početno rešenje a u metodi start() se dešava ćudo, sa namernim ć :) :

class SimulatorKaljenja(object):
def __init__(self, pocetno_resenje):
self.pocetno_resenje = copy.deepcopy(pocetno_resenje)
self.resenje = None
self.najbolje_resenje = None
self.duzine = []
self.temperature = []

def start(self, pocetna_temperatura=100, alpha=0.999, ponavljanja=10):
tekuce_resenje = copy.deepcopy(self.pocetno_resenje)
privremeno_resenje = copy.deepcopy(self.pocetno_resenje)
najbolje_resenje = copy.deepcopy(self.pocetno_resenje)
tek_temperatura = pocetna_temperatura
ZAUSTAVNA_TEMPERATURA = 0.4

while tek_temperatura > ZAUSTAVNA_TEMPERATURA:

# pamtimo kretanje duzine i temperature kroz iteracije
# da bismo ih kasnije prikazali na grafikonu
self.duzine.append(tekuce_resenje.duzina_puta)
self.temperature.append(tek_temperatura)

for i in xrange(ponavljanja):
privremeno_resenje.zameni_put()

delta_e = privremeno_resenje.duzina_puta - tekuce_resenje.duzina_puta

# ako je privremeno_resenje bolje (delta_e < 0) uzimamo ga
# a ako nije tada uz odredjenu verovatnocu po formuli:
# e^(-delta_e / tek_temperatura) odabiramo losije resenje kao tekuce
# ta verovatnoca je obrnuto srazmerna tek_temperaturi i sto je to novo
# resenje losije (sto je delta_e vece) ono se teze usvaja
if delta_e > 0: print delta_e, tek_temperatura, math.exp(-delta_e / tek_temperatura)
if delta_e < 0.0 or ( math.exp(-delta_e / tek_temperatura) > random.random() ):
tekuce_resenje = copy.deepcopy(privremeno_resenje)
if tekuce_resenje.duzina_puta < najbolje_resenje.duzina_puta:
najbolje_resenje = copy.deepcopy(tekuce_resenje)
else:
# resetujemo privremeno resenje i probamo dalje
privremeno_resenje = copy.deepcopy(tekuce_resenje)

tek_temperatura *= alpha

self.resenje = tekuce_resenje
self.najbolje_resenje = najbolje_resenje

return tekuce_resenje

Ovoj se klasi dodaje i metoda za prikaz rešenja u grafičkom obliku preko matplot biblioteke (i sve to besplatno i slobodno):

def stampaj_resenje(self):
if self.resenje:
try:
# ako je instalirana matplot biblioteka koristi je za graficki prikaz
import matplotlib.pyplot as plt

x1, y1 = self.pocetno_resenje.koordinate()
x2, y2 = self.resenje.koordinate()
x3, y3 = self.najbolje_resenje.koordinate()

# 1. tacku dodajemo ponovo da bi zatvorili konturu od poslednje ka 1. tacki
x1.append(x1[0])
y1.append(y1[0])
x3.append(x3[0])
y3.append(y3[0])

# prvi grafikon
plt.figure(1)
plt.title('prikaz pocetnog i nadjenog resenja')
plt.plot(x1, y1, 'k--', label='pocetno r.')
plt.plot(x3, y3, 'g', marker='.', markerfacecolor='r', markersize=20, label='najbolje r.')
plt.text(2, 95, 'pocetna duzina: %f.2' % self.pocetno_resenje.duzina_puta)
plt.xlabel('x koordinata')
plt.ylabel('y koordinata')
plt.text(2, 87, 'najbolja duzina: %f.2' % self.najbolje_resenje.duzina_puta)
plt.legend()

for grad in self.resenje.gradovi:
plt.text(grad.x, grad.y - 2.5, grad.label)

plt.axis([min(x1) -1, max(x1) + 1, min(y1) -1, max(y1) +1])

# drugi grafikon
plt.figure(2)

# 1. subplot
plt.subplot(2, 1, 1)
plt.title('kretanje temperature po iteracijama')
plt.plot(self.temperature, 'r')
plt.xlabel('iteracija')
plt.ylabel('temperatura')

# 2. subplot
plt.subplot(2, 1, 2)
plt.title('vrednost f. cilja po iteracijama (duzina puta)')
plt.plot(self.duzine)
plt.xlabel('iteracija')
plt.ylabel('duzina puta')

plt.show()
except ImportError:
# nema matplot biblioteke stampaj resenje
print "Pocetno resenje"
print self.pocetno_resenje.gradovi
print "duzina puta je: ", self.pocetno_resenje.duzina_puta
print "\nResenje"
print self.najbolje_resenje.gradovi
print "duzina puta je: ", self.najbolje_resenje.duzina_puta

I na kraju puštamo u pogon ovu skalameriju koja uz to i radi :)

def main():
try:
import psyco # ako je instalirana biblioteka koristi je (ubrzava izvsenje koda)
psyco.full()
except ImportError:
pass

gradovi = [Grad(1, 2, 'a'), Grad(29, 21, 'b'),Grad(100, 60, 'c'), Grad(12, 86, 'd'),
Grad(92, 46, 'e'), Grad(83, 38, 'f'), Grad(55, 36, 'g'), Grad(71, 99, 'h'),
Grad(12, 41, 'i'), Grad(34, 48, 'j'), Grad(69, 33, 'k'), Grad(78, 10, 'l'),
Grad(86, 68, 'm'), Grad(79, 27, 'n'), Grad(22, 69, 'o'), Grad(75, 55, 'p'),
Grad(51, 68, 'r'), Grad(91, 23, 's'), Grad(22, 42, 't'), Grad(47, 80, 'u'),
Grad(60, 10, 'w'), Grad(91, 79, 'v'), Grad(5, 66, 'x'), Grad(42, 90, 'y'),
Grad(23, 59, '1'), Grad(46, 83, '2'), Grad(93, 63, '3'), Grad(47, 17, '4'),
Grad(53, 79, '5'), Grad(76, 23, '6'), Grad(91, 62, '7'), Grad(44, 97, '8'),]

proizvoljno_resenje = Resenje(gradovi)

simulator = SimulatorKaljenja(proizvoljno_resenje)

t0 = time.time()
simulator.start(ponavljanja=10, alpha=0.999, pocetna_temperatura=90)
print "Izvrsenje trajalo", time.time() - t0, "sekundi"

simulator.stampaj_resenje()


if __name__ == "__main__":
main()

I jedan od rezultata može se prikazati graficima (rezultat metode stampaj_resenje):