Naloga

Imamo koordinate nekaterih slovenskih krajev.

koordinate = [('Piran', (0, 0)), ('Koper', (8, 2)),
              ('Ilirska Bistrica', (49, 5)), ('Postojna', (46, 29)),
              ('Nova Gorica', (2, 48)), ('Ajdovščina', (24, 42)),
              ('Idrija', (34, 54)), ('Logatec', (48, 46)),
              ('Cerknica', (60, 31)), ('Vrhnika', (57, 51)),
              ('Žiri', (39, 60)), ('Ljubljana', (68, 61)),
              ('Ribnica', (87, 26)), ('Kočevje', (95, 15)),
              ('Grosuplje', (82, 49)), ('Litija', (95, 61)),
              ('Kranj', (58, 82)), ('Kamnik', (78, 80)),
              ('Škofja Loka', (54, 73)), ('Trbovlje', (112, 71)),
              ('Novo mesto', (119, 32)), ('Krško', (162, 56)),
              ('Celje', (129, 80)), ('Maribor', (156, 117)),
              ('Velenje', (117, 94)), ('Slovenska Bistrica', (150, 97)),
              ('Murska Sobota', (196, 138)), ('Ptuj', (173, 102)),
              ('Ormož', (196, 100)), ('Ljutomer', (199, 112)),
              ('Gornja Radgona', (184, 139)),
              ('Jesenice', (36, 102)),
              ("Kot pri Dramlju", (119, -6))]

(H krajem, ki so bili sicer v domači nalogi sem tu dodal še dva. Jesenice: Slovenija brez Jesenic je videti čudna, razpotegnjena navzgor. Jesenice potrebujemo, da jo "napnejo" čez celo površino. Kot pri Dramlju: najjužnejši slovenski kraj. V originalnem seznamu je bil Piran hkrati najjužnejši in najzahodnejši, kar olajša nalogo, zato rešitev ni tako nazorna, kot bi bila sicer.)

Najprej si pripravimo funkcijo, ki nariše zemljevid, da bomo lahko gledali, kaj delamo. Uporabili bomo matplotlib. Zemljevid je razsevni diagram (plt.scatter, točke označimo s plt.annotate).

%matplotlib inline

import numpy as np
from matplotlib import pyplot as plt

def zemljevid():
    vid = plt.gcf()
    vid.set_size_inches(10, 7)   
    
    x = [koord[0] for _, koord in koordinate]
    y = [koord[1] for _, koord in koordinate]
    # Šlo bi tudi v enem zamahu, tako:
    # x, y = zip(*(koord for _, koord in koordinate))

    plt.scatter(x, y)

    for ime, (x, y) in koordinate:
        plt.annotate(ime, (x,y), textcoords="offset points", xytext=(1,5), ha='center')

zemljevid()

Očitnejša rešitev: preudarno dodajanje

Začeli bomo z najzahodnejšim krajem in si želeli iti navzgor. Ker to ne bo šlo, bomo zavili v kraj, ki bo najmanj odstopal od smeri navzgor (to bo Nova Gorica). V naslednjem koraku bomo morali zaviti nekoliko desno; želimo si čim manj desno, torej bomo izbrali tisti kraj, ki bo najbolj v isti smeri kot je bila smer Piran - Nova Gorica. (Nevarnosti, da bi zavili levo, ni, saj je Piran najbolj levi kraj in Nova Gorica najbolj leva možna smer.) Naslednji kraj bodo torej Jesenice. Zdaj nadaljujemo v istem duhu: spet bomo morali zaviti malo desno, vendar bi radi zavili čim manj. Izbrali bomo kraj, ki bo najbolj v isto smer kot je bila smer Nova Gorica - Jesenice; to nas bo pripeljalo naravno v Gornjo Radgono. Iz nje nadaljujemo v smer, ki je čimbolj podobna smeri Jesenice - Gornja Radgona in se znajdemu v Murski Soboti...

Zdaj, ko smo se domislili splošne strategije, moramo rešiti matematični problem: če imamo, recimo, smer Nova Gorica - Jesenice: kako med vsemi možnimi nadaljevanji izberemo tisto, ki je najbolj podobno tej smeri? Če si povezave med kraji predstavljamo kot vektorje, je "podobnost smeri" kar kot med vektorjema. Spomnimo se kosinus kota je enak normiranemu skalarnemu produktu:

$$cos(\phi) = \frac{\mathbf{a}\mathbf{b}}{|\mathbf{a}||\mathbf{b}|}$$

Vsi naši skalarni produkti bodo pozitivni, vsi koti manjši od 180 stopinj ... zato se ne bomo zafrkavali s kosinusom, temveč rekli kar, da sta si dve smeri tem bolj podobni, čim večji je normirani skalarni produkt med njima.

Pripravimo si funkcijo skalp, ki bo računala skalarni produkt dveh vektorjev, podanih s komponentami. Funkcija bo predpostavila, da je prvi vektor normiran, tako da bo vsoto produktov komponent delila samo z dolžino drugega vektorja. Zgodilo se bo tudi, da bo dolžina drugega vektorja enaka 0; v tem primeru bo delila z 1, da ne dobimo deljenja z 0.

from math import sqrt

def skalp(ax, ay, bx, by):
    return (ax * bx + ay * by) / (sqrt(bx ** 2 + by ** 2) or 1)

Zdaj pa imamo vse kar potrebujemo. Najprej izberemo prvi kraj. Ker smo se na prejšnjih predavanjih mimogrede naučili nekaj o lambdah, vemo, da lahko pokličemo kar min(koordinate), zraven pa kot ključ podamo funkcijo, ki ji bo min podajala elemente seznama koordinate (terke) in funkcija bo vrnila element z indeksom 1, torej dejanske koordinate.

V spremenljivkah kraj, tx in ty bomo hranili trenutni kraj in njegove koordinate, px in py pa bosta želena smer - v začetku je to navzgor. elastika bo vsebovala seznam krajev, ki se jih elastika dotakne.

Pred zanko damo v elastiko prvi kraj. Nato vrtimo zanko, dokler se ne vrnemo v izhodiščni kraj, torej dokler je elastika[0] različna od elastika[-1], pri čemer seveda popazimo na posebni primer v začetku, ko elastika vsebuje en sam kraj. (Morda je kdo opazil: če bi Python imel zanko do - while, bi lahko prestavili pogoj na konec, pa ne bi potrebovali len(elastika) == -1. Kot sem rekel na drugih predavanjih: Python takšne zanke nima, vendar nas to ne boli preveč, saj je možno tudi v situacijah, ko bi jo res potrebovali, brez večjih muk shajati brez nje.)

Znotraj zanke poiščemo tisti kraj in njegove koordinate, za katerega velja, da je smer do njega čimbolj poravnana z zadnjo prepotovano smerjo. Zadnja prepotovana smer je px, py, novo smer pa dobimo tako, da od koordinate posamičnega kraja odštejemo trenutne koordinate tx, ty. Zdaj bomo uporabili max(koordinate), ključ pa je lambda, ki prejema argumente k, ki so terke (kraj, (x, y)). Do koordinat torej pridemo s k[1][0] in k[1][1], kar ni najbolj pregledno, je pa žal neizogibno, če hočemo tole rešiti z lambdo.

Ko imamo naslednji kraj (kraj) in njegove koordinate (nx, ny), izračunamo novo trenutno smer (px, py), shranimo nove koordinate kot trenutne koordinate in dodamo kraj v elastiko.

To je vse:

kraj, (tx, ty) = min(koordinate, key=lambda t: t[1])
px, py = 0, 1

elastika = [kraj]
while len(elastika) == 1 or elastika[0] != elastika[-1]:
    kraj, (nx, ny) = max(koordinate, key=lambda k: skalp(px, py, k[1][0] - tx, k[1][1] - ty))
    px, py = nx - tx, ny - ty
    tx, ty = nx, ny
    elastika.append(kraj)

(Je kdo opazil, da vektor (px, py) v resnici ni normiran? Skalarni produkt pa misli, da bo? Nič hudega. Vsi skalarni produkti so zato preveliki, vendar so preveliki za enak faktor. Ker jih uporabljamo za urejanje kandidatov, zaradi tega ne bo nič narobe.)

Narišimo, da vidimo, da je rešitev pravilna.

Mimogrede strpamo koordinate v slovar, da bomo elastiko narisali tako, da bomo nabrali vse koordinate x in y za vsak kraj v elastika.

zemljevid()

koord = dict(koordinate)
plt.plot([koord[ime][0] for ime in elastika], [koord[ime][1] for ime in elastika])
[<matplotlib.lines.Line2D at 0x7f9800097040>]

Razmislek o rešitvi

Rešitev je kar preprosta in očitna, ni pa hitra. Recimo, da elastika vsebuje k krajev, vseh krajev pa je n. Rešitev tedaj izračuna kn skalarnih produktov. V najslabšem primeru so v elastiki vsi kraji; v tem primeru rešitev izračuna n2 skalarnih produktov. Če vzamemo desetkrat več krajev, bo trajalo stokrat dlje. To ni najhujša katastrofa, ki se lahko zgodi, obstajajo tudi hujši problemi s počasnejšimi rešitvami. Vseeno pa raje naredimo hitreje, če gre. In gre.

Hitrejša rešitev: postopno odstranjevanje

Spet izberimo najvzhodnejši kraj; če bi bila takšna kraja dva, bomo vzeli najjužnejšega. Vse kraj uredimo po kotu, ki ga oklepajo z navpično. Z drugimi besedami, razvrstimo jih v smeri urinega kazalca, ki bi bil pripet v najvzhodnejši kraj. Na začetek in konec prilepimo še ta kraj.

kraj0, (x0, y0) = el0 = min(koordinate, key=lambda t: t[1])

elastika = [el0] \
           + sorted(koordinate, key=lambda k: (k[1][1] - y0) / sqrt((k[1][0] - x0)** 2 + (k[1][1] - y0) ** 2 or 1)) \
           + [el0]

Ključ za funkcijo sorted je skalarni produkt z navpičnico.

Zakaj to počnemo, bomo videli, če pogledamo, kaj smo dobili.

zemljevid()

plt.plot([koord[ime][0] for ime, _ in elastika], [koord[ime][1] for ime, _ in elastika])
[<matplotlib.lines.Line2D at 0x7fb0020ef910>]

Namesto najvzhodnejšega bi lahko vzeli tudi najjužnejši kraj in kraje razvrščali glede na kot z vodoravnico.

kraj0, (x0, y0) = el0 = min(koordinate, key=lambda t: t[1][::-1])

elastika = [el0] \
           + sorted(koordinate, key=lambda k: (k[1][0] - x0) / sqrt((k[1][0] - x0)** 2 + (k[1][1] - y0) ** 2 or 1)) \
           + [el0]

zemljevid()

plt.plot([koord[ime][0] for ime, _ in elastika], [koord[ime][1] for ime, _ in elastika])
[<matplotlib.lines.Line2D at 0x7fb002287e50>]

Pravzaprav bi lahko vzeli poljuben kraj ali, recimo, težišče krajev ali kako drugo točko in razvrščali glede na kote. Vendar ostanimo pri najbolj južnem, ker je to najbolj praktično, pa tudi slika je kar lepa.

Kar smo dobili, je neka še kar pametno napeta elastika, ki pa vsebuje veliko odvečnih krajev.

Če gremo po vrsti po teh krajih - Piran, Koper, Ilirska Bistrica, Nova Gorica ... vemo, da bomo vedno zavijali samo desno, nikoli levo. To bo recept za odstranjevanje krajev.

Začnemo s Piranom.

In tako naprej.

Kako vemo, ali je ovinek levi ali desni? Prej smo imeli skalarni produkt, za določanje orientacije potrebujemo vektorskega. Vedno bomo opazovali zadnje tri točke; njihove koordinate naj bodo (ax, ay), (bx, by) in (cx, cy). Vektor iz prve v drugo je (bx - ax, by - ay) in vektor iz druge v tretjo je (cx - bx, cy - by). Zanima nas, v kakšno smer je potrebno zavrteti prvi vektor, da dobimo drugega. Če je predznak vektorskega produkta pozitiven, gre za vrtenje nasprotno smeri urinega kazalca (counter clockwise, ccw) - to je levi ovinek in to je tisto, česar ne maramo.

Sprogramirajmo funkcijo, ki vrne True, če so točke a, b, c razporejene tako, da je ovinek, ki ga naredimo v b, levi.

def ccw(ax, ay, bx, by, cx, cy):
    return (bx - ax) * (cy - by) - (cx - bx) * (by - ay) >= 0

Vse je pripravljeno. Sprogramirajmo. Dodali bomo še nekaj print-ov, da bo lažje slediti.

i bo indeks v elastiko. Vse točke do i-te so "pogojno sprejete*. V vsakem krogu zanke povečamo i, s čimer sprejmemo novo točko. Potem kličemo ccw s koordinatami točk z indeksi i, i - 1 in i - 2. Če tvorijo levi ovinek, odstranimo točko i - 1; zaradi praktičnosti to storimo tako, da zmanjšamo i in pobrišemo i-to. Tudi odstranjevanje je v zanki, saj lahko v istem koraku odstranimo več točk.

Funkciji ccw moramo poslati šest argumentov, ki jih dobimo v treh terkah. Klic poenostavimo tako, da damo pred terko zvezdico in Python bo vse elemente terke uporabil kot ločene argumente.

kraj0, (x0, y0) = el0 = min(koordinate, key=lambda t: t[1][::-1])

elastika = [el0] \
           + sorted(koordinate, key=lambda k: (k[1][0] - x0) / sqrt((k[1][0] - x0)** 2 + (k[1][1] - y0) ** 2 or 1)) \
           + [el0]

print("Imam", elastika[0][0])
i = 1
while i < len(elastika):
    print("Dodal sem", elastika[i][0])
    while i > 1 and ccw(*elastika[i - 2][1], *elastika[i - 1][1], *elastika[i][1]):
        i -= 1
        print("Odstranim", elastika[i][0])
        del elastika[i]
    i += 1
Imam Kot pri Dramlju
Dodal sem Piran
Dodal sem Koper
Dodal sem Ilirska Bistrica
Dodal sem Nova Gorica
Odstranim Ilirska Bistrica
Odstranim Koper
Dodal sem Postojna
Dodal sem Ajdovščina
Odstranim Postojna
Dodal sem Cerknica
Dodal sem Vrhnika
Odstranim Cerknica
Odstranim Ajdovščina
Dodal sem Idrija
Odstranim Vrhnika
Dodal sem Logatec
Dodal sem Žiri
Odstranim Logatec
Odstranim Idrija
Dodal sem Kočevje
Dodal sem Ribnica
Odstranim Kočevje
Dodal sem Škofja Loka
Odstranim Ribnica
Odstranim Žiri
Dodal sem Jesenice
Odstranim Škofja Loka
Dodal sem Ljubljana
Dodal sem Kranj
Odstranim Ljubljana
Dodal sem Grosuplje
Dodal sem Kamnik
Odstranim Grosuplje
Odstranim Kranj
Dodal sem Litija
Dodal sem Trbovlje
Odstranim Litija
Odstranim Kamnik
Dodal sem Velenje
Odstranim Trbovlje
Dodal sem Novo mesto
Dodal sem Kot pri Dramlju
Dodal sem Celje
Odstranim Kot pri Dramlju
Odstranim Novo mesto
Dodal sem Maribor
Odstranim Celje
Odstranim Velenje
Dodal sem Slovenska Bistrica
Dodal sem Gornja Radgona
Odstranim Slovenska Bistrica
Odstranim Maribor
Dodal sem Ptuj
Dodal sem Murska Sobota
Odstranim Ptuj
Dodal sem Ljutomer
Dodal sem Krško
Dodal sem Ormož
Odstranim Krško
Dodal sem Kot pri Dramlju

Pa še brez print-ov, da bo očitneje, kako preprosto je vse skupaj (ko enkrat veš, kaj bi rad sprogramiral).

kraj0, (x0, y0) = el0 = min(koordinate, key=lambda t: t[1][::-1])

elastika = [el0] \
           + sorted(koordinate, key=lambda k: (k[1][0] - x0) / sqrt((k[1][0] - x0)** 2 + (k[1][1] - y0) ** 2 or 1)) \
           + [el0]

i = 1
while i < len(elastika):
    while i > 1 and ccw(*elastika[i - 2][1], *elastika[i - 1][1], *elastika[i][1]):
        i -= 1
        del elastika[i]
    i += 1

Narišimo, da vidimo, da res deluje.

zemljevid()

plt.plot([koord[ime][0] for ime, _ in elastika], [koord[ime][1] for ime, _ in elastika])
[<matplotlib.lines.Line2D at 0x7fb0022733d0>]

Razmislek o rešitvi

Tale rešitev je bistveno hitrejša. Najpočasnejši del je urejanje: n točk je mogoče urediti v času sorazmernem nlog n, kar je bistveno manj od n2.

Sledita jim dvojni zanki, kar bi na prvi pogled lahko imelo zahtevnost n2, vendar ni tako: vsaka točka je lahko odstranjena le enkrat. Torej se bo notranja zanka skupno, v celotnem poteku algoritma, zavrtela največ tolikokrat, kolikor je točk.

Rešitve, ki so jih pošiljali študenti, so bili pogosto neka variacija te rešitve. Nekateri so napisali praktično enako rešitev, le na kak daljši način in z drugo začetno točko; tipično je bila to neka točka v sredini, na primer težišče. Drug pogost pristop je bilo, da so vzeli najbolj levi in najbolj desni kraj, ostale pa razdelili na tiste, ki so nad in pod povezavo med njima. Kraje so uredili po koordinati x in nato potovali po zgornji strani in izločali kraje, ki so bili med dvema višjima sosedoma, nato pa po spodnji strani in izločali kraje med dvema nižjima. Še druge variacije rešitve so razdelile zemljevid v štiri kvadrante.

Vse to je dobro in je na poti do te, optimalne rešitve.

Razmislek o nalogi - kaj smo počeli?

Problem, ki smo ga reševali, se imenuje iskanje konveksna ogrinjača (convex hull). Le-ta je, skupaj z nekaj sorodnimi problemi, pomembna reč v matematiki in računski geometriji. Konveksno ovojnice lahko iščemo tudi v tridimenzionalnem ali v kakem višje- dimenzionalnem prostoru, kjer so stvari seveda bolj zapletene. Takšnole preprosto urejanje po kotu si lahko privoščimo le na ravnini. Algoritem, ki smo ga sestavili, je Grahamovo preiskovanje.

Funkcijo, ki smo jo poimenovali ccw, običajno imenujemo ccw in je vlečni konj računske geometrije. Za primer: zanima nas, ali se daljici AB in CD sekata. Najzanesljivejši način, da to preverimo je tale: orientaciji (ccw) trojke ABC in ABD morata biti različni, prav tako morata biti različni orientaciji trojk CDA in CDB. Z drugimi besedami, C in D morata biti na različnih straneh AB, A in B pa na različnih straneh CD.