Ferenci Tamás (https://www.medstat.hu/)
- Összefoglaló
- Kérdésfelvetés
- Módszertan
- Eredmények
- Kis történelmi kitérő az izokrón térképek kapcsán
- A legtávolibb magyar települések
- Az elérési idők eloszlása
- Az egész ország átlagos közelsége egy településről
- Az ország középpontja, avagy hová telepítsünk kórházat?
- A magyar TSP megoldása: hogy lehet leggyorsabban körbejárni az összes magyar települést?
- Továbbfejlesztési ötletek, kutatási lehetőségek
- Az útvonaltervezés technikai részletei
- A kórháztelepítés megoldása lineáris programozással
- Számítógéppel letöltöttem Magyarország térképét, majd egy útvonaltervező programmal meghatároztam a legrövidebb utat bármely két magyar település között, közúton, autóval haladva. Ez az adatbázis számos érdekes kérdés vizsgálatát teszi lehetővé.
- Az egyik legkézenfekvőbb, hogy megnézhetjük, hogy adott településről mennyi idő alatt érhető el az ország többi része. Ebből szép színes térkép is gyártható, illetve megjeleníthetjük azt is, hogy hol van az a határ, amit adott időn belül (pl. 1 órán belül, 2 órán belül stb.) el lehet érni; ezt szokták izokrón térképnek nevezni. Ezt is kirajzoltam az összes magyar településre.
- Kiderül írásomból, hogy melyik a két „legtávolabbi” magyar település közúton (spoiler: Kishódos és Felsőszölnök, a forgalom nélküli, de szabályokat betartó közlekedéssel az elvi eljutási idő majdnem 7,5 óra).
- Kiderül, hogy mi az „ország középpontja”, azaz az a település, ahonnan az összes többi átlagosan a legrövidebb idő alatt érhető el. Ezt nézhetjük lélekszámmal súlyozva is, mondván, hogy a nagyobb településekre való eljutás ideje nagyobb súllyal essen latba. De akár azt is megkérdezhetjük, hogy melyik település az, amire igaz, hogy a tőle legtávolabbi település távolsága a lehető legkisebb…?
- A dolgot tovább is vihetjük: a „középpont” felfogható úgy is, mint ahová mondjuk egy kórházat érdemes telepíteni, ha egyetlen kórház van az országban. De hová rakjuk a kórházakat (továbbra is az eljutási időt minimalizálva), ha 2 van? Ha 3? 10? 50? Megmutatom e feladat megoldását is, melyhez operációkutatási eszközöket fogunk bevetni; itt hosszabban beszélek ezen eredmény limitációiról is.
- Ezek a vizsgálatok elvégezhetőek nem csak Magyarország településeire, hanem egy településen belül is; ezt egy egy másik írásom mutatja be Budapest példáján.
- Az összes számításomat és eredményemet teljesen transzparensen közlöm: mind a nyers adatokat megadom, mind a kódokat valamennyi számítás mögött, melyek együtt lehetővé teszik munkám teljes reprodukálását és ellenőrzését. A dolgozat végén több továbbfejlesztési lehetőségre, további érdekes kutatási irányra is utalok.
Sokan szinte nap mint nap használnak valamilyen útvonaltervező programot. Ilyenekkel pillanatok alatt megkereshetjük, hogy melyik a legrövidebb vagy leggyorsabb út két pont (legegyszerűbb esetben: két település) között, de ezt látva felvetődhet az emberben: akkor már miért nem bővítjük a vizsgálódásunk körét? Mi volna, ha minden magyar település-pár között megkeresnénk a legrövidebb utat…? Nem is azért, hogy a navigációt javítsuk, arra tökéletesek ezek az alkalmazások, hanem azért, hogy a rendszerről, jelen esetben a magyar közúti rendszerről megtudjunk valamit! Hogyan épül fel? Mely települések érhetőek el gyorsabban, melyek lassabban? Adott településből hova lehet 1 órán belül eljutni? Igaz, hogy Budapest-centrikus az útvonalhálózat? Mennyire erős ez a hatás? Az ember eleresztheti jobban is a fantáziáját. Melyik két magyar település van a legtávolabb – akár kilométerben, akár időben – egymástól közúton? Mennyi a minimális út vagy idő, ami ahhoz szükséges, hogy az összes magyar települést autóval felkeressük? (És mi az optimális útvonal?) Melyik az a település, amire igaz, hogy a tőle vett távolságok átlaga a legkisebb, azaz, ahonnan az átlagosan legkönnyebb bárhová eljutni az országba? Na és amire igaz, hogy a tőle legtávolabbi település távolsága a legkisebb? (A dolog persze szimmetrikus: akkor ide a legkönnyebb eljutni bárhonnan. Hova érdemes kórházat telepíteni, ha azt szeretnénk, hogy a leggyorsabban elérhető legyen, átlagosan, vagy legrosszabb esetben is?)
Ezekre a kérdésekre fogok a jelen esszében megpróbálni válaszolni.
(Természetesen ezek a kérdések nem csak a különböző települések közötti eljutás kapcsán tehetőek fel, hanem egy település különböző pontjai között eljutásra vonatkozóan is. Egy másik írásom ezt mutatja be Budapest példáján.)
A problémakört szép nevén közúti elérhetőségnek lehetne nevezni. Valójában nincs szó nagy újdonságról, sokaknak ismerős lehet a régi autóstérképekről a kis táblázat, mely a nagy városok között tüntette fel – táblázat alakban elrendezve az összes lehetséges pár között – a távolságot. Ugyanezt fogjuk most megtenni, csak kihasználva az informatika lehetőségeit sokkal alaposabban és átfogóbban, valamint a kapott eredményeknek többféle alkalmazását is ki fogjuk próbálni.
Több útvonaltervező alkalmazás van széles körű használatban, azonban nem mindegyik támogatja (pláne könnyen) a nagyon sok település-pár közötti tömeges lekérdezést, amire most szükségünk volna. Nem csak arról van szó, hogy ez technikailag hogy oldható meg, hanem arról is, hogy ha interneten kell a kéréseket átküldeni, az könnyen lehet, hogy nagyon lassú lesz. Mindezekre tekintettel a választásom az OpenStreetMap rendszerre esett. Ez több fontos előnnyel is bír:
- Ingyenes és nagyon szabad licensszel bír.
- Az egész ország térképe letölthető (erre a célra a Geofabrik oldalát használtam), így minden számítás elvégezhető off-line, a helyi számítógépen, internet szüksége nélkül.
- Elérhető hozzá olyan program, amely nagyon jó tulajdonságú és gyors útvonaltervezést végez a térkép alapján, ez az OSRM (Open Source Routing Machine), ami ráadásul ingyenesen és nyílt forráskódú.
A fenti rendszer felállításának és az útvonaltervezés végrehajtásának a technikai részleteit külön pontban tárgyalom meg.
Néhány megszorítást tegyünk. Az első, hogy most kizárólag közúti eljutással fogunk foglalkozni. Ugyanígy feltehető azonban a kérdés bármilyen más közlekedési módra, legyen az bicikli, gyaloglás, vagy akár tömegközlekedés. Ezek vizsgálata továbbfejlesztési lehetőség; annyit jegyzek meg, hogy az előbbiek vizsgálata nem nehéz, hiszen pontosan ugyanaz a térkép, és ugyanaz az útvonaltervező használható (lényegében csak az ún. profilt kell lecserélni az OSRM-ben), az utóbbi azonban zűrösebb, hiszen a menetrendre is szükség van a vizsgálatához.
A második fontos megjegyzés, hogy a számítások eredménye függhet attól, hogy mikor futtatjuk le. Ennek van egy kellemetlen, és van egy nagyon is érdekes aspektusa. Az utóbbi a hosszú távú trend: továbbfejlesztési lehetőségként felvethető, hogy ha mondjuk különböző években nézzük meg az eredményeket, akkor vajon mutatkozik-e javulás Magyarország közúti elérhetőségében? Ha igen, mekkora? Mikor volt nagyobb a javulás, mikor volt kisebb? Mely területeken volt nagyobb javulás és mely területeken kisebb? Ez tehát nem csak, hogy nem probléma, hanem ellenkezőleg, egy potenciálisan nagyon izgalmas kutatási kérdés, a másik vetület azonban káros: a rövid távú hatások. Egy útlezárás, felújítási munka stb. megváltoztathatja a menetidőket, adott esetben nagyon rövid időre, és ebből fakadóan minden érdemi jelentőség nélkül. Ennél sokkal fontosabb a dugók kérdése: ez nem kiküszöbölendő, hanem egy beépítendő tényező, hiszen ezeknek igenis van érdemi jelentősége, pláne, hogy ebben vannak állandó mintázatok: melyik úton van tipikusan dugó, melyiken nincs, mi a helyzet nappal, éjszaka, csúcsidőben stb. Ez tehát egy nagyon is vizsgálatra érdemes szempont, van-e ahol az állandó dugók miatt lassabb eljutással kell számolni, hiába rövidebb az út, mennyire más a helyzet éjszaka indulva mint csúcsidőben és így tovább. A probléma, hogy bár az OSRM ma már támogatja a forgalmi adatok felhasználását és figyelembevételét, én nem találtam adatforrást hozzá, ahonnan ezeket az információkat nyílt forrásból le lehetne tölteni. Ezt tehát továbbfejlesztési lehetőségként nyitva hagyom, addig is nagyon fontos rögzíteni, hogy a mostani eredmények konkrét forgalomtól független, átlagos sebességgel számolt értékek!
A harmadik szűkítés, hogy jelen esszében kizárólag eljutási idővel fogok foglalkozni, távolsággal nem. Ennek egy rettenetesen földhözragadt oka van, az, hogy az OSRM alapból időre optimalizál, és csak kézi barkácsolás árán lehet ezen változtatni. Ebbe most nem mentem bele, továbbfejlesztési lehetőségként hagyom nyitva ezt a kérdést.
Minden számítást az R statisztikai környezet alatt végeztem el, a kódokat – a transzparencia és a nyílt tudomány jegyében – ez a dokumentum teljes egészében tartalmazza. Ezeket az egyes elemzéseknél külön-külön is közölni fogom, ha valakit ez nem érdekel, a szürke hátteres részeket nyugodtan ugorja át.
Ha megtervezzük az útvonalat az összes lehetséges településpár között, akkor lényegében egy mátrixot kapunk, matematikai szóval élve. Ezeknek csakugyan van szép matematikája, de a mostani céljainkra túlnyomó többségében az is elég lesz, ha egyszerűen számok négyzet alakban elrendezett táblázataként gondolunk rá, ahol egy adott szám megadja, hogy a sorában szereplő településből mennyi idő alatt lehet eljutni az oszlopában szereplő településhez. Pont mint az autóstérképek végében szokott lenni. Olvassuk is be az adatokat, konvertáljuk az időtartamokat órába, mert eredetileg másodpercben van, az OSRM így szolgáltatja, és nézzük is meg egy kis részletét ennek a mátrixnak (táblázatnak):
durations <- fread("osrmdurations.csv")
durations <- as.matrix(durations)
rownames(durations) <- colnames(durations) <- locs$NAME
durations <- durations/60/60
durationsSymm <- (durations + t(durations))/2
kableExtra::column_spec(knitr::kable(durations[1:5, 1:5], format = "html"), 1, bold = TRUE)| Aba | Abádszalók | Abaliget | Abasár | Abaújalpár | |
|---|---|---|---|---|---|
| Aba | 0.000000 | 2.918028 | 2.219889 | 2.084722 | 3.471139 |
| Abádszalók | 2.872611 | 0.000000 | 4.668250 | 1.551722 | 2.243694 |
| Abaliget | 2.215917 | 4.693556 | 0.000000 | 3.940139 | 5.326556 |
| Abasár | 2.079500 | 1.548000 | 3.931167 | 0.000000 | 1.980611 |
| Abaújalpár | 3.468972 | 2.248750 | 5.320639 | 1.975000 | 0.000000 |
Ez természetesen csak egy pici részlet; az egész mátrixnak 3177 sora, és persze ugyanennyi oszlopa van.
Az rögtön látható, és értelemszerű, hogy a főátlóban csupa 0 elem található. A második észrevétel, hogy a mátrixunk nem tökéletesen szimmetrikus. Eléggé az, de nem pontosan. Mindazonáltal a nagyobb eltérések ritkák, és ahol elő is fordulnak, általában ott is valamilyen nagyon partikuláris ok van a háttérben. (Példának okáért vegyük a legnagyobb különbséget: Lórévből Rácalmás 24 perc, de Rácalmásról Lórév 59 perc. Megnézve azonban a belinkelt térképen a konkrét útvonalat, azonnal kiderül, hogy valamiért a tervező-program azt hiszi, hogy a Lórév és Adony közötti kompon csak az egyik irányban lehet átkelni, ami elég furcsa lenne, és persze nincs is így.) Bizonyos célokra célszerűbb lehet szimmetrizálni a mátrixot (kiátlagolni a két értéket; mátrixos nyelven szólva: összeadva a mátrixot a transzponáltjával és elosztva 2-vel); az előbbiekből látható, hogy nagy módosítást ez nem jelent, most ezt is kiszámoljuk.
A későbbi feldolgozáshoz célszerű lehet, ha átrendezzük egy olyan formátumba az adatokat, melyben minden sor egy pár (adattudósok úgy mondanák: long formátum). Ha ezt megtesszük, akkor ráadásul a két típusú adat kényelmesen tárolható egyetlen objektumban:
durationsLong <- as.data.table(reshape2::melt(durations, value.name = "Duration"))
durationsLong$Var1 <- as.character(levels(durationsLong$Var1))[durationsLong$Var1]
durationsLong$Var2 <- as.character(levels(durationsLong$Var2))[durationsLong$Var2]
durationsSymmLong <- as.data.table(reshape2::melt(durationsSymm, value.name = "Duration"))
durationsSymmLong$Var1 <- as.character(levels(durationsSymmLong$Var1))[durationsSymmLong$Var1]
durationsSymmLong$Var2 <- as.character(levels(durationsSymmLong$Var2))[durationsSymmLong$Var2]
saveRDS(durationsSymmLong, "durationsSymmLong.rds")
fwrite(durationsLong, "durationsLong.csv", dec = ",", sep = ";", bom = TRUE)
if(!file.exists("durationsLong.zip")) zip("durationsLong.zip", "durationsLong.csv")
knitr::kable(head(durationsLong))| Var1 | Var2 | Duration |
|---|---|---|
| Aba | Aba | 0.000000 |
| Abádszalók | Aba | 2.872611 |
| Abaliget | Aba | 2.215917 |
| Abasár | Aba | 2.079500 |
| Abaújalpár | Aba | 3.468972 |
| Abaújkér | Aba | 3.362333 |
Az adatállomány ezen a linken külön is elérhetővé tettem tömörített CSV formátumban, hogy bárki tetszőleges saját elemzéshez is felhasználhassa.
Most már nekiállhatunk a munkának! Első feladatként nézzük meg, hogy milyen messze vannak (menetidőben) a magyar települések egy bizonyos ponttól; legyen ez Budapest, 1. kerület (ha már a 0 kilométerkő itt van):
knitr::kable(head(durationsLong[Var1=="Budapest 01. kerület"]))| Var1 | Var2 | Duration |
|---|---|---|
| Budapest 01. kerület | Aba | 0.9066667 |
| Budapest 01. kerület | Abádszalók | 2.2096389 |
| Budapest 01. kerület | Abaliget | 2.7583333 |
| Budapest 01. kerület | Abasár | 1.2558333 |
| Budapest 01. kerület | Abaújalpár | 2.6422500 |
| Budapest 01. kerület | Abaújkér | 2.5355000 |
És most jön a lényeg. Hogyan ábrázoljuk ezt? Célszerűen térképen! Egyszerűen színezzük be az összes települést az itt látható idővel arányosan, így egy látványos, azonnal áttekinthető ábrává konvertáljuk ezt a nagyon hosszú adathalmazt:
ggplot(merge(geodata, merge(locs, durationsLong[Var1=="Budapest 01. kerület",
.(NAME = Var2, Duration)])),
aes(x = X, y = Y, fill = Duration)) + geom_sf(color = NA) +
labs(x = "", y = "", fill = "Eljutási idő [h]", caption = captionlab)
Az egyetlen probléma, hogy élesen látszódnak az egyes települések körvonalai. (A térkép a településeket nem pontként jeleníti meg, hanem egy kis területet rendel hozzájuk úgy, hogy átfedésmentesek legyenek, és együtt kiadják az egész országot.) Ez egyáltalán nem hibás, de pusztán esztétikai értelemben zavaró lehet, ha mi inkább egy összképet akarunk látni. Azt kellene ehhez elérni, hogy ne legyenek ilyen éles ugrások. Egy lehetséges megoldás a simítás, jobban mondva jelen esetben az interpoláció: azt elérni, hogy az értékek ne hirtelen ugorjanak át a két település határán, hanem folytonos legyen az átmenet. Ehhez a településeket most egyetlen pontnak fogjuk tekinteni (ami nem más lesz, mint a kis területek középpontja), ahhoz rendeljük az értékeket, és azok között fogunk lineárisan interpolálni. Illusztrálva mindezt, azt mondhatjuk, hogy e helyett a színezés helyett (a pontok a városok, a függőleges tengely a szín intenzitása):
interppelda <- data.frame(x = 1:5, y = c(3, 1, 2, 5, 0))
ggplot(interppelda, aes(x = x, y = y)) + geom_point(size = 3) + geom_step(direction = "mid")
A következő színezést használjuk:
ggplot(interppelda, aes(x = x, y = y)) + geom_point(size = 3) + geom_line()
A mi konkrét helyzetünk két szempontból nehezebb mint a fenti
illusztráció: egyrészt nem egy, hanem két dimenzióban vagyunk, tehát a
pontok a síkban helyezkednek el, és minden irányban lehetnek
szomszédaik, másrészt nem szép szabályos rácspontokban vannak, hanem
össze-vissza. Mindezek kezelésére vannak bevált algoritmusok, mi most
Akima módszerét fogjuk
használni; ezt az R-ben az akima csomag valósítja meg. (Interpoláció
helyett használhatnánk ügyesebb simítást is, például spline-okat, de
jelen esetben nincs sok értelme: a pontok nagyon sűrűn vannak, így egy
ilyen simítás csak jóval lassabb lenne, miközben nem igazán nyernénk
vele semmit.) Az akima::interp végrehajtja az interpolációt, majd az
eredményeket egyenletes rácspontokban adja vissza, melyek lefedik az
ország területét. Ez utóbbi miatt még egy lépésre szükség van: ki kell
szűrnünk azokat a pontokat, amik az ország határain kívül vannak, hiszen
a fenti rácsnak lesznek ilyen pontjai is.
Ezzel kapjuk a végleges ábránkat, az ország minden településének elérhetőségét a 0 kilométerkőtől indulva:
bpcontour <- merge(durationsLong[Var1=="Budapest 01. kerület", .(NAME = Var2, Duration)], locs)[
, with(akima::interp(X, Y, Duration, nx = 1000, ny = 1000),
cbind(CJ(Y = y, X = x), Duration = c(z))[!is.na(Duration)])]
bpcontour <- bpcontour[sapply(st_intersects(st_as_sf(bpcontour, coords = c("X", "Y"),
crs = st_crs(geodata)), geodata), length)==1]
p <- ggplot(bpcontour, aes(x = X, y = Y)) + geom_raster(aes(fill = Duration)) +
labs(x = "", y = "", fill = "Eljutási idő [h]", caption = captionlab) +
metR::scale_x_longitude(ticks = 1, expand = waiver()) +
metR::scale_y_latitude(ticks = 0.5, expand = waiver())
ggsave("BudapestKozutiEljutasiIdo.pdf", p, width = 16, height = 9, device = cairo_pdf)
p
Ha egy szubjektív kiszólást tehetek, ez az ábra szerintem önmagában is szép: mint ahogy az erek hálózzák be a testet. (Ha más is így gondolná, letöltheti jó minőségű PDF formátumban is ugyanezt az ábrát!)
Mindazonáltal az ábra még tovább is fejleszthető. Vegyünk most egy másik települést példának, a Borsod-Abaúj-Zemplén megyei Csobádot. Az előbbi ábrát egészítsük ki egy ponttal, ahol a település van, de ami még fontosabb: jelöljük meg azokat a pontokat, amelyek adott időn belül, például 1 órán belül, 2 órán belül stb. elérhetőek a vizsgált településünkből. E területek határát szokás izokrón görbének nevezni (kb. azonos idő alatt elérhető): adott idő alatt eddig lehet eljutni a településből. Avagy, fordítva megfogalmazva – minthogy a dolog szimmetrikus – innen lehet adott idő alatt eljutni a településre. Technikailag egy dologra kell figyelni, ez pedig az, hogy az eredeti interpoláció nagyon finom térbeli felbontású, ami jó, hiszen ettől szép sima a térkép, viszont most nem szerencsés, mert ha ez alapján húzzuk meg ezeket a kontúrvonalakat, akkor azok nagyon szaggatottak lesznek, hiszen ide-oda fog ugrálni, hogy pontosan hol van az 1 óra határa (és a finom felbontás miatt ezt jól tudja követni). Úgyhogy ehhez érdemes egy durvább felbontású interpolációt készíteni, és az alapján behúzni az izokrón vonalakat:
contourplotter <- function(datain, geodatain, lab, telepules = NA, hascontour = FALSE, title = "",
filllims = NULL) {
telepulescontour <- datain[, with(akima::interp(X, Y, value, nx = 1000, ny = 1000),
cbind(CJ(Y = y, X = x), value = c(z))[!is.na(value)])]
telepulescontour <- telepulescontour[sapply(st_intersects(st_as_sf(telepulescontour,
coords = c("X", "Y"),
crs = st_crs(geodatain)),
geodatain), length)==1]
telepulescontour2 <- datain[, with(akima::interp(X, Y, value, nx = 40, ny = 40),
cbind(CJ(Y = y, X = x), value = c(z))[!is.na(value)])]
telepulescontour2 <- telepulescontour2[sapply(st_intersects(st_as_sf(telepulescontour2,
coords = c("X", "Y"),
crs = st_crs(geodatain)),
geodatain), length)==1]
ggplot(telepulescontour, aes(x = X, y = Y)) + geom_raster(aes(fill = value)) +
{if(!is.na(telepules)) stat_sf_coordinates(data = geodatain[geodatain$NAME==telepules,],
inherit.aes = FALSE, color = "red")} +
{if(hascontour) metR::geom_contour2(data = telepulescontour2, aes(x = X, y = Y, z = value,
label = after_stat(level)),
inherit.aes = FALSE, breaks = 1:10, skip = 0, color = "red")} +
labs(x = "", y = "", fill = lab, caption = captionlab, title = title) +
{if(!is.null(filllims)) scale_fill_gradient(limits = filllims)} +
metR::scale_x_longitude(ticks = 1, expand = waiver()) +
metR::scale_y_latitude(ticks = 0.5, expand = waiver())
}
contourplotter(merge(durationsLong[Var1=="Csobád", .(NAME = Var2, value = Duration)], locs), geodata,
"Eljutási idő [h]", "Csobád", TRUE)
Érdemes megnézni, hogy – bár jóval halványabban és elmosódottabban – de az autópályák itt is látszanak: sok település eléréséhez az a jó stratégia Csobádról indulva, ha előbb feljövünk a fővárosba.
Most, hogy a fenti térkép megvan, és minden ezzel kapcsolatos nehézséget
leküzdöttünk, gondolkozhatunk eggyel nagyobban is: miért nem csináljuk
meg ezt az összes magyar településre? A dolog innen már tényleg nem
nehéz, lényegében csak automatizálni kell a fentieket. Mindenben
követjük a fent látott algoritmust, csak egyetlen kiegészítést érdemes
tenni: célszerű egyetlen, nagy, long formátumú fájlba összefogni az
összes lehetséges paramétert, beleértve a fenti data.table
harmadik (csoportosító) argumentuma így jól kihasználható: pontosan
ugyanazt a függvényt kell meghívni, melt-eljük, és beírjuk a harmadik argumentumok közé
a csoportosításban.)
Tehát, első lépésben kikeressük adott kiinduló településhez (Var1) a
céltelepülések (Var2) koordinátáit:
contourdata <- merge(durationsLong[, .(Var1, NAME = Var2, Duration)], locs, sort = FALSE)Ezt követően megduplázzuk az adatbázist, egyszer a 40-es (durva
felbontású), egyszer pedig a finom felbontású interpolációhoz szükséges
contourdata <- rbind(cbind(contourdata, n = 40), cbind(contourdata, n = 500))Ez után jöhet az előbb említett előny, a data.table segítségével
egyetlen lépésben elintézhető az interpoláció: minden
kiindulótelepülésen és minden akima::interp eredményét
megfelelően alakítsuk át: ez alapváltozatában egy vektorban adja vissza
az CJ függvénnyel, majd a mátrixot
vektorrá alakítjuk. Ez utóbbit a c megteszi, mégpedig oszloponként,
ezért ha a CJ-ben az NA-t
fog visszaadni (lévén, hogy extrapolációt nem végzünk). Ezeket kidobjuk,
de megjegyzendő, hogy ezzel a problémával még kell foglalkoznunk, mert
így is lesznek az ország határain kívüli pontok. Végeredményben tehát a
szükséges kód:
contourdata <- contourdata[, with(akima::interp(X, Y, Duration, nx = n, ny = n),
cbind(CJ(Y = y, X = x), Duration = c(z))[!is.na(Duration)]),
.(Var1, n)]Most jöhet a határokon kívüli pontok kidobása. Mivel a fenti tábla nagyon hosszú, és rengeteg ismétlődő koordináta van benne, így nem érdemes minden sorát külön szűrni. Jobban járunk, ha kiszedjük az egyedi koordinátákat, azokat szűrjük le aszerint, hogy az ország határain belül vannak, majd az így kapott táblázatot visszaegyesítjük a nagy táblánkkal (hozva a jelölést, hogy melyik koordináta van belül). Ezt használva aztán kiszűrhetjük a szükséges pontokat, sokkal gyorsabban, mintha minden ponton végigmennénk:
uniquecoords <- unique(contourdata[, .(X, Y)])
uniquecoords$coordOK <- sapply(st_intersects(st_as_sf(uniquecoords, coords = c("X", "Y"),
crs = st_crs(geodata)), geodata), length)
contourdata <- merge(contourdata, uniquecoords, by = c("X", "Y"))
contourdata <- contourdata[coordOK==1]
contourdata <- contourdata[, .(Var1, n, X, Y, Duration)]
saveRDS(contourdata, "contourdata.rds")Innen már nincs más dolgunk, mint elvégezni a térképezést minden egyes településre. Elkészítünk egy izokrón vonalakat tartalmazó, és egy a nélküli térképet is minden településre, és az eredményeket PNG és PDF formátumokban is kimentjük, kinek mi lesz majd a célszerű:
for(loc in locs$NAME) {
p <- ggplot(contourdata[Var1==loc&n==500], aes(x = X, y = Y)) + geom_raster(aes(fill = Duration)) +
stat_sf_coordinates(data = geodata[geodata$NAME==loc,], inherit.aes = FALSE, color = "red") +
labs(x = "", y = "", fill = "Eljutási idő [h]", caption = captionlab) +
metR::scale_x_longitude(ticks = 1, expand = waiver()) +
metR::scale_y_latitude(ticks = 0.5, expand = waiver())
p1 <- p + metR::geom_contour2(data = contourdata[Var1==loc&n==40],
aes(x = X, y = Y, z = Duration, label = after_stat(level)),
inherit.aes = FALSE, breaks = 1:10, skip = 0, color = "red")
ggsave(paste0("./maps/", stringi::stri_replace_all_fixed(iconv(loc, to = "ASCII//TRANSLIT"),
c(" ", "-", "."), "",
vectorize_all = FALSE),
".png"), p, width = 16, height = 9)
ggsave(paste0("./maps/", stringi::stri_replace_all_fixed(iconv(loc, to = "ASCII//TRANSLIT"),
c(" ", "-", "."), "",
vectorize_all = FALSE),
".pdf"), p, width = 16, height = 9, device = cairo_pdf)
ggsave(paste0("./maps/", stringi::stri_replace_all_fixed(iconv(loc, to = "ASCII//TRANSLIT"),
c(" ", "-", "."), "",
vectorize_all = FALSE),
"_izokron.png"), p1, width = 16, height = 9)
ggsave(paste0("./maps/", stringi::stri_replace_all_fixed(iconv(loc, to = "ASCII//TRANSLIT"),
c(" ", "-", "."), "",
vectorize_all = FALSE),
"_izokron.pdf"), p1, width = 16, height = 9, device = cairo_pdf)
}Az eredményeket egy külön oldalon (a település nevének első betűje szerint: A-K és L-Zs) tettem elérhetővé: itt bárki letöltheti tehát minden magyar településre a fenti közúti eljutási idő térképeket!
Talán érdekes egy pillanatra kitérni az izokrón térképek történetének néhány fejezetére is. Az első ilyet Sir Francis Galton konstruálta (igen, az a Galton, lásd Galton-deszka, Galton-Watson folyamat) 1881-ben. Nagyon tanulságos megnézni, alul láthatóak az időigények, a térkép London bázissal készült:
Érdekes módon Magyarországról sem sokkal később készült izokrón térkép! Albrecht Penck német térképész 1887-ben közzétette Isochronenkarte der österreichisch-ungarischen Monarchie című cikkét, melynek térképmelléklete a vasúton történő eljutás izokrón görbéit tartalmazta a Monarchián belül, többek között Budapest (azaz Buda-Pest) központtal is.
1914-ben John G. Bartholomew publikálta izokrón térképének új kiadását (ugyanis már korábban is készített ilyet), szintén London középponttal:
A dolgot csak azért említem meg, mert egy izgalmas projekt azt célozta meg 2016-ban, azaz nagyjából 90 évvel később, hogy újraalkosson egy pontosan ugyanilyen stílusú és kinézetű térképet – csak épp az aktuális közlekedési lehetőségek mellett. Nem minden nehézség nélkül, de megoldották a feladatot, íme az eredmény:
A dolog egyébként olyan sikert aratott, hogy egy térképgyártó felkarolta, úgyhogy ha valaki szeretné, a mai napig rendelhet magának ilyen térképet a szobája falára, e pillanatban a papírra nyomott 90 x 60 cm-es térkép 70 angol fontba kerül.
Térjünk egy pillanatra vissza Magyarországra! A korszerű keresési eszközökkel, gondolok itt mindenekelőtt az Arcanumra, igazi gyöngyszemeket lehet találni.
A Budapesti Hírlap 1912. február 17-én kelt száma például beszámol arról, hogy a „magyar Adria-egyesület földrajzi szakosztályának tegnapi ülésén Prinz Gyula dr. egyetemi magántanár tartott előadást a közlekedésföldrajzi térképekről”. Ennek során az előadó „bemutatta maga szerkesztette Magyarország izokrontérképét”. Ezen a térképen „vonallal vannak összekötve az összes pontok, melyek Budapestről ugyanannyi idő alatt érhetők el”. És mit tudtunk meg minderről? „Magyarország izokrontérképe sok tanulságos közlekedési jelenséget ábrázol igen áttekinthető módon. Kimutatja, hogy a Duna ma, a hidak kevés száma miatt, rendkívüli közlekedési akadály, továbbá, hogy nemcsak Erdély, hanem a tengerpart felé irányuló közlekedésünk berendezése is rendkívül hiányos. Ha megépülne a Kulpavölgy és a Grobnik-mező közötti 20 km. hosszú Karszt-alagút, e vidék izokronjai rögtön rendkívül megváltoznának, ábrázolva azt a forradalomszerű változást, melyet a tengerpartnak az Alföldhöz közelebb hozatala okozott.” (A Grobniki karsztmező Fiumétől északra található, a Kulpa pedig a horvát-szlovén határfolyó, a jelek szerint a 20. század elején ezt akarták átfúrni alagúttal, hogy a hajón kirakott árukat ne kelljen felvinni a hegyre, majd a túloldalon le. A csel az, hogy nem közúti vagy vasúti alagútat akartak építeni, amire az ember elsőre gondolna, hogy arra rakják át a fiumei kikötőben az árut, ilyen puhányoknak való megoldásnál jóval ambíciózusabb célt tűztek ki: hajóknak való alagutat terveztek! Magyarán az óceánt akarták bevezetni a hegy alá, hogy a hajók egyszerűen behajózzanak az alagútba a fiumei kikötőben, majd a hegy túloldalán kibukkanjak a hajó-alagútból és egy zsiliplépcsővel átlépjenek a Kulpába. A Kulpa a Száva mellékfolyója, az pedig a Dunáé, azaz a tervvel közvetlen összeköttetést teremtettek volna a Duna és az Adriai tenger között. Ha lehet hinni a Wikipedianak, akkor ilyen hajó-alagút még ma, 2022-ben sem létezik egy sem a világon; a norvég Stad Hajó-alagút van a legközelebb a megvalósuláshoz, az építésének a megkezdését 2023-ra tervezik. Még ha el is készül, sehol sem lesz a tervezett Karszt-alagúthoz képest, ugyanis a hossza mindössze 1800 méter…)
Kicsit előbbre ugorva, 1917-ből már térképként rajzolt hazai izokrónt is találunk: a Bátky Zsigmond szerkesztette „Zsebatlasz naptárral és statisztikai adatokkal az 1918. évre” (Magyar Földrajzi Intézet R-T. 1917) közli Kogutowicz Károly „Hazánk izokrón-térképei” című írását, mely vasúti eljutás alapján közöl izokrónokat. Például a tényleges menetrenden alapuló, de a legnagyobb elérhető sebesség alapján számolt térkép így néz ki:
(Fellelnem nem sikerült, de úgy tűnik hogy ezt már az 1913-as zsebatlasz is tartalmazta.)
Az 1930-as évekből aztán már mindenféle specializált térképet is lehet találni, például Debrecenről vagy akár a Balatonról.
Végezetül megjegyzem, hogy léteznek „visszamenőleges” izokrón térképek is: Czére Béla egy 1991-es közleményében 1847-ig visszamenő időpontokra vonatkozóan szerkesztett magyar izokrón-térképeket, az egykorú menetrendek és jármű-paraméterek alapján. Így néz ki például az 1867-es évre rekonstruált budapesti bázisú izokrón térkép:
Most hogy ismerjük az elérési időket, egy kérdés nagyon adja magát: mi ezek közül a legnagyobb, az összes lehetséges település-párt vizsgálva? (A legkisebb nem annyira izgalmas, hiszen vannak nagyon közeli települések, pláne, hogy a budapesti kerületek is külön településnek számítanak a mostani vizsgálatban.) Azaz, melyik két település között vesz a legtöbb időt igénybe az eljutás? És mennyi ez az idő?
knitr::kable(durationsLong[which.max(Duration)])| Var1 | Var2 | Duration |
|---|---|---|
| Kishódos | Felsőszölnök | 7.426778 |
A válasz tehát: Felsőszölnök és Kishódos. Ez a két legtávolabbi magyar település közúton és menetidőről beszélve; majdnem 7 és fél óra az egyikről átjutni a másikra! (Fontos feleleveníteni, hogy ez forgalom nélküli, de szabályokat betartó, gépjárművel történő eljutási idő.) Érdemes ezt az eredményt a biztonság kedvéért más módszerrel is megerősíteni; így néz ki a Google Maps útvonalterve erre a település-párra:
Ezt, amit most kiszámoltunk, a gráfelméletben úgy szokták hívni, hogy mekkora a gráf átmérője. Itt érdemes egy rövid kitérőt tenni. A matematikusok leegyszerűsítve gráfnak neveznek pontokat (ezek szép neve: csúcs), amelyek közül bizonyosak egy vonallal (szép néven: él) vannak összekötve. Az élekhez esetleg számokat is rendelhetünk, melyek az adott összeköttetést jellemzik, szokás ezt néha súlynak is nevezni. Azonnal látszik, hogy ez a matematikai objektum kiválóan megfelel a mostani problémánk leírására, sőt, lényegében közvetlenül alkalmazható: a csúcsok a települések, él pedig azon települések között van, amelyek között el lehet jutni közúton (jelen esetben minden csúcs minden csúccsal össze van kötve, hiszen valahogy minden magyar településről el lehet jutni minden településre), a súly pedig az eljutási idő.
Ahogy említettem, az átmérő egy általános gráfelméleti tulajdonság, mi most a fentiekben ezt számoltuk ki. Természetesen egy sor más (és nem kevésbé izgalmas, sőt) gráfelméleti vizsgálati lehetőség van; érdekes továbbfejlesztési lehetőség kipróbálni, hogy ezek mit fednek fel a magyar közúthálózatról.
Itt csak egyetlen szempontra hívnám fel a figyelmet: az a gráf, amit mi most használunk, nem a fizikai infrastruktúra gráfja (ahol a közvetlenül elérhető, közút által közvetlenül, tehát más települést nem érintő módon összekötött pontok vannak a gráfon is összekötve). Ezt mi kiegészítettük, jobban mondva áttranszformáltuk azáltal, hogy kiszámítottuk a leggyorsabb utakat az összes lehetséges pont között. Amiről itt szó van, az lényegében az, hogy ugyanazt az infrastruktúrát hogyan reprezentáljuk gráffal, a fizikailag összekötött pontok legyenek összekötve, az elérhető pontok legyen összekötve stb. Amint láttuk, ezek között átjárás is lehetséges. Nem arról van szó, hogy valamelyik ezek közül „jó”, más meg „rossz”, egyszerűen más szempontból írják le az infrastruktúrát, így mást árulnak el az országról. (Ennek tömegközlekedési hálózatok kapcsán kiterjedt irodalma van, L terű reprezentáció az, ahol az útvonalon egymás után következő megállók vannak összekötve, P terű reprezentáció pedig az, ahol egy járat által érintett összes megálló össze van kötve egymással.) E kérdéskör alaposabb vizsgálata szintén egy lehetséges továbbfejlesztési irány.
Érdekes kérdés lehet az is, hogy mi az elérési idők eloszlása egy adott településről kiindulva. (Lényegében a térkép színei mögött lévő számokat vesszük ki, és vizsgáljuk meg.) Például Csobád esetében ez így néz ki:
ggplot(durationsLong[Var1=="Csobád"], aes(x = Duration)) +
geom_histogram(boundary = 0, bins = ceiling(log2(nrow(durationsLong[Var2 > Var1])))+1) +
scale_y_continuous(label = scales::comma) +
labs(x = "Elérési idő [h]", y = "Települések száma", caption = captionlab)
Ez egy hisztogram: azt mutatja meg, hogy hány település van, ami a vizsgált kiindulóponttól egy adott elérési időre van, pontosabban szólva az elérési idők egy adott tartományába esik (például hány olyan van, aminek az elérési ideje 10 és 20 perc között van). Ez tehát egy sűrűség-jellegű ábra, hiszen azt mutatja, hogy adott elérési idő környékén milyen sűrűn fordulnak elő települések. A másik lehetséges ábrázolás az eloszlásfüggvény-jellegű ábra, ami nem azt mutatja, hogy hány település van adott elérési idő környékén, hanem azt, hogy mennyi van adott elérési idő alatt. Érdekes módon, bár általában adatvizualizációra jobb az előbbi, mert szabad szemmel sokkal felismerhetőbben mutatja a mintázatokat, itt most az utóbbi is jól megragadható: hány települést tudunk elérni adott időhatáron belül? Avagy, mivel a dolog szimmetrikus, hány településről vagyunk elérhetőek adott időkorláton belül. (Lényegében tehát azt kérdezzük, hogy hány település van a vizsgált település körüli adott sugarú körön belül, azzal, hogy itt nem „körzővel felmért” valódi körről beszélünk, hanem az izokrón vonalról, tehát, ahol nem a távolság állandó adott ponttól, hanem az elérési idő.) Íme a válasz Csobád esetében:
ggplot(durationsLong[Var1=="Csobád"], aes(x = Duration)) +
stat_ecdf(geom = "step", pad = FALSE) +
scale_y_continuous(label = scales::percent,
sec.axis = sec_axis(~ . * length(unique(durationsLong$Var1)),
name = "Elért települések száma [db]",
label = scales::label_number())) +
labs(x = "Elérési idő [h]", y = "Elért települések aránya [%]", caption = captionlab)
Van azonban egy más jellegű problémakör is ezekkel az ábrákkal: nem igazán az a releváns mutató, hogy „darabra” számoljuk meg az adott időn belül elérhető településeket. Sokkal kézenfekvőbb lenne azt kérdezni: hány lakost tudunk elérni adott időn belül? Mivel a települések lélekszáma ismert, így ez is könnyen kiszámolható! Íme a válasz, továbbra is Csobád esetére:
durationsLongPop <- merge(durationsLong[, .(Var1, Helység.megnevezése = Var2, Duration)], HNTdata,
sort = FALSE)[, .(Var1, Helység.megnevezése, Duration, Pop = Lakó.népesség)]
ggplot(durationsLongPop[Var1=="Csobád"], aes(x = Duration, weights = Pop)) +
StatCompLab::stat_ewcdf(geom = "step", pad = FALSE) +
scale_y_continuous(label = scales::percent,
sec.axis = sec_axis(~ . * sum(HNTdata$Lakó.népesség),
name = "Elért lakosok száma [fő]",
label = scales::label_number())) +
labs(x = "Elérési idő [h]", y = "Elért lakosok aránya [%]", caption = captionlab)
Ez az ábra ráadásul nagyon jól összehasonlíthatóvá teszi a településeket. Vessük például össze Csobádot a korábban szintén emlegetett, központi elhelyezkedésű budapesti 1. kerülettel:
ggplot(durationsLongPop[Var1%in%c("Csobád", "Budapest 01. kerület")],
aes(x = Duration, weights = Pop, group = Var1, color = Var1)) +
StatCompLab::stat_ewcdf(geom = "step", pad = FALSE) +
scale_y_continuous(label = scales::percent,
sec.axis = sec_axis(~ . * sum(HNTdata$Lakó.népesség),
name = "Elért lakosok száma [fő]",
label = scales::label_number())) +
labs(x = "Elérési idő [h]", y = "Elért lakosok aránya [%]", color = "", caption = captionlab) +
theme(legend.position = "bottom")
Érdemes megfigyelni egyrészt, hogy az ábra milyen jól mutatja, hogy az ország valóban sokkal „közelebb van” az 1. kerületből, másrészt, hogy bár a színezett térkép nagyon látványos, bizonyos szempontból az ilyen ábra jobb összehasonlíthatóságot teremt meg.
Végezetül, ha nagyon belelendülünk, akkor akár azt is megtehetjük, hogy ezt a görbét az összes településre kirajzoljuk. Ilyenkor érdemes a vonalakat erősen átlátszóvá tenni, különben – a nagyon sok görbe miatt – lényegében egyetlen fekete masszát kapnánk. Az átlátszóság révén viszont érzékelhetővé válik, hogy hol húzódnak sűrűbben a görbék, és hol ritkábban:
ggplot(durationsLongPop,
aes(x = Duration, weights = Pop, group = Var1)) +
StatCompLab::stat_ewcdf(geom = "step", pad = FALSE, alpha = 0.01) +
scale_y_continuous(label = scales::percent,
sec.axis = sec_axis(~ . * sum(HNTdata$Lakó.népesség),
name = "Elért lakosok száma [fő]",
label = scales::label_number())) +
labs(x = "Elérési idő [h]", y = "Elért lakosok aránya [%]", caption = captionlab)
Most, hogy tudjuk, hogy adott településről mennyi idő alatt érhető el az összes többi, elég kézenfekvően adja magát a kérdés, hogy honnan érhető el átlagosan a legkönnyebben az ország? Azaz vesszük egy adott településnél az összes többi település elérési idejének az átlagát, és azt kérdezzük, hogy ez a szám melyik magyar településre lesz a legkisebb. Esetleg megtehetjük, hogy nem az átlagot vesszük, hanem a maximumot, ekkor arra a kérdésre válaszolunk, hogy melyik az a település amire a tőle a legtávolabbi település elérése a legkisebb? (Pesszimisták vagyunk, nem az érdekel minket, hogy átlagosan mennyi a többi település elérése, hanem, hogy legrosszabb esetben mennyi.) Ez utóbbi szokták worst-case számításnak hívni, a matematikusok inkább azt mondják, hogy minimax szabály.
És még egy szempont ehhez: joggal vetődik fel ugyanaz a kérdés, ami az előbb is, tehát, hogy a fenti számítások során nem ugyanúgy számítanak a különböző települések – könnyen lehet, hogy valaki számára nem ugyanolyan fontos, hogy egy adott települéről mennyi idő Budapestre eljutni, és mennyi Csobádra. Ez persze részben szubjektív (éppenséggel lehet valakinek a Budapestre eljutás kevésbé fontos), de egy objektív, és nagyon sok esetben előkerülő, sok mércével korreláló mérőszáma azért van a fontosságnak: a lélekszám. Ha azt mondjuk, hogy a nagyobb lélekszámú településekre fontosabb az eljutási idő, akkor egyszerű a megoldás: sima átlag helyett használjunk – lélekszámmal – súlyozott átlagot.
Ennyi felvezetés után nézzük az eredményeket!
A súlyozatlan átlagos közúti elérési idő az egyes településekről:
AccessMetrics <- durationsLongPop[, .(unweighted = mean(Duration),
weighted = weighted.mean(Duration, Pop),
max = max(Duration)), .(NAME = Var1)]
contourplotter(merge(AccessMetrics[, .(NAME, value = unweighted)], locs, by = "NAME"), geodata,
"Súlyozatlan átlagos\neljutási idő [h]")
A legjobb 10 elérési idejű település e metrika szerint:
knitr::kable(head(AccessMetrics[order(unweighted),
.(`Település` = NAME,
`Súlyozatlan átlagos elérési idő (h)` = unweighted)],
10), digits = 2)| Település | Súlyozatlan átlagos elérési idő (h) |
|---|---|
| Tárnok | 2.14 |
| Törökbálint | 2.15 |
| Budapest 09. kerület | 2.16 |
| Budapest 11. kerület | 2.16 |
| Budapest 07. kerület | 2.17 |
| Budaörs | 2.17 |
| Budapest 06. kerület | 2.17 |
| Budapest 22. kerület | 2.17 |
| Diósd | 2.17 |
| Szigetszentmiklós | 2.18 |
A legrosszabb 10:
knitr::kable(tail(AccessMetrics[order(unweighted),
.(`Település` = NAME,
`Súlyozatlan átlagos elérési idő (h)` = unweighted)],
10), digits = 2)| Település | Súlyozatlan átlagos elérési idő (h) |
|---|---|
| Tiszabecs | 4.77 |
| Zajta | 4.79 |
| Uszka | 4.81 |
| Milota | 4.81 |
| Kispalád | 4.82 |
| Méhtelek | 4.84 |
| Magosliget | 4.87 |
| Nagyhódos | 4.89 |
| Kishódos | 4.90 |
| Garbolc | 4.90 |
A lélekszámmal súlyozott átlagos közúti elérési idő:
contourplotter(merge(AccessMetrics[, .(NAME, value = weighted)], locs, by = "NAME"), geodata,
"Lélekszámmal súlyozott\nátlagos eljutási idő [h]")
A legjobb 10 elérési idejű település e metrika szerint:
knitr::kable(head(AccessMetrics[order(weighted),
.(`Település` = NAME,
`Lélekszámmal súlyozott átlagos elérési idő (h)` = weighted)],
10), digits = 2)| Település | Lélekszámmal súlyozott átlagos elérési idő (h) |
|---|---|
| Budapest 09. kerület | 1.51 |
| Budapest 07. kerület | 1.52 |
| Budapest 06. kerület | 1.52 |
| Budapest 05. kerület | 1.54 |
| Budapest 14. kerület | 1.55 |
| Budapest 23. kerület | 1.55 |
| Budapest 01. kerület | 1.55 |
| Budapest 13. kerület | 1.56 |
| Budapest 11. kerület | 1.56 |
| Szigetszentmiklós | 1.57 |
A legrosszabb 10:
knitr::kable(tail(AccessMetrics[order(weighted),
.(`Település` = NAME,
`Lélekszámmal súlyozott átlagos elérési idő (h)` = weighted)],
10), digits = 2)| Település | Lélekszámmal súlyozott átlagos elérési idő (h) |
|---|---|
| Tiszabecs | 4.24 |
| Zajta | 4.26 |
| Uszka | 4.28 |
| Milota | 4.28 |
| Kispalád | 4.30 |
| Méhtelek | 4.31 |
| Magosliget | 4.34 |
| Nagyhódos | 4.36 |
| Garbolc | 4.37 |
| Kishódos | 4.37 |
A worst case (minimax) elérési idő:
contourplotter(merge(AccessMetrics[, .(NAME, value = max)], locs, by = "NAME"), geodata,
"Legrosszabb\neljutási idő [h]")
A legjobb 10 elérési idejű település e metrika szerint:
knitr::kable(head(AccessMetrics[order(max),
.(`Település` = NAME,
`Legrosszabb elérési idő (h)` = max)],
10), digits = 2)| Település | Legrosszabb elérési idő (h) |
|---|---|
| Budapest 15. kerület | 3.76 |
| Budapest 14. kerület | 3.83 |
| Budapest 04. kerület | 3.84 |
| Budapest 16. kerület | 3.86 |
| Ecser | 3.86 |
| Mogyoród | 3.86 |
| Fót | 3.86 |
| Maglód | 3.86 |
| Budapest 06. kerület | 3.87 |
| Budapest 13. kerület | 3.87 |
A legrosszabb 10:
knitr::kable(tail(AccessMetrics[order(max),
.(`Település` = NAME,
`Legrosszabb elérési idő (h)` = max)],
10), digits = 2)| Település | Legrosszabb elérési idő (h) |
|---|---|
| Kispalád | 7.35 |
| Méhtelek | 7.36 |
| Szakonyfalu | 7.39 |
| Magosliget | 7.40 |
| Szalafő | 7.40 |
| Kétvölgy | 7.41 |
| Nagyhódos | 7.42 |
| Felsőszölnök | 7.42 |
| Garbolc | 7.42 |
| Kishódos | 7.43 |
Látható, hogy az összkép a használt mérőszámtól nagyban függetlenül meglehetősen egységes!
Az előbbi probléma, illetve a megoldása elvezet minket egy általánosabb kérdéshez. Hiszen mit jelent a legkisebb átlagos közelségű pont? Azt, ahonnan minden legjobban elérhető, illetve, ezzel egyenértékűen, azt, ami a legjobban elérhető mindenhonnan. Az egyszerűség kedvéért most válasszuk ki a lélekszámmal súlyozott átlagos elérési időt metrikaként; a továbbiakban végig ezt fogom használni, de minden elmondható lenne a többi mutatóval is.
Miért izgalmas ez? Mondok egy hipotetikus kérdést: az országban egyetlen kórház sincs, most akarjuk az elsőt létesíteni. Hová telepítsük, hogy a betegek átlagosan a leggyorsabban beérjenek a kórházba? Ez a kérdés persze nyilvánvaló marhaság (nem egyetlen kórház lesz az országban, azokat nem úgy telepítik, hogy egyszerűen egy ilyen átlagos elérési időt minimalizáljanak, léteznek orvosi szakterületek is a világon, progresszivitási szintek, számíthatnak a történeti adottságok, nem ugyanolyan könnyű bárhová kórházat telepíteni, van meglevő infrastruktúra, a kórházaknak van kapacitás-korlátjuk stb. stb.), de mégis, ennek a megértése az alaplépés a valóéletbeli problémák megoldásához. Ami igenis előjöhet egy kórház telepítésekor, de előjöhet számos más területen, hulladékgyűjtéstől kezdve a felnyitó bányatérségek optimális telepítéséig.
Na de mi a válasz a kérdésre? Természetesen az előző pontban elért eredmény: súlyozott átlagos elérési időt használva Budapest 9. kerületébe érdemes a kórházat telepíteni.
Eddig tehát nincs nagy újdonság, de tegyük a kérdést érdekesebbé: két kórházat akarunk telepíteni! Hová helyezzük ezeket? Az az eljárás, hogy az elsőt lerögzítjük a 9. kerületben, és a második helyét a maradék 3176 település között keressük, nézve, hogy melyik mellett minimális az eljutási idő (vigyázat, most már az egyes települések elérési ideje úgy értendő, hogy az adott településhez közelebbi kórházba mennyi az elérési idő!) sejthetőleg nem lesz túl szerencsés, legalábbis általában. Gondoljunk például egy képzeletbeli hosszú téglalap alakú országra, nagyjából hasonló lélekszámú településekkel. Ekkor az egy kórházas optimum a téglalap felében van, a két kórházas optimum viszont a harmadolópontokban – azaz itt nem kell kórház oda, ahol az egy kórházas optimum volt. Ilyenkor a fenti logikájú (szokták úgy is hívni: mohó) algoritmusok tévútra visznek.
De akkor mi az általában is jól működő módszer?
Két kórház esetén talán még megpróbálhatjuk az összes lehetséges kombinációt végigpróbálni, bár már ekkor is hatalmas lesz a műveletigény: 3177 településünk van, ebből 5045076 párt lehet képezni. Három vagy több kórház esetén ez praktikusan reménytelen: a lehetséges kórháztelepítések száma exponenciálisan nő.
Élhetünk azzal az egyszerűsítéssel, hogy kórházat csak járási székhelyre telepítünk (Budapest esetén minden kerületet ilyennek tekintünk). Ez valószínűleg nem rossz közelítés, és nagyon fontos hangsúlyozni, hogy ezt most nem úgy értem, hogy az egyéb szempontok miatt, hanem pusztán matematikailag, csak az elérési idő minimalizálást tekintve is: a járási székhelyeknek mindig jobb az elérhetőségük, tehát valószínű, hogy a kórház egy járáson belül akkor is odakerülne, ha egyébként megengednénk bármely más járási településre is az elhelyezést. Így már csak 197 lehetőség van. Párból is csak 19306 van, de az exponenciális növekedést nem kerültük el. Kitérőként megjegyzem, hogy ez egy nagyon jó példa arra, ahogy az emberek általában alábecslik az exponenciális növekedés súlyát: első ránézésre úgy tűnhet, hogy az, hogy kb. 3000 helyett csak 200 települést kell végigpróbálgatni, az egészen drámai egyszerűsítés (nagyságrendi). Igen ám, de ha csak eggyel több kórházt akarunk telepíteni, az máris 200-szorosára növeli a végigpróbálandó kombinációk számát – tehát két nagyságrendet ront! Magyarán: a látszólag drasztikus egyszerűsítés annyit jelentett, hogy jó esetben eggyel több kórház optimális elhelyezésére van ilyen módszerrel reális esélyünk. Mondjuk nem 2, hanem 3 helyét tudjuk meghatározni. De 4 kórházhoz már vagy szuperszámítógép, vagy nagyon sok türelem kell, 5 meg annál több kórház még ezekkel is nehéz lesz… (Az illusztráció kedvéért: egy viszonylag erős gépen, 8 magra párhuzamosítva 1 kórház 1 percen belül megvan. 2 kórházhoz picit több mint 5 perc kell, 3 kórház viszont már több mint 8 óra, 4-et ezek után meg sem mertem próbálni.)
Akkor tehát reménytelen a helyzet? Pláne, ha 10, 20 vagy mondjuk 50 kórházat kell optimálisan telepítenünk? Vagy muszáj egzakt megoldás helyett valamilyen közelítéssel beérnünk?
Nem. Szerencsére nem, meghatározhatjuk az optimális telepítést, mégpedig egzaktan, mindenféle közelítés nélkül, és akár 50 kórházra is – csak épp a fenti „brute force” algoritmusnál okosabban kell eljárnunk. (Jó példa arra, hogy megoldhatatlannak tűnő problémák is megoldhatóvá válhatnak, ha jobb eszközt választunk.) Jelen esetben nem is az a fontos, hogy mi magunk kitaláljuk a szuper-algoritmust, hiszen az ilyenek elmélete egy komplett tudományterület, valószínűleg nem sok jó sül ki abból, ha ezt a nulláról megpróbáljuk magunk kitalálni. A fontos e helyett az, hogy átírjuk a problémát egy olyan keretre, aminél ez a – kész – tudás bevethető, amire vannak nagy tudású, rendkívül optimalizált algoritmusok. Ha ugyanis a problémát el tudjuk ilyen keretben mondani, akkor nincs más dolgunk, mint ezeket bevetni a megoldásra.
Jelen esetben a célravezető az, ha a problémát átfogalmazzunk matematikai optimalizálási (lineáris programozási, még pontosabban egészértékű programozási) problémává, majd bevetjük az arra a területre ismert – nagyon hatékony! – algoritmusok valamelyikét. Ennek a matematikai részleteit egy külön pontban mesélem el. Ami számunkra fontos, hogy ez a módszer drámaian gyorsabb, miközben továbbra is egzakt eredményt szolgáltat! Nézzük a futásidőket ezzel a módszerrel (idézzük fel a brute force megoldás számait!):
m <- nrow(durations)
n <- sum(JarasInd)
LocationResult <- readRDS("LocationResult.rds")
ggplot(data.frame(n = 1:length(LocationResult), time = sapply(LocationResult, function(x) x$time/60)),
aes(x = n, y = time)) + geom_point() + geom_line() +
labs(x = "Telepített kórházak száma", y = "Futási idő [min]")
Látható, hogy a dolog teljesen vállalható, még nagyobb kórházszám esetén is! (Sőt, érdekes módon egy ponton túl ez még jót is tesz. Ennek az okát én sem tudom – sajnos az operációkutatáshoz nagyon kevéssé értek – talán a probléma struktúrája válik olyanná, ami segíti az algoritmust. Külön érdekes, hogy mi történik 26 kórháznál: kétszer is lefuttattam két külön gépen, és így nézett ki az eredmény, szóval ez valószínűleg nem valamiféle hiba, hanem ez is a probléma struktúrájából adódik.)
Ezen felbuzdulva akkor nézzük is meg az eredményeket! Piros ponttal jelölöm a kórházak helyét, a háttérben lévő szín – a szokásos simítással – mutatja az elérhetőséget adott kórház-telepítés mellett (itt ez tehát a legközelebbi kórház alapján megy!). A jobb oldalon lévő sáv mutatja, hogy az ország egészének milyen az elérhetősége adott kórház-konfiguráció mellett, mérve ezt az összes kórház (súlyozott) átlagos elérési idejével. Nagyon szépen látszik, ahogy az újabb és újabb kórházak telepítése színezi át az ország térképét (és az is látszik, hogy a mohó megoldás nem optimális, néha le is kell venni meglevő kórházat amikor több kórházas optimumot keresünk):
kkorhazplot <- function(result, k) {
optlocs <- sapply(LocationResult, function(x)
HNTdata[JarasInd,][x$sol[(m*n+1):(m*n+n)]==1, ]$Helység.megnevezése)
p <- contourplotter(merge(melt(
Reduce(function(...) merge(..., by = "NAME"), lapply(1:k, function(i)
setNames(durationsLong[Var1==HNTdata[JarasInd,][which(
result[[k]]$sol[(m*n+1):(m*n+n)]==1),]$Helység.megnevezése[i],
.(NAME = Var2, value = Duration)], c("NAME", paste0("value", i))))),
id.vars = "NAME")[, .(value = min(value)) , .(NAME)], locs), geodata,
"Súlyozott eljutási idő,\nlegközelebbi kórház [h]", title = paste0(k, " kórházas megoldás"),
filllims = c(0, 4))
for(i in 1:k) p <- p + stat_sf_coordinates(data = geodata[geodata$NAME==optlocs[[k]][i],],
inherit.aes = FALSE, color = "red")
p
}
if(!file.exists("korhazplots.rds")) {
korhazplots <- lapply(1:length(LocationResult), function(k) kkorhazplot(LocationResult, k))
saveRDS(korhazplots, "korhazplots.rds")
} else korhazplots <- readRDS("korhazplots.rds")
korhazbars <- lapply(1:length(LocationResult), function(k)
ggplot(data.frame(y = sapply(LocationResult, function(x) x$obj/sum(HNTdata$Lakó.népesség))[k]),
aes(y = y, x = 1)) + geom_col() + scale_y_continuous(limits = c(0, 2)) +
theme(axis.text.x=element_blank(), axis.ticks.x=element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),) +
labs(x = "", y = "Összesített átlagos eljutási idő [h]"))
for(i in 1:length(LocationResult)) egg::ggarrange(korhazplots[[i]], korhazbars[[i]], ncol = 2,
widths = c(7, 1))
Érdemes külön is ábrázolni az országos átlagos elérési idő javulását a kórházak számának függvényében:
ggplot(data.frame(n = 1:length(LocationResult),
weighted = sapply(LocationResult, function(x) x$obj/sum(HNTdata$Lakó.népesség))),
aes(x = n, y = weighted)) + geom_point() + geom_line() +
labs(x = "Telepített kórházak száma", y = "Súlyozott átlagos elérési idő [h]")
Az eredményeket a fenti vizualizációkon túl konkrétan is letölthetővé tettem (a pontos futásidőkkel, elérési időkkel, és a településlistákkal).
Megjegyzem, hogy ez egy – nagyon primitív és végletekig leegyszerűsített, de mégis csak – példa arra, amit szükséglet-alapú tervezésnek szoktak nevezni: a kórházakat az alapján tervezzük (jelen esetben térben), hogy milyen szükséglet van rájuk (térben). A fenti számítás azért nagyon leegyszerűsített, mert egyéb információ híján a „szükséglet”-et egész egyszerűen a lélekszámmal méri (hiszen ez jelent meg a súlyozásban az átlagolásnál). Ez nem irreális, elvégre ahol többen laknak, ott valószínűleg a kórház igénybevétele is nagyobb, de azért jócskán finomítható lenne: egy diabetológiai osztálynál használhatnánk a cukorbeteg lakosok számát, egy kardiológiai osztálynál a szívbetegek számát, egy gyerekosztálynál csak a gyerekek számát stb. Én most azon egyszerű költői oknál fogva maradtam a lélekszámnál, mert egyedül erre volt megbízható, nyilvános, település-szintű információ. (Esetleg a korcsoportos szűkítéssel lehet játszani, ha erre valaki érez késztetést, ilyen adatok ugyanis letölthetőek település-szinten is a KSH-tól.) Mindenesetre a módszertant ez is demonstrálja – ezt lehetne használni a komplexebb, szükségletek alapján történő racionális tervezés esetén is.
(Itt muszáj közbeszúrnom egy megjegyzést. A fentiekben magától értetődő kérdésnek vettem, hogy tudjuk hány cukorbeteg van adott településen, de sajnos ez közel nem ilyen egyszerű. Az ilyen jellegű információk egyik alapvető forrását az ún. betegségregiszterek jelentik: ezek a teljesség igényével gyűjtik egy országban előforduló eseteit adott betegségnek. Magyarországon a két legismertebb ilyen regiszter a Nemzeti Szívinfarktusregiszter és a Nemzeti Rákregiszter. Ezek előnye a teljeskörűségen túl az, hogy részletgazdag klinikai információk érhetőek el a betegekről, hiszen az orvos mindegyik esethez kitölt egy kérdőívet, valamint, hogy – pont ebből fakadóan, hiszen orvos tölti ki, az általa épp akkor ellátott betegről – jó az adatminőség is. (Bár közel nem tökéletes, még a nagy hagyományú Rákregiszterben is érdemleges zűrök vannak.) A dolog ára, hogy mindez rettenetesen munkaigényes: nagyon terheli az orvosokat az említett adatszolgáltatás, amit egyesével, kézzel kell végrehajtani minden betegre, így ilyen regisztert működtetni csak a legfontosabb néhány betegségre lehet. Éppen ezért vált egyre népszerűbbé az utóbbi időben az ún. adminisztratív/finanszírozási adatok használata. Ehhez tudni kell, hogy minden betegellátási epizódot elektronikus úton jelenteni kell a finanszírozó NEAK-nak, márpedig ebben a jelentésben – a beteg adatain túl – szerepel a diagnózis, vagy az elvégzett beavatkozás kódja is. Akkor, kérdezhetné valaki, miért nem gyűjtjük ki ezekből – például – a cukorbetegség kódjait, és meg is vagyunk. A dolog azért nagyon csábító, mert az adatbázis eleve rendelkezésre áll, így nulla plusz munkával, minden külön beruházás, vagy adatszolgáltatási plusz-terhelés nélkül tudunk válaszolni, lényegében ügyesen felhasználva a más cél miatt amúgy is begyűjtött adatokat. És valóban, a cukorbetegség kódjainak legyűjtése pár perc, sőt, tulajdonképpen az sem túl nehéz, hogy ezt akár az összes betegségre megtegyük. Mindez nagyon jól hangzik, csak két probléma van. Az egyik, hogy az adatminőség rossz. Egyfelől jóhiszemű okok miatt, a kódokat nagyon könnyű elgépelni, másrészt rosszhiszemű okok miatt, ugyanis a jelentéseket a kórházak orrba-szájba hamisítják (bocsánat: optimalizálják) épp azért, mert a finanszírozásuk múlik ezen. Ez ellen persze lehet tenni: a hamisítás módszereit a túloldalon is jól ismerik, be lehet lőni, hogy mi „kemény” adat, ami azért nem valószínű, hogy nagyon odébbmászna, a véletlen tévedések ellen is lehet védekezni, például megkövetelni, hogy egy kód többször is előforduljon, esetleg adott időszakon belül, külön lehet vizsgálni a fekvő- és járóbeteg ellátást, sőt, akár a gyógyszerkasszát is hozzá lehet venni, ha mondjuk a betegséghez biztosan kell szedni valamilyen gyógyszert, de még olyanra is van példa, ahol úgy próbálták körbebástyázni a dolgot, hogy azt nézték, hogy adott diagnosztikai eljárást elvégeztek-e a betegnél, mert ha tényleg a vizsgált baja volt, akkor az kötelező. Ennek milliónyi részlete van, most az egésszel csak annyira akartam kilyukadni, hogy nagyon-nagyon nem annyi a dolog, hogy beírjuk, hogy „kód = X”, majd ütünk egy entert és a kezünkbe hullik a válasz, a jó minőségű lekérdezéshez a területet orvosilag ismerő szakértőkkel kell konzultálni, finanszírozás rejtelmeit ismerő szakértőkkel kell konzultálni, és könnyen lehet, hogy a végén egy kilométer hosszú lekérdezés fog kijönni. Ami ráadásul csak az adott konkrét betegségre lesz jó, értelemszerűen, egy másik betegségnél lehet az egészet elölről kezdeni. A másik probléma, hogy az adminisztratív/finanszírozási adatok nem tartalmaznak klinikai részleteket. Nem lehet megtudni, hogy az alany dohányzik-e, sportol-e, továbbá bár azt látjuk, hogy készült mondjuk egy röntgenje, de hogy azon mi látszódott, arról fogalmunk sem lesz. Cserében azonban ott van az a hatalmas előny, hogy az ország praktikusan minden lakosáról, hosszú időbeli átfogással, többlet-adatszolgáltatás nélkül lehet információt szerezni. Összességében véve az adminisztratív/finanszírozási adatok felhasználása egy csábító lehetőség, de sok csapdával, amiket tudatosan kell végiggondolni. És még egyetlen dolog a végére: akármilyen megoldást is választunk, alapvető kérdésnek tartom a transzparenciát. Mind a használt módszertanra vonatkozóan, melyet széleskörű, nyilvánosan lefolytatott szakmai egyeztetés alapján kell meghatározni, mind a végeredmény tekintetében: a kapott eredményeket, természetesen az adatvédelmi szabályoknak megfelelő módon – például kizárólag járási incidenciákat közölve – bárki számára elérhetően, kutathatóan, részletes magyarázattal, meta-adatokkal ellátva közölni kell. Kérdés, hogy ennek megvalósulása a mai magyar helyzetben mennyire csak illúzió a részemről…)
És a végére még egy gondolat. Valakiben felmerülhet az ötlet, hogy akkor miért nem határozzuk meg a szükséges kórházak optimális darabszámát magát is egy ilyen matematikai optimalizációs modellel? A dolog minden további nélkül megvalósítható (sőt, tulajdonképpen az optimalizációs feladat egyszerűsítését jelenti, ugyanis a fentiekben plusz korlátként kellett belerakni, hogy adott számú kórház legyen).
Az előbbi grafikon már sejteti erre a választ: akkor azt kapnánk, hogy minden településre (járási székhelyre) kell egy kórházat telepíteni! De ha meggondoljuk, ez az eredmény tökéletesen logikus: ha kizárólag az a célfüggvény, hogy minél jobb legyen az átlagos elérési idő, akkor tényleg az a jó, ha mindenhová építünk egy kórházat! Ezen a ponton muszáj egy pillanatra megállnunk, és ezt a kérdést kicsit közelebbről megvizsgálnunk. A probléma ugyanis első ránézésre csak anyagi: persze, hogy ez az eredmény jött ki, hiszen modellünkben ingyen van a kórház telepítése – akkor még szép, hogy azt mondja, hogy mindenhová építsünk egyet. Valójában azonban bonyolultabb a helyzet: ha egy jótündér kifizetné az országnak tetszőleges számú kórház megépítését, akkor sem biztos, hogy az lenne a jó ötlet, ha mind a 3177 településre építenénk egy kórházat. Egy másik írásomban részletesen megtárgyaltam ezt a kérdést, ezt bátran merem ajánlani az érdeklődőeknek, röviden összefoglalva három baj van. Egyrészt a rengeteg kórház fenntartása nyilvánvalóan nem gazdaságos. Ez a legkisebb probléma, ezt figyelembe vehetnénk a költségben, vagy az egyszerűség kedvéért mondjuk most azt, hogy a jótündér ezt is fizeti. De még ekkor sem jó ötlet rengeteg kórházat telepíteni! A második gond, hogy az orvosok és ápolók szintén szét fognak töredezni, a méretgazdaságosság nem csak pénzben üt vissza: nem lesz elég ember ennyi kórházat működtetni. Sőt, valójában nem csak az emberekről van szó, hanem a berendezésekről, műszerekről: ha rengeteg kórház van, akkor nem fogjuk tudni mindegyikbe megvenni a legkorszerűbb eszközöket (vagy, ha mégis megtesszük, akkor meg ott fognak állni kihasználatlanul). De ha ezt is félrerakjuk, akkor jön a harmadik és legfontosabb probléma: az, hogy ez esetben nagyon kevés esetet fog ellátni egy kórház. De ez miért baj? – kérdezhetné valaki. Nem direkt jó? Hiszen akkor nincs kapkodás, nyugodtan lehet dolgozni az orvosoknak, nem kell sokat várnia a betegeknek…! Igen, ezek jogos szempontok, és lesznek esetek, amikor ezek dominálnak. Csak sajnos van egy ellentétes szempont is, és az esetek egy jelentős részében ez lesz a (sokkal) fontosabb: hogy a személyzet így nem tud gyakorlatot szerezni! Különböző beavatkozások egész során dokumentálták empirikus kutatásokkal, hogy a nagyobb kórházi ellátási volumen javítja a kimenetet. Illusztratív példaként egy részlet egy ma már klasszikus kutatás eredményeiből:
Az ábrán az látszik, hogy négy különböző műtéttípus esetén (ez a négy csoportja az oszlopoknak) hogyan alakul a – 30 napon belüli vagy kórházi – halálozás a szerint, hogy az adott kórház mennyi műtét végez a kérdéses típusból (minél világosabb az oszlop, annál többet). Az eredmény nem igényel sok kommentárt: a növekvő gyakorlattal csökken, esetenként egészen drámai mértékben a halálozás. (Természetesen felvethető, hogy a különböző kórházak más összetételű, például más súlyosságú vagy életkorú betegeket láttak el, de az eredmények erre már korrigálva vannak.)
Szó szerint több ezer hasonló vizsgálat készült az idők folyamán, különböző országokban, különböző időintervallumokban, a legkülönfélébb beavatkozásokra, lényében a fentivel egyező eredményekkel. Ez ma már praktikusan egy önálló kutatási terület („volume-outcome”, azaz volumen-kimenet kapcsolat). Van ilyen eredmény szívmegállástól kezdve térdműtéten át a szülészetig. (Tulajdonképpen minden beavatkozásra, műtétre megvizsgálható ez a kérdés, olyannyira, hogy nemzetközileg több helyen már az is megtörtént, hogy a betegbiztonság érdekében minimális szükséges esetszámot határoztak meg kórházra, vagy akár orvosra is, beavatkozás-típusonként.)
Egy szó mint száz, egy sor betegség, szakterület esetén a túl sok kórház káros a betegeknek, mert rontja az ellátás minőségét, így még akkor is kerülendő, ha minden létező anyagi szempontot félreteszünk.
Valójában az a gond, hogy túlságosan teátrális példa volt a részemről a kórház telepítése. Bár tényleg sokaknak ez jut az eszébe a gyors eljutásról, de az igazság az, hogy a kórházak egy halom dolgot csinálnak, amihez egyáltalán nem fontos gyorsan odajutni. Igen, természetesen léteznek sürgősségi területek, de ez csak egy aspektus. (Talán jobb is lett volna mentőállomást, vagy épp tűzoltóságot hozni példának. Ott egyébként gyakran a minimax-elv a fontosabb: azt kell garantálni, hogy még a legtávolabbi településre is eljusson adott időkorláton belül a mentő, nem csak azt, hogy átlagban gyorsan kiérjen.) Ezen kívül azonban legtöbbször sokkal fontosabb egy kórház esetében az, hogy jól felszerelt, legkorszerűbb műszerekkel rendelkező, a szűkebb szakterületet is jól ismerő, gyakorlott orvosokkal és ápolókkal rendelkező centrum lássa el a betegeket, mint az, hogy könnyen be lehessen jutni. Persze, értem, hogy nem jár a busz, rossz az út stb., de talán akkor az utat meg a buszmenetrendet kell megjavítani, és nem kórházat építeni minden második sarokra… Vegyük figyelembe, hogy a dolog alternatívája az, hogy könnyű bejutni, csak ott valószínűbb, hogy rosszabb ellátást kapunk (drámaibban: valószínűbb, hogy meghalunk). Mi a jó választás? A mai magyar helyzet igazából a szovjet mintájú szemaskó-i szervezés öröksége, kombinálódva azzal a helyzettel, hogy politikai okokból – kevés, a lakosság érdekeit a politikai kár vállalása árán is érvényesítő kivételtől eltekintve – nem lehet kimondani, hogy ágyszámot csökkentünk, pláne kórházat zárunk be, miközben mindenki tudja, hogy erre lenne szükség. De ez már végképp az említett másik írásom témája.
Ha már itt tartunk, valakiben felmerülhet az a gondolat is, hogy akkor miért nem építjük be mindezt a célfüggvénybe? Végső soron ez nem is lenne nehéz, kiszámítjuk, hogy hány beteg jut egy kórházra, majd korábban említett (volumen–eredmény) irodalmi adatok alapján ráeresztünk egy függvényt, ami bünteti azt, ha túl kevés a beteg, megjelenítve, hogy a gyakorlat hiánya miatt ott rosszabb lesz a teljesítmény. (Hangsúlyozom, hogy itt végig csak az orvosi szempontokat vizsgálom, az anyagi aspektus nélkül.) Ez, ha elvileg nem is lehetetlen, gyakorlatilag nagyon zűrősnek tűnik nekem. Hogyan konstruáljuk meg pontosan ezt a függvényt? Milyen az alakja? Mik a paraméterei? A probléma az, hogy valójában rengeteg választási lehetőség van, amiket kérdés, hogy mennyire lehet pontosan megmondani, de közben nagyon fog függeni tőlük az eredmény. Az is felmerülhet, hogy ha ezt beépítjük, akkor már ezer más szempont beépítése is ugyanennyire – vagy még sokkal inkább – indokolt lenne, akkor meg miért pont ezt ragadjuk meg? Összességében azt gondolom, hogy ha valaki nagyon bátor/elszánt, akkor gondolkozhat ilyen modellekben, de ezt a területet eléggé aknamezőnek érzem.
Utazóügynök problémának (travelling salesman problem, TSP) szokták a matematikában azt a kérdést nevezni, hogy egy adott gráfon melyik az az összes pontot érintő körút, melynek a súlya, értve ez alatt az útvonalat alkotó élek súlyainak az összegét, a lehető legkisebb. Jelen esetben ennek egy rendkívül szemléletes tartalma van: milyen útvonalon tudjuk körbejárni a lehető legkevesebb idő alatt az összes magyar települést? (Körutat keresünk: a végén ugyanoda kell érkeznünk, ahonnan indultunk.) Mi az „itiner”, ami a legjobb? Hány órára van ehhez szükség? (Érdemes a megoldás elolvasása előtt tippelni!)
A TSP megoldása általában véve bonyolult számítási szempontból (vájtfülűek kedvéért: NP-teljes). Ez azt jelenti, hogy bár elvi szinten semmilyen problémát nem jelent a megoldása (ha semmi más nem jut az eszünkbe: egyszerűen nézzük végig az összes lehetséges útvonalat), de gyakorlatban a kiszámítása nagyon lassú lehet. Az például rögtön látszik, hogy az előbbi zárójeles megjegyzés rendkívül rossz algoritmus: az összes lehetséges útvonal darabszáma eszement módon nő a gráf méretével, ezért már közepes méretű gráfoknál is teljesen reménytelen így megkeresni az optimumot, bármilyen szuperszámítógépünk is legyen. Ennél eljárhatunk sokkal okosabb módokon is, de az NP-teljesség pont azt jelenti, hogy még a legokosabb megoldás sem tud általában véve gyors lenni (ehhez persze precízen definiálni kellene, hogy mit értünk „gyors” alatt, de ez most már nem fontos nekünk ennyire pontosan). Ezért nagyon gyakran ún. heurisztikus algoritmusokat használnak, amelyek nem feltétlenül az optimális útvonalát adják meg, de egy olyat, ami jellemzően azért súlyban (összidőben) közel van hozzá, cserében viszont sokkal gyorsabban lefut. Szerencsére a mi esetünk, 3177 pont a mai számítástechnikai lehetőségek mellett nem olyan nagy, hogy ne lehetne kivárható idő alatt megoldani egzaktan.
Erre a célra a Concorde TSP
Solver nevű programot használtam,
ami bár kissé koros, ám a mai napig a legjobb algoritmusnak tarják nagy
TSP-problémák egzakt kezelésére. Ennek van egy saját adatformátuma, de
szerencsére a TSPLIB nevű R könyvtárral könnyen kimenthetjük az
adatainkat ebben a formában. Két dologra kell csak figyelni. Az egyik,
hogy érdemes a szimmetrizált távolságmátrixot használni: mivel, mint
láttuk is, a mátrix közel szimmetrikus, így ez sokat nem oszt vagy
szoroz, cserében viszont sokkal-sokkal egyszerűbb lesz a probléma
számítástechnikai szempontból. A másik egy informatikai korlát: a
Concorde kizárólag egész számokat hajlandó megenni távolságként. Ha a
súlyok (idők) nem egész számok, mint ahogy most természetesen nem azok,
akkor úgy oldhatjuk meg a problémát – ezt a TSPLIB beépítetten el is
végzi – hogy a tizedesvesszőt odébbrakjuk (minden számot felszorzunk
TSP::write_TSPLIB(TSP::as.TSP(durationsSymm), "SymmDurationMatTSP.tsp", precision = 4)Ha ez megvan, akkor le is futtathatjuk a Concorde-ot ezen az adatbázison:
./concorde SymmDurationMatTSP.tsp > SymmDurationMatTSP.log
Az eredményül kapott log fájlt elérhetővé tettem; ebből mellesleg az is látszik, hogy a teljes futási idő nálam 69221 másodperc, azaz kicsit több mint 19 óra volt.
De mi az optimális útvonal? Ez a .sol fájlból olvasható ki, egyedül
arra kell figyelni, hogy a Concorde 0-tól kezdi a pontok számozását
(magyarán 1-et hozzá kell adni minden számhoz, hogy megkapjuk a
bemenetben, jelen esetben a locs-ban szereplő elemekre a hivatkozást).
Nem csigázom tovább a kedélyeket, íme a válasz:
TSPsol <- read.table("SymmDurationMatTSP.sol", header = FALSE, skip = 1, fill = TRUE)
TSPsol <- c(t(TSPsol))
TSPsol <- TSPsol[!is.na(TSPsol)]
TSPsol <- TSPsol + 1
TSPsolItinary <- data.table(Var1 = locs$NAME[TSPsol],
Var2 = c(locs$NAME[TSPsol][-1], locs$NAME[TSPsol][1]))
TSPsolItinary <- merge(TSPsolItinary, durationsLong, by = c("Var1", "Var2"), sort = FALSE)
fwrite(TSPsolItinary, "TSPsolItinary.csv", dec = ",", sep = ";", bom = TRUE)
ggplot(geodata) + geom_sf(color = NA) +
geom_point(data = locs[TSPsol,], aes(x = X, y = Y), size = 0.3) +
geom_path(data = rbind(locs[TSPsol,], locs[TSPsol[1],]), aes(x = X, y = Y)) +
labs(x = "", y = "", caption = captionlab) +
metR::scale_x_longitude(ticks = 1, expand = waiver()) +
metR::scale_y_latitude(ticks = 0.5, expand = waiver())
A teljes itinert innen letölthetővé tettem – az érdeklődők ez alapján már neki is vághatnak hazánk felfedezésének!
Mivel körútról van szó, természetesen bárhol be lehet kapcsolódni, az esetleges, hogy az itiner pont Abán kezdődik. De ha valaki itt indul, akkor látható, hogy először Sárkeresztúrra kell vezetnie, ez kevesebb mint 10 perc, utána mehet Sárszentágotára (5 perc), és így tovább, míg végül, 3177 település érintése után Tácról fog újra visszaérni Abára egy bő 20 perces vezetés után.
Na de mennyi idő mindez? (Dobpergés:) 541.7 óra. Ennyi idő alatt lehet, a legjobb útvonalat választva hazánk valamennyi települését érinteni. (Az érték egy nagyon kicsit eltér a log fájl végén olvashatótól, ez a szimmetrizálás miatt van.) Mivel a gráf itt nem a fizikai infrastruktúrát jelképezi, így elvileg lehet, hogy egy településen többször is végigmegyünk, de nem hinném, hogy ennek nagy jelentősége lenne (könnyen lehet, hogy alaposabb végiggondolással azt is be lehetne bizonyítani, hogy ez elő sem fordulhat).
Összességében tehát azt kaptuk, hogy bő 500 óra alatt hazánk összes települését felkereshetjük autóval!
Felsorolok néhány olyan lehetőséget, amelyek a fentiekben nyitva maradtak, de szerintem a vizsgálatuk érdekes és izgalmas lehet.
- A közúti elérhetőség hosszú távú trendjének vizsgálata: elérhetőek régebbi térképek is az OSM-től? (Ha igen, meddig lehet visszamenni?) Ezzel vizsgálható, hogy javult-e a közúti elérhetőség, mennyit javult, mikor javult, mely területeken javult stb.
- Útvonaltervezés (és az eredmények vizsgálata) nem leggyorsabb, hanem legrövidebb utazásra optimalizálva.
- Más közlekedési módok (gyalog, kerékpár, tömegközlekedés, vasút stb.) vizsgálata. Ezek egy részét az OSRM minden további nélkül támogatja, mert van külön profil, így elvileg pillanatok alatt lecserélhető az egész számítás. Más esetek zűrösebbek, hiszen a vasúthoz vagy városi tömegközlekedéshez a menetrendre is szükség volna.
- Forgalmi adatok figyelembevétele.
- További gráfelméleti lehetőségek vizsgálata (pl. centralitások!).
- Más térben történő gráfreprezentációk vizsgálata.
- Szcenárióelemzések: mely összeköttetés megteremtése (új út építése) jelentené a legnagyobb javulást elérhetőségben?
- Szcenárióelemzések: mely összeköttetés elvesztése (baleset, katonai támadás stb.) jelenté a legnagyobb érvágást?
- Heurisztikus TSP-megoldók tesztelése és egybevetése az egzakt optimummal.
- Átlagos eljutási idő vizsgálatok úgy, hogy csak a saját megyét (vagy saját járást) nézzük.
- Egyéb szükséglet-mutatók, például korcsoportos lélekszámok használata a kórház-telepítési probléma megoldásában.
- Egyéb célfüggvények vizsgálata az optimális kórház-telepítési probléma megoldásában.
- A kórházak vizsgálata kapcsán külön kérdés, hogy itt a valóságban egy, helyenként ráadásul elég komplex beutalási rend érvényesül – ennek figyelembevétele szintén érdekes kérdés lehet. A kórházak területi ellátási kötelezettségének alapvonalairól egy másik írásomban értekeztem dióhéjban.
Magyarország térképét az OpenStreetMap magyar oldaláról tölthetjük le:
if (!file.exists("kozighatarok.zip")) {
unlink("./kozighatarok/", recursive = TRUE)
download.file("http://data2.openstreetmap.hu/hatarok/kozighatarok.zip", "kozighatarok.zip")
unzip("kozighatarok.zip")
}A kapott fájlok ESRI
shapefile
formátumban tartalmazzák a térképeket; ezek beolvasásához az sf
könyvtárat használhatjuk.
A 8-as kód tartalmazza a településeket. Ezt beolvassuk, de mivel Budapest kerületei külön kódon vannak (a 8-as egyben tartalmazza egész Budapestet), így azt ebből eltávolítjuk, ezt majd később kipótoljuk, de kerület szinten:
geodata1 <- st_read("./kozighatarok/admin8.shp", quiet = TRUE)
geodata1 <- geodata1[geodata1$NAME!="Budapest",]Ezt követően beolvassuk (9-es kód) a budapesti kerületeket. Itt a neveket átírjuk rómairól arab számra, majd még egy problémát megoldunk: a 9-es kódú kerületek nem tartalmazzák a Margitszigetet (ezért aztán ha egyszerűen ezt a térképet egyesítenénk az előbbivel, a Margitsziget kimaradna). Szerencsére a 10-es kódú adatbázisban megvan a Margitsziget, így azt egyesítjük a 13. kerülettel, majd az így kapott 13. kerülettel lecseréljük a kerületi szintű térképben szereplő 13. kerületet, magyarán a Margitszigetet „hozzácsatoljuk” a 13. kerülethez. (Szigorúan véve a Margitsziget ma már nem a 13. kerülethez tartozik, de mivel a Helységnévtár – szemben a kerületekkel! – nem tartalmazza önálló soron a Margitszigetet, ezért célszerűbb így kezelni.) A megvalósítás tehát:
geodata2 <- st_read("./kozighatarok/admin9.shp", quiet = TRUE)
geodata2$NAME <- paste0("Budapest ", sprintf("%02d", as.numeric(as.roman(sapply(strsplit(
geodata2$NAME, ".", fixed = TRUE), `[`, 1)))), ". kerület")
geodata3 <- st_read("./kozighatarok/admin10.shp", quiet = TRUE)
geodata2 <- rbind(geodata2[geodata2$NAME!="Budapest 13. kerület",],
st_union(geodata2[geodata2$NAME=="Budapest 13. kerület",],
geodata3[geodata3$NAME=="Margitsziget",])[, c("NAME", "ADMIN_LEVE",
"geometry")])Ez után már egyesíthetjük az így kapott térképet az elsővel, ezzel létrehozva Magyarország település-szintű térképét, melyben Budapest kerületekre van bontva:
geodata <- rbind(geodata1, geodata2)A térkép összesen 3196 településnevet tartalmaz, ezek közül azonban csak 3177 az egyedi. A jelenség magyarázata, hogy bizonyos települések területe több részből áll (és ezek így több soron is szerepelnek a térképfájlban). Például Balatoncsicsó így néz ki:
plot(geodata[geodata$NAME=="Balatoncsicsó", "NAME"], main = "", key.pos = NULL)
Megoldásként egyesítsük ezeket egyetlen objektumba:
geodata <- rbind(geodata[!geodata$NAME%in%geodata$NAME[duplicated(geodata$NAME)],],
do.call(rbind, lapply(geodata$NAME[duplicated(geodata$NAME)], function(nam)
Reduce(st_union, split(geodata[geodata$NAME==nam,],
seq(nrow(geodata[geodata$NAME==nam,]))))[
, c("NAME", "ADMIN_LEVE", "geometry")])))
geodata <- geodata[order(geodata$NAME),]Ellenőrzésként megnézhetjük a térképet amit így kaptunk (itt a települések nem pontok, hanem mindegyikhez egy kis terület van rendelve, hogy együtt kiadják Magyarország egész területét):
plot(geodata["NAME"], main = "")
A térkép koordinátarendszerét átállítjuk WGS84-re (hogy szokásos koordinátákat lássunk; az OSRM is ilyen formában fogja majd várni a koordinákat), és utána el is mentjük a későbbi gyors feldolgozáshoz:
geodata <- st_transform(geodata, crs = "WGS84")
saveRDS(geodata, "geodata.rds")A fenti térkép birtokában összeszedhetjük egy listába Magyarország összes települését (Budapest esetén pedig a kerületeket). Erre azonban van egy másik adatforrásunk is: a Központi Statisztikai Hivatal ún. Helységnévtára! (Régi nevén Helységnévkönyv.) Töltsük ezt le:
if (!file.exists("hnt_letoltes_2022.xlsx"))
download.file("https://www.ksh.hu/docs/helysegnevtar/hnt_letoltes_2022.xlsx",
"hnt_letoltes_2022.xlsx", mode = "wb")Majd olvassuk be:
HNTdata <- readxl::read_excel("hnt_letoltes_2022.xlsx", .name_repair = "universal", skip = 2)
HNTdata <- HNTdata[!HNTdata$Helység.megnevezése%in%c("Összesen", "Budapest"),]
HNTdata$Helység.megnevezése <- ifelse(grepl("Budapest", HNTdata$Helység.megnevezése),
paste0(substring(HNTdata$Helység.megnevezése, 1,
nchar(HNTdata$Helység.megnevezése)-1), "ület"),
HNTdata$Helység.megnevezése)
saveRDS(HNTdata, "HNTdata.rds")Ellenőrizzük le, hogy a két adatforrásból ugyanarra jutottunk-e:
nrow(geodata)## [1] 3177
nrow(HNTdata)## [1] 3177
sum(geodata$NAME!=HNTdata$Helység.megnevezése)## [1] 0
Szerencsére minden stimmel, ugyanaz a listánk, megvan Magyarország mind 3177 települése (Budapest helyett a kerületeket számolva).
A biztonsági ellenőrzésen túl azért volt hasznos a Helységnévtár letöltése, mert így majd később hasznos települési adatokat nyerhetünk ki (például a lélekszámokat, és a járási, megyei hovatartozásokat).
Zárásként nyerjük ki a települések koordinátáit, hiszen az útvonaltervezéshez majd ehhez lesz szükségünk. A térkép, mint fent láttuk is, minden településhez egy kis területet rendelt; legyen a koordináta ennek a középpontja:
locs <- data.frame(NAME = geodata$NAME, st_coordinates(st_centroid(geodata)))
saveRDS(locs, "locs.rds")
write.csv(locs[, c("X", "Y")], "osrmlocs.csv", row.names = FALSE)
knitr::kable(head(locs))| NAME | X | Y |
|---|---|---|
| Aba | 18.51912 | 47.06457 |
| Abádszalók | 20.62460 | 47.45214 |
| Abaliget | 18.09622 | 46.14251 |
| Abasár | 20.01678 | 47.81171 |
| Abaújalpár | 21.25843 | 48.30330 |
| Abaújkér | 21.19252 | 48.31008 |
Az OSRM installálásának a részletei a https://github.com/Project-OSRM/osrm-backend/wiki/Building-OSRM címen elérhetőek. Tételesen, Ubuntu 20.04.4 alatt a következő parancsokat kell lefuttatni:
sudo apt install build-essential git cmake pkg-config \
libbz2-dev libstxxl-dev libstxxl1v5 libxml2-dev \
libzip-dev libboost-all-dev lua5.2 liblua5.2-dev libtbb-dev
Ezt követően szükség van a node.js környezet
installálására is. Arra figyelni kell, hogy olyan verziót telepítsünk,
ami összhangban van azzal, hogy az OSRM mit támogat. Például az OSRM
5.27.1-es
verziója által támogatott legfrissebb node.js a 108-as (ezt a
fájlnévben láthatjuk a -node- után következő számból). Megnézve azt
látjuk, hogy a 108-as szám az
a 18-as verzió, tehát ezt kell
telepítenünk
is:
curl -fsSL https://deb.nodesource.com/setup_18.x | sudo -E bash - &&\
sudo apt-get install -y nodejs
Ez után letöltjük az osrm csomagot a node.js-hez:
npm install osrm
Ezen előkészületek után nekiállhatunk az OSRM telepítésének. Célszerű a
legfrissebb
release-t
használni (a https://github.com/Project-OSRM/osrm-backend.git az éppen
aktuális, fejlesztési verzió, ami lehet jóval kevésbé tesztelt), az
aktuális helyzet esetében ez így néz ki:
wget https://github.com/Project-OSRM/osrm-backend/archive/refs/tags/v5.27.1.tar.gz
tar xzvf v5.27.1.tar.gz
cd osrm-backend-5.27.1/
mkdir -p build
cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
cmake --build .
sudo cmake --build . --target install
cd ..
Az útvonaltervezés gépidő- és memóriaigényes; elképzelhető, hogy swap fájl beállítása is szükséges lesz a későbbiekhez.
A térképet a Geofabrik oldaláról tölthetjük le, egész Magyarország egyben is elérhető:
wget https://download.geofabrik.de/europe/hungary-latest.osm.pbf
Ezt követően néhány előfeldolgozási lépésre van szükség a térképen, hogy alkalmas legyen a későbbi útvonaltervezéshez (most a gépjármű-profilt használjuk, hiszen az egész vizsgálódás erről szól, de pontosan ugyanígy:
osrm-extract hungary-latest.osm.pbf -p ./profiles/car.lua
osrm-partition hungary-latest.osrm
osrm-customize hungary-latest.osrm
Ha mindezek megvannak, akkor elindíthatjuk az OSRM-et. A dolgot úgy kell elképzelni, hogy egy szervert indítunk el, ami folyamatosan várja a kéréseket (a megtervezendő útvonal paramétereit), és válaszként elküldi a megoldást. A szerver elindítása:
osrm-routed hungary-latest.osrm
Ez már eleve egy jobb megoldás, hiszen egy helyi szervertől kérjük le a
válaszokat, és nem interneten keresztül küldjük át. Valójában azonban
nekünk ez sem lesz elég, mert még egy helyi szervernek sem feltétlenül
könnyű ekkora
méretű útvonaltervezést végrehajtania. Éppen ezért egy teljesen más
megoldást választunk: a szervert el sem kell indítani, e helyett más
módon,
egy node.js szkript segítségével futtatjuk le az útvonaltervezést:
node osrmtable.js
Itt az osrmtable.js tartalma:
const OSRM = require('osrm');
const osrm = new OSRM('hungary-latest.osrm');
const fs = require('fs');
const { convertCSVToArray } = require('convert-csv-to-array');
const { convertArrayToCSV } = require('convert-array-to-csv');
var coordinatesin = fs.readFileSync('locs.csv', 'utf8')
coordinatesin = convertCSVToArray(coordinatesin, {
header : false,
type: 'array',
separator: ',',
});
osrm.table({
coordinates: coordinatesin,
annotations: ['duration', 'distance']
}, function(err, response) {
fs.writeFileSync('osrmdurations.csv', convertArrayToCSV(response.durations));
fs.writeFileSync('osrmdistances.csv', convertArrayToCSV(response.distances));
});
Ez a szkript egy locs.csv nevű fájlban várja a koordinátákat, majd az
osrmdurations.csv és az osrmdistances.csv fájlokban adja vissza a
megoldásokat (útvonalak idejei és távolságai minden pár között).
Vigyázat, ez mindenképp időre van optimalizálva, tehát a távolság úgy
értendő, hogy a legrövidebb idejű utak távolsága, nem úgy, hogy a
legkisebb távolság.
A kórházak optimális helyének megválasztása tulajdonképpen nem nagy újdonság az operációkutatásban: telepítési probléma néven tárgyalja az irodalom, egy klasszikus fejezet, amit lényegében közvetlenül alkalmazhatunk a mi mostani problémánk megoldására.
Vezessünk be a következő jelöléseket. Jelölje
Mik a korlátozó feltételeink? Először is, valahogy biztosítani kell,
hogy adott település lakói pontosan egy kórházhoz tartozzanak. Mivel az
A második a kórházak léte. Valahogy nyomon kell követnünk, hogy hol van
kórház: jelölje
Az utolsó dolog, amire figyelnünk kell, hogy az országban telepített
kórházak darabszáma az adott, előírt érték legyen, jelöljük ezt
És ezzel meg is vagyunk! A fentiekben egy lineáris programozási feladatot írtunk elő, pontosabban, mivel minden változónk csak egész értéket vehet fel, egy egészértékű programozási feladatot. (Sőt, binárisat.) Innentől kezdve használhatjuk a megoldására az ilyen feladatok megoldására kidolgozott jól bejáratott, messzemenőkig optimalizált algoritmusokat.
Az első feladat a fenti program átírása olyan formára, amit a
megoldóalgoritmusok megesznek. Ez egyetlen okból nem teljesen direkt: a
döntési változóink (egy része) mátrixban van elrendezve, miközben az
optimalizáló programok csak vektort fogadnak. A megoldás tehát, hogy az
c vektort
csinál egy mátrixból, egyetlen baj, hogy oszloponként, ezért a fentihez
először transzponálnunk kell. Még azt ne felejtsük el, hogy a végére
hozzá kell fűznünk
c(c(t(durs)), rep(0, n))Itt a durs tartalmazza az eljutási időket (durs jelentése lesz más: az eljutási
idők helyett az eljutási idők, szorozva a lélekszámmal (
Mi a helyzet a feltételekkel? Ezeket egy mátrixba kell rendezni;
előrebocsátom, hogy mindenhol ritka mátrix formátumot fogok használni a
megadáshoz (mégpedig a Matrix csomag használatával), ugyanis bizonyos
későbbi megoldó algoritmusok explicite támogatják ezt (és szükség is
lesz rá, hiszen hatalmasak a mátrixok, viszont tényleg ritkák), ha meg
nem, veszteni akkor sem vesztettünk semmit. A mátrix a fent felsorolt
három dolgot kell, hogy tartalmazza: minden település pontosan egy
kórházhoz tartozik, ahová egy települést bekötünk ott legyen kórház, a
kórházak száma
Az első feltétel a
Matrix::sparseMatrix(i = rep(1:m, each = n), j = 1:(m*n), x = 1, dims = c(m, m*n + n))A jobb oldal pedig egyszerűen
A második feltétel az
cbind(Matrix::Diagonal(m*n, 1), Matrix::sparseMatrix(i = 1:(m*n), j = rep(1:n, m), x = -1))A jobb oldal
Végezetül a harmadik feltétel, hogy
Matrix::sparseMatrix(i = rep(1, n), j = ((m*n)+1):(m*n + n), x = 1)A jobb oldal
Ezzel gépi úton feldolgozható formára hoztuk a problémát; kérdés, hogy hogyan oldjuk meg.
Mivel az operációkutatáshoz nem igazán értek, így pusztán egy hevenyészett internetes keresés alapján négy lehetőséget találtam az egészértékű optimalizálás tényleges, számítógépes megoldására. (A kritériumom annyi volt, hogy ingyenes és nyílt forráskódú legyen a program, illetve, hogy R alatt használható legyen.)
Az első feladat annak meghatározása, hogy melyik programcsomagot
alkalmazzuk élesben. A kérdés nem nyilvánvaló, mert félő, hogy a
probléma nagy mérete miatt lesznek olyan algoritmusok, amik
elfogadhatatlanul lassan és/vagy nagy memóriaigénnyel futnak csak le.
Éppen ezért az első feladat ennek a kipróbálása. (Tudom, hogy vannak
mindenféle benchmark-ok, de mivel ehhez nem igazán értek, úgy döntöttem,
hogy empirikusan, a konkrét problémán mérem le a teljesítményeket.)
Azért, hogy feleslegesen ne bonyolítsuk az életünket, az
Nézzük a lehetőségeket!
A GLPK használatát R alatt az Rgplk és a glpkAPI csomag is támogatja; mivel az utóbbi kicsit nyakatekertnek tűnik, az előbbit használtam.
Íme a kód:
solGLPK <- function(durs, k) {
n <- nrow(durs)
result <- Rglpk::Rglpk_solve_LP(
obj = c(c(t(durs)), rep(0, n)),
mat = rbind(Matrix::sparseMatrix(i = rep(1:n, each = n), j = 1:(n*n), x = 1, dims = c(n, n*n + n)),
cbind(Matrix::Diagonal(n*n, 1),
Matrix::sparseMatrix(i = 1:(n*n), j = rep(1:n, n), x = -1)),
Matrix::sparseMatrix(i = rep(1, n), j = ((n*n)+1):(n*n + n), x = 1)),
dir = c(rep("==", n), rep("<=", n*n), "=="),
rhs = c(rep(1, n), rep(0, n*n), k),
bounds = list(lower = list(ind = 1:(n*n + n), val = rep(0, n*n + n)),
upper = list(ind = 1:(n*n + n), rep(1, n*n + n))),
types = rep("B", n*n + n),
max = FALSE)
list(obj = result$optimum, sol = result$solution)
}Látszik, hogy mindez lényegében a fenti gondolatok direkt implementálása (egyedül a változók határainak megadása lehet meglepő, mert azt is ritka mátrix stílusban kell közölni).
Az lp_solve használatát R alatt az
lpSolveApi
és az
lpSolve
csomag teszi lehetővé, ahogy néztem, ez utóbbi nem támogatja a
A kód:
solLPSOLVE <- function(durs, k) {
n <- nrow(durs)
mat1Sparse <- Matrix::sparseMatrix(i = rep(1:n, each = n), j = 1:(n*n), x = 1, dims = c(n, n*n + n))
mat2Sparse <- cbind(Matrix::Diagonal(n*n, 1),
Matrix::sparseMatrix(i = 1:(n*n), j = rep(1:n, n), x = -1))
mat3Sparse <- Matrix::sparseMatrix(i = rep(1, n), j = ((n*n)+1):(n*n + n), x = 1)
Amat <- rbind(mat1Sparse, mat2Sparse, mat3Sparse)
lpmodel <- lpSolveAPI::make.lp(n*n + n + 1, n*n + n)
for(i in 1:ncol(Amat)) lpSolveAPI::set.column(lpmodel, i, Amat[, i])
lpSolveAPI::set.objfn(lpmodel, c(c(t(durs)), rep(0, n)))
lpSolveAPI::set.bounds(lpmodel, rep(0, n*n + n), rep(1, n*n + n))
lpSolveAPI::set.constr.type(lpmodel, c(rep("=", n), rep("<=", n*n), "="))
lpSolveAPI::set.constr.value(lpmodel, c(rep(1, n), rep(0, n*n), k))
lpSolveAPI::set.type(lpmodel, 1:(n*n + n), "binary")
lpSolveAPI::solve.lpExtPtr(lpmodel)
list(obj = lpSolveAPI::get.objective(lpmodel),
sol = lpSolveAPI::get.variables(lpmodel))
}Itt már bonyolultabb kicsit a helyzet, ráadásul – számomra elég érthetetlen módon – nincs olyan utasítás, amivel a megkötések mátrixát egyben át lehetne adni, így kénytelenek vagyunk oszloponként betölteni azt.
A Cbc használatát az R alatt az rcbc csomag támogatja.
A kód:
solCBC <- function(durs, k) {
n <- nrow(durs)
result <- rcbc::cbc_solve(
obj = c(c(t(durs)), rep(0, n)),
mat = rbind(Matrix::sparseMatrix(i = rep(1:n, each = n), j = 1:(n*n), x = 1, dims = c(n, n*n + n)),
cbind(Matrix::Diagonal(n*n, 1),
Matrix::sparseMatrix(i = 1:(n*n), j = rep(1:n, n), x = -1)),
Matrix::sparseMatrix(i = rep(1, n), j = ((n*n)+1):(n*n + n), x = 1)),
row_lb = c(rep(1, n), rep(-Inf, n*n), k),
row_ub = c(rep(1, n), rep(0, n*n), k),
col_ub = rep(1, n*n + n),
col_lb = rep(0, n*n + n),
is_integer = rep(TRUE, n*n + n),
max = FALSE,
cbc_args = list(log = 0)
)
list(obj = result$objective_value, sol = result$column_solution)
}Ez ismét meglehetősen közvetlen átültetés, egyedül arra kell figyelni,
hogy a jobb oldalaknál nem lehet relációt megadni, csak alsó és felső
értéket, így az egyenlőséget két azonos számmal, a kisebb-egyenlőt pedig
az alsó határt
A HiGHS az R alatt a highs csomag segítségével hívható meg.
A kód:
solHIGHS <- function(durs, k) {
n <- nrow(durs)
result <- highs::highs_solve(
L = c(c(t(durs)), rep(0, n)),
lower = rep(0, n*n + n),
upper = rep(1, n*n + n),
A = rbind(Matrix::sparseMatrix(i = rep(1:n, each = n), j = 1:(n*n), x = 1, dims = c(n, n*n + n)),
cbind(Matrix::Diagonal(n*n, 1),
Matrix::sparseMatrix(i = 1:(n*n), j = rep(1:n, n), x = -1)),
Matrix::sparseMatrix(i = rep(1, n), j = ((n*n)+1):(n*n + n), x = 1)),
lhs = c(rep(1, n), rep(-Inf, n*n), k),
rhs = c(rep(1, n), rep(0, n*n), k),
types = rep("I", n*n + n), maximum = FALSE
)
list(obj = result$objective_value, sol = result$primal_solution)
}Ez mindenben az előbbihez hasonló implementáció.
A vizsgálat két szempontra kell hogy kiterjedjen: futásidő és memóriaigény. Érdemes először nagyon kis problémákkal kezdeni, és ránézni az eredményre, és csak ezután menni a nagyobb problémák felé, hátha valamelyik algoritmus már ekkor láthatóan nagy rossz (és ezért a végtelenségig várhatnánk, ha rögtön nagyobb problémákat is próbálnánk). A probléma méretét úgy állítottam be, hogy a teljes magyar mátrixból a meghatározott számú települést véletlenszerűen kiválasztottam, és így futtattam az optimalizálást.
A megvalósításhoz a bench csomagot használtam, ami három, a mostani szempontból fontos előnyt is kínál:
- Egyszerre mér idő- és memóriaigényt.
- Könnyen automatizálható: van függvénye, ami megadott paraméterek összes lehetséges kombinációjára automatikusan lefuttatja a méréseket, és összefogja az eredményeket.
- Több függvény összehaslítása során automatikusan ellenőrzi, hogy mindegyik ugyanazt az eredményt adja-e vissza, és hibát ad ha nem; így egyúttal tartalmilag is tudjuk ellenőrizni, hogy nem rontottuk-e el az implementációkat (vagy legalábbis az összeset azonosan rontottuk el…).
A bench használata során nem kell megadni a szükséges futtatások
számát, adaptív módon határozza meg az időigény alapján.
Amint mondtam, a településeket véletlenszerűen választottam ki az összes
település közül. Igen ám, de egy teljesítménymérés, bár többször
futtatja az optimalizálót, de mindig ugyanazokkal a településekkel.
Mivel azonban az eredmény függhet a települések megválasztásától is,
ezért érdemes az egészet még egy külső ciklusba rakni, és többször
újrafuttatni, de immár különböző településekkel is, hogy ennek a hatását
is lássuk; a lenti kódokban ezt szolgálja a rep változó.
Elsőként mind a négy optimalizáló – egészértékű programozást megoldó – eljárást használjuk, de csak szolid, legfeljebb 90 települést tartalmazó esetekben. Nézzük először az időigényt:
if(!file.exists("benchres1.rds")) {
benchres1 <- bench::press(
n = c(seq(10, 90, 20)),
k = c(1, 5, 10),
rep = 1:3,
{
locs <- sample(1:nrow(durations), n, replace = FALSE)
durs <- durations[locs, locs]
bench::mark(GLPK = solGLPK(durs, k),
LPSOLVE = solLPSOLVE(durs, k),
CBC = solCBC(durs, k),
HIGHS = solHIGHS(durs, k))
}
)
saveRDS(benchres1, "benchres1.rds")
} else benchres1 <- readRDS("benchres1.rds")
benchres1DF <- with(benchres1, data.frame(n, k, rep, median = as.numeric(median),
memory = as.numeric(mem_alloc)/1024/1024/1024,
solver = bench:::as.character.bench_expr(expression)))
ggplot(benchres1DF, aes(x = n, y = median, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k) + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Időigény [s]", color = "Megoldó")
Látszik, hogy az lp_solve eljárás érezhető lassabb mint a másik három.
Vizsgáljuk meg a memóriaigényt is:
ggplot(benchres1DF, aes(x = n, y = memory, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k) + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Memóriaigény [GB]", color = "Megoldó")
Az lp_solve memóriaigénye valami egészen elmebeteg módon elszáll,
tehát biztosan nem érdemes vele számolni a továbbiakban.
Nézzük tehát a maradék három eljárást, immár egy picit nagyobb problémákra is, elmenve 200 településig. Szokásosan kezdjük a futás időigényével:
if(!file.exists("benchres2.rds")) {
benchres2 <- bench::press(
n = c(seq(30, 120, 30), 150, 200),
k = c(1, 5, 10, 30),
rep = 1:3,
{
locs <- sample(1:nrow(durations), n, replace = FALSE)
durs <- durations[locs, locs]
bench::mark(GLPK = solGLPK(durs, k),
CBC = solCBC(durs, k),
HIGHS = solHIGHS(durs, k))
}
)
saveRDS(benchres2, "benchres2.rds")
} else benchres2 <- readRDS("benchres2.rds")
benchres2DF <- with(benchres2, data.frame(n, k, rep, median = as.numeric(median),
memory = as.numeric(mem_alloc)/1024/1024,
solver = bench:::as.character.bench_expr(expression)))
ggplot(benchres2DF, aes(x = n, y = median, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k) + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Időigény [s]", color = "Megoldó")
Látszik, hogy a GLPK időigénye némileg nagyobb a többinél. Az is érzékelhető, hogy érdekes módon a nagy számú kórház optimalizálása gyorsabban fut.
Nézzük most a memóriaigényt:
ggplot(benchres2DF, aes(x = n, y = memory, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k) + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Memóriaigény [MB]", color = "Megoldó")
Itt viszont a GLPK már nagyon egyértelműen kilóg, így ha az előbbi után volt is kétségünk, ez bizonyossá teszi, hogy tovább már csak a maradék két optimalizálóval érdemes menni. Jól látszik egyébként, hogy a különböző futtatások között memóriaigényben nem volt különbség (a három görbe egybeesik), de ez logikus is, ha meggondoljuk, hiszen a memóriaigény szempontjából mindegy, hogy konkrétan mely településeket választjuk ki.
A harmadik és egyben utolsó körbe, a „döntőbe” tehát a Cbc és a HiGHS jutott be. Itt már végképp bátrak vagyunk, és a pontos feltérképezés érdekében egészen az 500-as településszámig elmegyünk. Nézzük a futásidőt:
if(!file.exists("benchres3.rds")) {
benchres3 <- bench::press(
n = seq(100, 500, 100),
k = c(1, 10, 30, 50),
rep = 1:3,
{
locs <- sample(1:nrow(durations), n, replace = FALSE)
durs <- durations[locs, locs]
bench::mark(CBC = solCBC(durs, k),
HIGHS = solHIGHS(durs, k))
}
)
saveRDS(benchres3, "benchres3.rds")
} else benchres3 <- readRDS("benchres3.rds")
benchres3DF <- with(benchres3, data.frame(n, k, rep, median = as.numeric(median),
memory = as.numeric(mem_alloc)/1024/1024,
solver = bench:::as.character.bench_expr(expression)))
ggplot(benchres3DF, aes(x = n, y = median, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k) + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Időigény [s]", color = "Megoldó")
A most már nagyon eltérő időigények miatt érdemes az ábrát úgy is megnézni, hogy az egyes panelek függőleges tengelyének beosztása eltérhet:
ggplot(benchres3DF, aes(x = n, y = median, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k, scales = "free") + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Időigény [s]", color = "Megoldó")
A helyzet most már nem egyértelmű, de egy leheletnyit a HiGHS tűnik jobbnak.
Nézzük meg a memóriaigényt is:
ggplot(benchres3DF, aes(x = n, y = memory, color = solver, group = interaction(solver, rep))) +
facet_wrap(~k) + geom_line(alpha = 0.5) +
labs(x = "Városok száma", y = "Memóriaigény [MB]", color = "Megoldó")
Ebben gyakorlatilag semmilyen értékelhető különbség nem mutatkozik.
A fenti vizsgálatok egyik konklúziója nagyon egyértelmű: sem a GLPK, sem az lp_solve nem alkalmas nagyobb méretű problémák megoldására. Az, hogy a Cbc és a HiGHS közül melyik a jobb, már közel nem ilyen egyértelmű, egy nagyon kicsit talán a HiGHS tűnik jobbnak, így a továbbiakban azt fogom használni, de ismétlem, hogy nagy különbség minden bizonnyal nincsen.
Tulajdonképpen már csak egy módosítást kell tennünk: meg kell engedni az
általános esetet, amikor az
m <- nrow(durations)
n <- sum(JarasInd)
if(!file.exists("LocationResult.rds")) {
LocationResult <- lapply(1:50, function(k) {
tim <- bench::system_time(res <- highs::highs_solve(
L = c(c(t(durations[, JarasInd] * outer(HNTdata$Lakó.népesség, rep(1, n)))), rep(0, n)),
lower = rep(0, m*n + n),
upper = rep(1, m*n + n),
A = rbind(Matrix::sparseMatrix(i = rep(1:m, each = n), j = 1:(m*n), x = 1, dims = c(m, m*n + n)),
cbind(Matrix::Diagonal(m*n, 1),
Matrix::sparseMatrix(i = 1:(m*n), j = rep(1:n, m), x = -1)),
Matrix::sparseMatrix(i = rep(1, n), j = ((m*n)+1):(m*n + n), x = 1)),
lhs = c(rep(1, m), rep(-Inf, m*n), k),
rhs = c(rep(1, m), rep(0, m*n), k),
types = rep("I", m*n + n),
maximum = FALSE
))
list(time = tim[["real"]], obj = res$objective_value, sol = round(res$primal_solution))
})
saveRDS(result, "LocationResult.rds")
}Amint látható, elmentettem a futtatáshoz szükséges időt is (ez alapján készült a főszövegben a megfelelő ábra).
Mentsük ki a konkrét eredményeket is:
LocRes <- data.table(k = 1:length(LocationResult),
obj = sapply(LocationResult, function(x) x$obj/sum(HNTdata$Lakó.népesség)),
time = sapply(LocationResult, function(x) x$time))
LocRes <- cbind(LocRes, rbindlist(lapply(LocationResult, function(x)
data.table(t(HNTdata[JarasInd,][x$sol[(m*n+1):(m*n+n)]==1, ]$Helység.megnevezése))), fill = TRUE))
fwrite(LocRes, "LocRes.csv", dec = ",", sep = ";", bom = TRUE)