Realistyczny model wielokrotnych pożarów lasu cz. 2

Nie wiem ilu z czytelników podjęło własną próbę implementacji. Mam nadzieję że wielu, bo bez własnych prób programowania nauczyć się nie sposób.
Ale teraz moja implementacja - możecie zobaczyć na ile nasze myślenie podobnymi/odmiennymi drogami.

Po pierwsze, ponieważ mamy plik Log musimy zdefiniować procedurę exit() która ma głównie za zadanie zapisać bufor danych pliku na dysk, a potem zamknąć plik (linie 186-187):

Warto jednak wziąć pod uwagę, że deklarując taką procedurę "przykrywamy" istniejącą procedurę domyślną, która prawdopodobnie też wykonuje jakąś użyteczną pracę. Dlatego musimy tą starą procedurę wywołać, co odbywa się za pomocą super.exit()  (linia 189) .

No to przechodzimy do procedury doMonteCarloStep() zaczynając od ewentualnego zapalenia lasu. Najprostszym sposobem byłoby po prostu umieszczenie w głównej pętli instrukcji

 if(LigtP>random(0,1))
World[i][j]=coś tam...



Problem w tym, że taka instrukcja if z wywołaniem kosztownej funkcji random wykonywałaby się w każdym kroku Monte Carlo dla każdego żywego drzewa, a ze względu na to że LigtP jest, i powinno być bardzo małe, to bardzo rzadko warunek byłby spełniony.
W naszym programie oszukujemy więc trochę i bardzo przyśpieszamy ograniczając te losowania do niezbędnego minimum:

W każdym kroku M C doliczamy do zmiennej Burn sumaryczne prawdopodobieństwo zapłonu całego lasu (linia 76). Potem wykonujemy pętlę, która trwa tak długo jak zmienna Burn jest większa od 0 (linia 77). W pętli tej wykonujemy losowanie komórki świata, i JEŚLI TRAFIMY W DRZEWO to je podpalamy, co polega na obliczeniu czasu pożaru z dzielenia rozmiaru drzewa przez parametr FireTimeDiv. Wynik zmieniamy na ujemny co jest dla nas sygnałem, że drzewo już płonie. Jeśli podpalenie się uda to zmniejszamy zmienną Burn o jeden. No i ustawiamy "magiczną" flagę dla wizualizacji is_burning na true.

Teraz już możemy przystąpić do właściwej pętli kroku Monte Carlo:

Tak jak w poprzedniej wersji programu losujemy N*N komórek i w zależności od stanu wylosowanej komórki. Jeśli jest to pusta komórka (czyli o wartości 0) to próbujemy zasiać drzewo (linie 95-100).  Jeśli wartość jest dodatnia to po prostu drzewo trochę rośnie (linie 101-106). Wreszcie, gdy wartość jest negatywna to próbujemy zapalenia sąsiadów w sąsiedztwie Moora (linie 109-121) i dodajemy jeden do wartości komórki co razem z obcięciem części ułamkowej przy podpalaniu (linia 83 i 117) gwarantuje nam że każde płonące drzewo w końcu gaśnie stając się pustą komórką o wartości 0.  No i każde zapalenie drzewa ustawia nam flagę is_burning na potrzeby wizualizacji.

No i wreszcie dochodzimy do samej wizualizacji:

 

 ... której główna część prawie się nie zmieniła, poza tym że zliczamy puste, żywe oraz płonące komórki. Za to na koniec dodajemy wydruki statystyk (linia 164-178) na ekranie i zapis ich do pliku logu (linia 180-181).

No i to by było na tyle... Rezultat działania tego programu możecie znaleźć na Youtube, ale lepiej obejrzyjcie to sami. Efekty mogą was zaskoczyć.

Realistyczny model wielokrotnych pożarów lasu

Las rośnie niezwykle powoli, a płonie szybko. Wzrost lasu liczymy latami, jego pożar w godzinach. Jeden rok to 365.5 * 24 h czyli... A może policzcie sami. Na serwetce. Ostatecznie na kalkulatorze w komórce ;-)
To oczywiste, jednak stanowi spory problem dla kogoś, kto chciałby stworzyć model biorący pod uwagę oba te procesy...
Jest to jednak wykonalne i w tym rozdziale właśnie taki model zrobimy.

    

Będziemy w zasadzie modyfikować poprzedni program, ale zmian będzie na tyle dużo, że równie dobrze możecie zacząć od początku. Czyli od deklaracji zmiennych:

 

Najważniejsze  jak zwykle na górze - najpierw obliczamy sobie dwie wartości week i year przeliczone na godziny, które będą nam się później przydawać. Godzina będzie naszą jednostką czasu. Przyjmiemy że jeden krok Monte Carlo modelu to jedna godzina.

W porządniejszych językach programowania zrobilibyśmy z nich wartości const, ale w Processingu to nie działa.
Następnie "oczywiście" długość boku "macierzy świata" czyli N, co oznacza że będziemy mieli 900000 drzew. A dalej...  Będzie już mniej oczywiście.

Czas w jakim pali się drzewo zależy od jego różnych czynników (jak na przykład  gatunek drzewa czy wilgoLightPtność), których tu nie będziemy wprowadzać, choćby dlatego, że na razie zakładamy las jednogatunkowy - monokulturę (która wbrew pozorom może być też naturalna). To co jednak na pewno wpływa na ten czas to masa drzewa. Dla uproszczenia założymy że jest to zależność liniowa. Parametr FireTimeDiv mówi nam jak wielkość drzewa w naszych umownych jednostkach przelicza się na czas, w którym drzewo płonie i może zapalić sąsiednie. Wartość 50 oznacza że drzewo ważące 100 umownych jednostek pali się 2 godziny, a ważące 200 pali się 4 godziny (czyli kroki M C). Warto od razu zwrócić uwagę na deklaracje MatureT=220 określającą maksymalną dopuszczalną masę drzewa. Moglibyśmy użyć jakiegoś bardzie zaawansowanego modelu wzrostu niż liniowy, ale na razie liniowy z "tresholdem" musi nam wystarczyć

IgnitionP czyli prawdopodobieństwo zapłonu od płonącego sąsiada, oraz InitP czyli początkowa liczba to parametry znane z poprzedniej wersji modelu. Kolejne są już nowe. GrowS określa o ile umownych jednostek drzewo przyrasta średnio na godzinę. Raczej robi to bardzo powoli (0.0005), jak widać... Z kolei SeedP określa jakie jest prawdopodobieństwo że w danej godzinie na pustym polu wykiełkuje nowe drzewo. To też mała wartość bo nie powinna być większa niż "kilka na rok". No i wreszcie LightP - prawdopodobieństwo że w danej godzinie dane drzewo się zapali - np. od uderzenia piorunem. To już jest naoprawdę malutkie.

Następnie mamy dwie zmienne związane z wizualizacją. Znane wam już z poprzednich programów S oraz is_burning , którego ważna rola w optymalizacji wyświetlania zostanie wyjaśniona dalej.

Pozostałe zmienne - Step, empty, alives etc... posłużą nam do zbierania podstawowych statystyk modelu.
No i jeszcze Log. To uchwyt do pliku w którym będziemy skłądować dane statystyczne do późniejszej obróbki.
Teraz nieco zmodyfikowany  setup():

Zmieniamy frameRate() tak żeby zmaksymalizować liczbę kroków Monte Carlo obrabianych w ciągu sekundy. Chcielibyśmy żeby czas w modelu płynął co najmniej z prędkością miesiąc modelu na sekundę naszego czasu. Na niektórych komputerach będzie to trudne, ale na innych być może "wyciśniecie" nawet dwa miesiące. Wywołania noSmooth() i noStroke() maksymalnie upraszczają wyświetlanie, co jak już wiecie z wcześniejszych przykładów (zwłaszcza noSmooth() ) mocno przyśpiesza działanie programu.
Wreszcie tworzymy nazwę pliku zawierającą nasze główne zmienne kontrolne, tworzymy plik Log o takiej nazwie, i zapisujemy do niego pierwszy wiersz będący nagłówkami kolumn danych.
Procedura draw() będzie miała jedną, za to bardzo istotną modyfikacje, dzięki której będziemy mogli obserwować zarówno pożary jak i wzrost drzew, nie nudząc się przy tym zbytnio.

Jak widzicie w każdym wywołaniu draw() uruchamiamy doMonteCarloStep(), ale doVisualisation() uzależniamy od warunku (w liniach 63-64).
Ten warunek pozwala nam wyświetlać każdy krok gdy trwa pożar i zmienna is_burning jest ustawiona na true w poprzednim kroku Monte Carlo, ALBO (||) gdy nie ma pożaru obserwować jedną klatkę na "miesiąc" wzrostu drzew.
Tyle informacji powinno wam już wystarczyć do podjęcia własnej próby implementacji modelu.

CIĄG DALSZY W NASTĘPNYM WPISIE

Pożar lasu w Monte Carlo ;-)

Oczywiście w Monte Carlo czyli stolicy Monako lasów już nie ma. Na 1 km kwadratowym ledwo znalazło się miejsce dla "paru drzewek", a parkingi dla samochodów wydrążone są w skale. Ale ja dziś nie o tym chciałem... ;-)

Dziś będzie o modelu pożaru lasu ("forest fire"), którego klasyczna wersja będąca typowym automatem komórkowym pochodzi z końca lat 80tych XX wieku i jest jednym z przykładów dla

self-organized criticality oraz dla perkolacji.

Nasz model tym będzie się różnić od klasycznego, że zamiast synchronicznie stosować reguły, użyjemy algorytmu uaktualnienia Monte Carlo.

Powoduje to że w porównaniu z modelem klasycznym nasze drzewa muszą palić się nieco dłużej niż przez 1 krok czasu. Uznajemy że czas ten jest PROPORCJONALNY do wielkości czyli wieku drzewa.

Obok więc klasycznego parametru N (linia 6) czyli długości boku "macierzy świata" (World) mamy parametr FireTimeDiv , który określa czas "płonięcia" drzewa przez podzielenie jego wieku. Domyślnie ma on wartość 10, co możemy rozumieć tak, że drzewo stuletnie płonie 10 godzin (czyli 10 kroków M C modelu).
Ponadto mamy parametr InitT (linia 9) wskazujący jak gęsty jest las w porównaniu z lasem wypełnionym maksymalnie. Domyślnie 0.75 oznaczający że w lesie jest 75% możliwej maksymalnie liczby drzew.
Ostatni parametr to IgnitionP (linia 8) czyli prawdopodobieństwo zapalenia w danym kroku jednego z sąsiednich drzew przez drzewo, które już płonie. Wartości tego parametru nie przekładają się bezpośrednio na model klasyczny, w którym drzewo płonęło zawsze tylko jeden krok czasu.
W jaki sposób budujemy las? Ponieważ nie chcemy (na razie) modelować jego wzrostu zasiewamy go wg. jakiegoś rozkładu (linie 26-33). Możemy stworzyć monokulturę w której wszystkie drzewa mają po 100 lat (linia 29), ale możemy też użyć rozkładu płaskiego w którym tylko najstarsze drzewa mają 100 lat, średnia wynosi 50, i każdy wiek pomiędzy 0 a 100 lat jest równie prawdopodobny.
No prawie ;-) To co napisałem nie do końca zgadza się z tym co jest w programie. Co należy zrobić żeby linia 30 robiła dokładnie to?
Linia 31 tworzy nam las w którym wiek drzew ma rozkład zbliżony do normalnego, który najbardziej lubią statystycy i fizycy. Czy to lepiej odzwierciedla leśną rzeczywistość? Trochę wątpię, ale musimy zostawić to zagadnienie na później.

Procedura draw() służy nam tylko do wywołania dwóch typowych dla symulacji zadań - wizualizacji i uaktualnienia świata modelu.

Wizualizacja klasycznie już rysuje kwadraciki (choć łatwo by je było zastąpić kółkami), a może trójkątami czy nawet "choinkami") różniące się kolorem:

• Jeśli w tablicy World jest 0 to oznacza że jest ona pusta. Albo nigdy nie było tam drzewa, albo się ono już spaliło. 
• Liczba dodatnia oznacza żywe drzewo, którego intensywność zielonego koloru jest proporcjonalna do wieku. Choć nasze drzewa mają najwyżej 100 lat, zabezpieczamy się na przyszłość i nigdy drzewo nie jest bardziej zielone niż 255 (linia 55).
• Jeśli w komórce macierzy jest liczba ujemna to oznacza że drzewo właśnie płonie i jego kolor jest losową mieszanką składowej czerwonej i zielonej co daje nam głównie różne odcienie koloru żółtego.

A skąd się bierze pożar i jak przebiega? 

Cóż, pożar wywoła sam operator programu. Naciskając dowolny przycisk na klawiaturze aktywuje obsługę zdarzenia keyPressed(), a w tej obsłudze losowana jest jedna z komórek świata (linie 69 i 70) i jej wartość jest zmieniana na ujemną (linia 71), wynikającą z podzielenia wieku drzewa przez FireTimeDiv. Ponieważ jest to operacja przebiegająca na liczbach całkowitych, to drzewa młodsze niż FireTimeDiv spalałyby się w czasie 0 kroków. Odejmujemy więc jeszcze 1, czyli czas pożaru pojedynczego drzewa nie może być krótszy niż 1 krok M C.
Sam krok Monte Carlo przebiega w sposób następujący:

1. Ustalamy liczbę losowań M jako N do kwadratu.
2. Wykonujemy M nawrotów pętli (linia 78-100). Licznikiem pętli jest małe m, co jak najbardziej działa, jako że język Processing, tak jak praktycznie wszystkie języki C-podobne jest wrażliwy na wielkość liter (case sensitive) czyli zmienne M i m są różne. Ale nie polecam wam na przyszłość tej sztuczki. Pochodzi z arsenału konkursów na "obfuscated code" :-)
3. W pętli oczywiście losujemy kolejne komórki.
4. Jeśli komórka jest ujemna (test w linii 83) to 

• losujemy jednego z jej sąsiadów stosując znaną już metodę "zapinania" świata symulacji w torus za pomocą operacji modulo N (linie 86-87)
• jeśli wylosowany sąsiad jest żywym drzewem to z prawdopodobieństwem IgnitionP próbujemy go podpalić w taki sam sposób jak w przypadku pierwszego zapłonu (linie 89-94)
• Natomiast wartość komórki aktualnej powiększamy o jednostkę (linia 97) co gwarantuje nam że w końcu osiągnie 0 i przestanie się palić. Jeśli wszystko dobrze zrobimy to nigdy nie ujrzymy na konsoli znaku ? pochodzącego z testu kontrolnego (linia 98) .

Jak już zadziała zaprezentowany powyżej kod to proponuje się pobawić parametrami gęstości lasu (InitT) oraz prawdopodobieństwa zapłonu, które możemy rozumieć jako odwrotność wilgotności.
Zauważyliście już zapewne, że nasz las nie odrasta. Cóż... Las rośnie lata, a pożar trwa godziny. Pod tym względem klasyczny model "forest fire" jest BARDZO ODLEGŁY OD RZECZYWISTOŚCI.