środa, 10 marca 2021

Dywan Sierpińskiego

Dywan Sierpińskiego – to fraktal otrzymany z kwadratu za pomocą podzielenia go na dziewięć (3x3) mniejszych kwadratów, usunięcia środkowego kwadratu i ponownego rekurencyjnego zastosowania tej samej procedury do każdego z pozostałych ośmiu kwadratów. Nazwa pochodzi od nazwiska Wacława Sierpińskiego.


Co to dla nas oznacza? Jak się to ma do zbioru Cantora?

Zapewne widzicie, że taki dywan to coś do zbioru Cantora bardzo podobnego, tyle że dwuwymiarowego. Zamiast wycinać odcinek o długości ⅓ wycinamy kwadrat o boku ⅓ .

Odpowiednio do tego musimy więc zmodyfikować nasze procedury:


W programie brakuje kilku wywołań rekurencyjnych (od linii 28), ale myślę, że już sobie z tym poradzicie.

Ponieważ procedura rect() robi trochę roboty za nas, obliczenia są nieco prostsze niż dla poprzedniego ćwiczenia. Warto jednak przyjrzeć się, czy aby znowu nie jesteśmy trochę oszukiwani przez grafikę rastrową, zwłaszcza przy bardzo małych wartościach limitu.

Pośrednim etapem uzupełniania wywołań mogą być takie “ułomne” dywany jak poniżej:

Autor: o Brak komentarzy:

środa, 15 kwietnia 2020

Zbiór Cantora

Zbiór Cantora – podzbiór prostej rzeczywistej opisany w 1883[1] przez niemieckiego matematyka Georga Cantora. Zbiór ten odkrył w 1875 Henry John Stephen Smith.
Zbiór Cantora jest najprostszym przykładem fraktala.

Klasyczny zbiór Cantora (zwany także trójkowym zbiorem Cantora) to podzbiór przedziału domkniętego [0,1]. Jego konstrukcja polega na usuwaniu z odcinka jego środkowej 1/3 i aplikowaniu tej samej zasady rekurencyjnie do dwu pozostałych pod-odcinków. W świecie idealnym, po nieskończonej liczbie iteracji ;-) powstaje nam bardzo rozproszony podzbiór punktów z zadanego zakresu. W świecie grafiki komputerowej nie ma oczywiście sensu zajmować się czymś co jest poniżej rozdzielczości ekranu, musimy więc konstrukcję naszego zbioru zatrzymać na rozmiarze pojedynczego piksela.
Organizacja programu jest podobna do poprzedniego (testującego rekurencyjną procedurę rysowania linii kropkowanej) . Procedura setup() ustawia parametry okna i wywołuje testowaną procedurę rekurencyjną. Tyle że w tym wypadku dwie nieco odmienne.

Pokazana procedura cantorSetHor1() każdą iterację zbioru rysuje na innej linii ekranu, posługując się wartością d czyli długością odcinka dzielonego na danym etapie. Zbiór jest zdefiniowany w zakresie liczb rzeczywistych, więc posługujemy się ich przybliżeniem - typem float. Będziemy mieć z tym pewien problem, bo całość mapujemy na CAŁKOWITE współrzędne pikseli okna. Do tego dochodzi "skłonność" procedury line() do włączania końcowych punktów, co przy granicznych długościach powodowałoby asymetryczne nakładanie się rysowanych linii. Stąd jawnie zabieramy obu bocznym liniom punkty graniczące z linią środkową, czyli ta formalnie wycinaną, a praktycznie kolorowaną na zielono. Na wszelki wypadek kolorujemy też punkt środkowy odcinka. Będzie to widoczne tylko wtedy gdy linia zrobi się bardzo, bardzo krótka ;-)

Po wykonaniu rysowania (linie 26-33) wywołujemy znowu funkcję cantorSetHor1() dla prawego i lewego pod-odcinka (linia 35).

Alternatywna procedura cantorSetHor2() zbudowana jest niemal identycznie. Jedyna różnica to współrzędna Y okna używana w rysowaniu. Zawsze jest to height/2 , co powoduje że wszystkie efekty rysowania trafiają w tą samą linie. No i kolory zielony został zamieniony na 'cyan', a 'magenta' na czerwony.


Możemy się temu lepiej przyjrzeć, jeśli zmienimy w setup'ie grubość linii np. na 5. MUSIMY TEŻ WTEDY ZMIENIĆ  SPOSÓB KOŃCZENIA LINI!

strokeWeight(5); strokeCap(SQUARE);


Autor: o Brak komentarzy:

Rekurencja

"Rekurencja, zwana także rekursją (ang. recursion, z łac. recurrere, przybiec z powrotem) – odwoływanie się np. funkcji lub definicji do samej siebie.
W logice wnioskowanie rekurencyjne opiera się na założeniu istnienia pewnego stanu początkowego oraz zdania (lub zdań) stanowiącego podstawę wnioskowania (przy czym, aby cały dowód był poprawny, zarówno reguła, jak i stan początkowy muszą być prawdziwe). Istotą rekurencji jest tożsamość dziedziny i przeciwdziedziny reguły wnioskowania, wskutek czego wynik wnioskowania może podlegać tej samej regule zastosowanej ponownie."
"Rekurencja jest podstawową techniką wykorzystywaną w funkcyjnych językach programowania. Należy jednak zachować ostrożność przy używaniu rekurencji w rzeczywistych programach. Ryzyko istnieje szczególnie przy przetwarzaniu dużej ilości głęboko zagnieżdżonych danych."
Z Wikipedii
W praktyce rekurencja sprowadza się do wielokrotnego wywoływania funkcji przez samą siebie, coraz bardziej "w głąb stosu", dla coraz mniej skomplikowanego zadania, aż dochodzimy do trywialnej jego wersji i dalsze zagłębianie się nie ma już sensu. 
Rekurencja jest przez teoretyków programowania uznawana za samo sedno inteligencji algorytmicznej. Wg. nich nie jest programistą ten kto rekurencji nie rozumie ;-)
Ja nie byłbym tak ortodoksyjny, ale jednak uważam, że warto tą ideę wytłumaczyć, bo bez tego trudno o zrozumienie grafik fraktalnych, które są po prostu ładne :-)

Temat rekurencji wprowadza się zwykle za pomocą funkcji silni czy jakiś ciągów - np. Fibbonacciego. Ale to czysta matematyka, a ja w tym kursie staram się używać raczej grafiki, jako trafiającej do większej liczby odbiorców.
W grafice komputerowej raczej unika się algorytmów rekurencyjnych, jako mało efektywnych, ale istnieje jeden, historycznie istotny. To "wypełnianie przez sianie" czy też, wg. bardziej popularnej nazwy angielskiej "flood fill".
Niestety algorytm ten wymaga możliwości odczytania koloru punktu, który już jest na ekranie (czy w oknie), której jednak w Processingu nie znalazłem. Jest taka, która pozwala czytać kolor punktu załadowanego obrazka, co może kiedyś jeszcze wykorzystamy.

Dla ilustracji prostej graficznej procedury rekurencyjnej użyjemy procedury rysującej linię przerywaną. Jest ona oczywiście bardzo nieefektywna i do wymagających obliczeniowo zadań się nie nadaje, ale dydaktycznie powinna być wystarczająca.

Wynikiem działania programu będą czarne kropki naniesione wzdłuż dowolnego odcinka prostej. Ich gęstość może być różna. Na obrazki limit gęstości wynosi 10. Żółte tło pod kropkami stanowi "kontrolę" narysowaną zwykłą funkcją line()  Processingu.




Moglibyśmy te kropki oczywiście narysować też za pomocą pętli. Większość algorytmów rekurencyjnych da się zmienić na klasyczne, choć zapis rekurencyjny jest zazwyczaj dużo prostszy.
Oto kompletny program:


Funkcja setup() służy nam jak zwykle do ustawień, oraz do pierwszego wywołania naszej rekurencyjnej funkcji bline() (linia 12.) i narysowania kontrolnej linii klasycznym algorytmem (linia 11.). Funkcji draw() nie implementujemy bo jest w tym teście niepotrzebna.

Sama funkcja jest dosyć prosta. Najpierw sprawdzamy czy zadany odcinek jest wystarczająco długi żebyśmy chcieli postawić w nim kropkę (linie 17-20). Długość liczymy po prostu jako odległość Euklidesa między punktami.
Potem obliczamy współrzędne punktu leżącego pośrodku odcinka wyznaczonego przez parametry wywołania. Kolejność dalszych operacji jest dowolna. U mnie najpierw następują rekurencyjne wywołania dla obu połówek odcinka, a potem dopiero zaznaczenie punktu środkowego za pomocą point(), ale równie dobrze można zrobić odwrotnie.
Każde zatem wywołanie rysuje jeden punkt i próbuje narysować kolejne używając wywołania tej samej funkcji. Za zakończenie tego procesu odpowiada zmienna limit (linia 4.)

Jeśli zrozumiecie jak to działa, kolejne programy, rysujące znane fraktale będą dla was dużo bardziej zrozumiałe.
Autor: o Brak komentarzy:

środa, 25 marca 2020

Implementacja komórkowa modelu SIR

 Zaczniemy od NAJPROSTSZEJ WERSJI MODELU gdzie CHOROBA JEST BARDZO KRÓTKA i BARDZO MAŁO ZARAŹLIWA. Formalnie będzie to dwuwymiarowy, probabilistyczny (kroki MC) automat komórkowy z regułą SIR.

Setup odpowiada za zasiewanie tablicy z zadaną gęstością początkową zdrowymi komórkami, oraz jedną pojedynczą komórką zarażoną na środku.

Procedura draw() odpowiada za wizualizacje, i za zmianę stanu modelu. Zaczynamy od wizualizacji bo jest ona po prostu odziedziczona po poprzednich automatach komórkowych.

...

Ciąg dalszy w postaci PDF bo kopiowanie z Google docs na Google blogera okazuje się całkiem nie działać :-D 

Autor: o Brak komentarzy:

sobota, 21 marca 2020

Udostępnienie pierwszych rozdziałów książki o processingu na okoliczność epidemii koronawirusa

Zdecydowałem się udostępnić on line pierwszą część przygotowywanego skryptu (plik PDF), wraz z wcześniej już dostępnymi przykładami.
Jest to materiał przeznaczony dla osób, które jeszcze nigdy nie miały styczności z Processingiem, a może nawet nigdy jeszcze nie programowały.

Tak się śmiesznie złożyło, że to akurat 50 post na tym blogu 😋







Materiały dostępne są pod tym adresem:
https://github.com/borkowsk/sym4processing/tree/master/SimpleEdu

Skok do początku bloga! 

Skąd wziąć Processing?
Autor: o Brak komentarzy:

piątek, 20 marca 2020

SIR - Podstawowy model epidemii

Nazwa tego modelu to skrót od trzech stanów osoby zarażonej. S oznacza “susceptible” czyli „podatny”, I to “infected” czyli „zainfekowany/chory”, wreszcie R to “recovered” czyli “wyleczony” (i chwilowo odporny). Tą odporność możemy uznać za trwałą lub nietrwałą, przy czym z biologicznego punktu widzenia może to oznaczać dwie rzeczy:
  1.     Zanik “pamięci immunologicznej” dla danego drobnoustroju, co raczej zdarza się w normalnych warunkach rzadko (choć niektóre “czynniki zakaźne” wymagają kilku kolejnych szczepień, żeby odporność była trwała)
  2.     Mutacje drobnoustroju, powodującą że jego nowe szczepy uciekają przed odpornością nabytą już przez żywicieli. To znacznie częstszy przypadek, gdyż ewolucja bakterii, a zwłaszcza wirusów przebiega o kilka rzędów wielkości szybciej niż ewolucja człowieka.
Dla naszego prostego modelu mechanizm tego zaniku nie będzie jednak istotny. Bardzo ważne będą natomiast prawdopodobieństwa zarażenia (Infection), sposób wychodzenia zarażonych z populacji (wyzdrowienie, śmierć), i prawdopodobieństwo lub czas utraty odporności - o ile będziemy chcieli się tym zająć.
Ponieważ stany SIR możemy zakodować jako liczby całkowite, nasz model może być wykonany w konwencji probabilistycznego automatu komórkowego. Poniżej przedstawiłem bardzo ogólny algorytm. Występuje tam słowo “sieć”, ale nie przejmujcie się tym. Automat komórkowy to też rodzaj sieci - bardzo regularnej.
Zwróćcie uwagę że jeśli zamiast automatu synchronicznego użyjecie wersji Monte Carlo to krok uaktualnienia stanów automatu nie będzie poza wewnętrzną pętlą, tylko uaktualnienie będzie zachodzić po prostu w trakcie wykonywania pętli po losowych agentach. Pytanie też kiedy kończymy. Można do tego użyć jakiejś statystyki - np. Liczby chorych. Ale w Processingu możemy też kończyć po prostu gdy uznamy że już nic ciekawego się nie dzieje.
Macie wybór, jeśli chodzi o sposób uaktualniania modelu. Sugeruję jednak, żeby komórki wybierać metodą Monte Carlo. To bardziej realistyczne. Istoty żywe rzadko działają synchronicznie, a jeśli to robią, to dużym nakładem sił, używając odpowiednich mechanizmów synchronizacyjnych.
Liczba interakcji z sąsiadami komórki automatu może być różna - to dobry parametr modelu. Interakcja polega na „przekazaniu wirusa”, co jest możliwe tylko między aktualnie chorym, a aktualnie zdrowym osobnikiem. Reszta interakcji z punktu widzenia modelu jest „bezpłodna”.
Początkowy przyjmiemy dla ułatwienia że choroba trwa około jeden dzień, czyli do następnego wylosowania agenta.
Musimy mieć jednak “z tyłu głowy”, że czas trwania infekcji może być bardzo różny, więc powinniśmy mieć możliwość odliczania upływu czasu od początku infekcji agenta.
Upływ czasu może polegać na zmianie stanu komórki na liczbę dni od początku choroby i wreszcie zmianie stanu komórki chorej na stan odporności, po upływie zadanej liczby dni. To zgodnie z zasadą, że “katar leczony trwa siedem dni, a nie leczony tydzień”.
Alternatywnie możemy też przyjąć, nieco nierealistycznie, że w każdym dniu choroby mamy pewne prawdopodobieństwo wyleczenia (albo śmierci!) i wtedy czas trwania choroby będzie zmienny.

Teraz pomyślcie sami…
Może uda wam się zaimplementować model, zanim opublikuje moje rozwiązanie.

CDN




Autor: o Brak komentarzy:

czwartek, 19 marca 2020

Wstęp do modelu epidemii


“W języku potocznym termin epidemia używany jest jako synonim masowych zachorowań wywołanych chorobami zakaźnymi. Często definiuje się ją jako wystąpienie na danym obszarze zakażeń lub zachorowań na chorobę zakaźną w liczbie wyraźnie większej niż we wcześniejszym okresie albo wystąpienie zakażeń lub chorób zakaźnych dotychczas niewystępujących. Przypadki globalnych epidemii nazywa się pandemią.”
Oto kilka przykłady najważniejszych chorób epidemicznych (i pandemicznych):

  • grypy
    • grypa hiszpanka (1918-1919) – ponad 50 mln ofiar śmiertelnych na całym świecie
    • grypa azjatycka (1957) – ok. 1 mln ofiar śmiertelnych na całym świecie
    • grypa Hong-Kong (1968) – ok. 1 mln ofiar śmiertelnych na całym świecie
    • Pandemia grypy A/H1N1 (od 11 czerwca 2009) - ok. 12799 ofiar na całym świecie
  • AIDS – masowe zachorowania; zwłaszcza na kontynencie afrykańskim (wszystko od początku rozdziału do tego miejsca zaczerpnięte z tekstu “Modele epidemii“ ). Epidemia ta jest o tyle nietypowa, że wirus HIV bezpośrednio atakuje komórki układu odporności człowieka, upośledzające jego odpowiedź immunologiczną, a co więcej jest tzw. “retrowirusem”, co oznacza, że kopia jego kodu genetycznego jest włączana do DNA komórki, więc wirus ginie tylko wraz z zakażoną komórką, a jego produkcja może być tak niewielka, że zakażone komórki są dla układu immunologicznego nieodróżnialne od zdrowych.
  • Koronawirusy - SARS, MERS, SARS-2 czyli COVID-19
  • Stare choroby epidemiczne dorosłych, w większości zwalczone już higieną, szczepionkami i/lub antybiotykami: dżuma, czarna ospa (ostatni raz w Polsce w 1963 we Wrocławiu), cholera, kiła, gruźlica i trąd 
  • Opanowane choroby epidemiczne dzieci: polio, płonica (szkarlatyna), błonica (krup), krztusiec (koklusz), oraz jeszcze nie opanowane, mimo istnienia szczepionki: odra, świnka i różyczka.
  • Ospa wietrzna - jedna z bardziej zaraźliwych chorób wirusowych, ale na szczęście nie obarczona śmiertelnością. Chorują zarówno dzieci jak i dorośli, przy czym dorośli ciężej.
  • Chorobami epidemicznymi są też “katary”, czyli “nieżyty nosa” zwane też “przeziębieniami”. Wywołuje je wiele szczepów wirusów, które ze względu na długą koewolucję z człowiekiem nie są już niebezpieczne dla ludzi z normalnym systemem odpornościowym. W większości przypadków za powstanie kataru wirusowego odpowiadają rhinowirusy, wirusy paragrypy, wirus grypy i adenowirusy. Wielość wirusów wywołujących “przeziębienie”, ich względna nieszkodliwość i szybkie mutacje powodują, że przygotowanie szczepionek przeciw nim jest uznawane za nieopłacalne czy wręcz niewykonalne.
Różnorodne modele epidemii to jedna z ważniejszych użytkowo klas modeli matematycznych i komputerowych modeli symulacyjnych. Różnorodne modele epidemii chorób zakaźnych stosowane od lat w epidemiologii - niekiedy z niezłymi skutkami. Ostatnio coraz częściej są to modele oparte na sieciach społecznych – istotne szczególnie dla epidemii grypy, AIDS czy SARS (tym tematem zajmiemy się później). Zdarzają się także bardzo złożone, przypominające gry w SIMSy, modele epidemii w konkretnych miastach (np. Model grypy w Los Angeles), coś nieco zbliżonego możemy jeszcze zrobić ćwicząc w następnym rozdziale model epidemii w konwencji ABM.



EPIDEMIA jest także głównym modelem rozprzestrzeniania się informacji i innowacji w społeczeństwie, chociaż to nie do końca słuszne, i do czego też powinniśmy kiedyś wrócić.

CD w następnym wpisie.

Autor: o Brak komentarzy:

PROCESSING W EDUKACJI I MODELOWANIU

Dawno nie pisałem na tym blogu, co nie znaczy, że nic w tej kwestii nie robiłem.
Doraźne ciekawostki umieszczałem w tym czasie na stronie na Facebooku.
Powoli powstawała też książka pod analogicznym tytułem. Ma już nieco ponad 150 str. w Googlowym brudnopisie.
Jeszcze trochę i będę szukał wydawcy.

Pojawiła się też nowa wersja Processingu, a materiały pomocnicze do kursu umieściłem na GitHubie.

Poniżej aktualny spis treści powstającej książki. Może chcielibyście żebym umieścił w niej jeszcze jakieś tematy?


Autor: o Brak komentarzy:

środa, 26 lipca 2017

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ć.


Autor: o Brak komentarzy:

poniedziałek, 24 lipca 2017

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
Autor: o Brak komentarzy:

poniedziałek, 10 lipca 2017

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.



Autor: o Brak komentarzy:
Subskrybuj: Posty (Atom)