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.

Histogram liczb pseudolosowych (cz. II)

Mając już program wykonujący losowania i wizualizujący histogram trafień, możemy przetestować różne sposoby uzyskiwania rozkładów innych niż płaski.

 

W powyższym kodzie mamy jeszcze 8 różnych wariantów, z których możemy skorzystać. Poza bezpośrednim użyciem generatora z języka JAVA (linia 79), mamy kilka przykładów dwu podstawowych modyfikacji:

1. Mnożenie kilku liczb losowych z zakresu 0..1 (linie 80-82) dające zamknięty do zakresu 0..1 rozkład skośny. Tym bardziej skośny, im więcej liczb pomnożymy.
2. Uśrednianie kilku liczb losowych z zakresu 0..1 (linie 83-87) dające zamknięte w zakresie 0..1 rozkłady dzwonowe naśladujące na potrzeby symulacji rozkład normalny.

Przykład rozkładu dzwonowego był pokazany w poprzednim rozdziale, a jednego ze skośnych mamy poniżej:
 

No fajnie. Ale dlaczego robimy to w ten sposób? Przecież istnieją transformacje płaskiego rozkładu zmiennej losowej w inne rozkłady, np. transformacja Boxa-Mullera?
I na to kiedyś przyjdzie pora :-) To co tutaj pokazuje to najprostsze sposoby uzyskania rozkładów użytecznych do celów symulacji agentowych i gier, w których sam rozkład jest tylko środkiem do celu, a nie celem samym w sobie ;-)

Histogram liczb (pseudo-)losowych

W programowaniu gier czy symulacji, ale także w wielu bardziej specjalistycznych zastosowaniach potrzebujemy często sekwencji liczb losowych. Jednak uzyskiwanie prawdziwych liczb losowych wymaga specjalistycznego hardwaru, który nie jest rozpowszechniony, stąd programiści "od niepamiętnych czasów" zastępują je liczbami pseudo-losowymi uzyskiwanymi za pomocą "generatorów [liczb] pseudolosowych".
Każdy szanujący się język programowania posiada w standardowym wyposażeniu przynajmniej jeden najprostszy generator dający liczby "o rozkładzie płaskim", czyli teoretycznie równie prawdopodobne. Taki generator jest zazwyczaj funkcją biblioteczną o nazwie rand(), random() lub zbliżonej. Oczywiście w Processing nie jest ułomny i też takie funkcje o nazwie random() ma. Ta funkcja zwraca liczbę typu float z zadanego zakresu.
Ma też odziedziczoną po macierzystym języku JAVA funkcję biblioteczną Math.random() zwracającą liczbę typu double z zakresu 0..1

Na razie nie będziemy wnikać w szczegóły działania tych funkcji, natomiast zajmiemy się sprawdzeniem co to znaczy że generator daje rozkład płaski.
Czyli zrobimy program który za pomocą histogramu wizualizuje wyniki (pseudo-)losowania z takiego generatora.



Generator z Processingu ukrywamy w funkcji MyRandom() (linie 4-8) po to,  żeby nie modyfikując już potem programu móc testować inne generatory. Zakładamy że funkcja ta zawsze będzie zwracać liczbę typu double z zakresu 0..1. Reszta programu to wizualizacja rozkładu uzyskanego z generatora, która dla każdego generatora spełniającego powyższe wymagania będzie taka sama.
Także dla takiego jak poniżej, który ewidentnie nie daje rozkładu płaskiego.

Zatem musimy zadeklarować tablicę która będzie służyć za koszyki histogramu (nazwaną Basket[]) oraz zmienne reprezentujące liczbę koszyków i zliczające losowania (linie 10-12). Tablicę dla porządku i przypomnienia sposobu użycia zerujemy w funkcji setup() (linie 23-24).

Potrzebujemy też określić częstość wywoływaniu funkcji draw() (linia 15 oraz 22) oraz liczbę losowań wykonywanych w każdym wywołaniu draw(). Te zmienne wpływają tylko na szybkość działania programu, ale nie zmieniają końcowego rezultatu.





Odpowiednio do tych ustawień implementujemy więc funkcję draw() :


Główna pętla (w liniach 30-45) wykonuje odpowiednią liczbę losowań (linia 32). Każda wylosowana liczba jest sprawdzana pod względem spełnienia założeń, czyli znajdowania się w zakresie 0..1. Następnie mnożąc wylosowaną liczbę przez liczbę koszyków ustalamy indeks i koszyka do którego liczba "wpada" (linia 39 i 43). Tu czyha na nas jednak pułapka. Jeżeli trafimy na rzadką sytuację że funkcja MyRandom() zwróci dokładnie 1 to wynik będzie dokładnie równy NumOfBaskets. Podobnie gdy generator nie będzie spełniał założeń i da nam wynik większy od 1. Na takie wypadki mamy w tablicy Basket[] dodatkowy koszyk (zobacz linie 11), do którego wpadają wszystkie takie wyniki losowań. 

• W normalnych warunkach zabezpieczenie (linie 40-41) nie jest potrzebne, ale kiedyś się jeszcze ucieszymy że je zrobiliśmy :-)

 Na koniec funkcji draw() czyścimy okno i rysujemy histogram (linie 47-48), ale stosując się do paradygmatu proceduralnego to ostatnie zadanie zlecamy kolejnej funkcji:

Początki programowania obiektowego

Dawno się do was nie odzywałem... Cóż, "taki life" :-) Zwłaszcza że ten rozdział będzie trudny i nie mogłem się zebrać do jego napisania...
Ale wiosna przyszła i pora wrócić do nauki.

Wchodzimy dziś w zupełnie nowy obszar, albo jak mówią informatycy "zmieniamy paradygmat". Poprzednio trzymaliśmy się programowania strukturalnego i  programowania proceduralnego, choć tak naprawdę było w tym drobne oszustwo wynikające z głębokiej semantyki języka Processingu.
Teraz zaczniemy jawnie korzystać z możliwości programowania obiektowego.

Co to oznacza? Otóż zaczniemy rozumieć, że używamy predefiniowane obiekty,  a także nauczymy się je sami definiować. 
Czym są obiekty? Pewnie jeszcze nie skorzystaliście z linków do Wikipedii to wam ją przytoczę :-)

• "Obiekt to podstawowe pojęcie wchodzące w skład paradygmatu programowania obiektowego w analizie i projektowaniu oprogramowania oraz w programowaniu. W tym paradygmacie programy definiuje się za pomocą obiektów – elementów łączących stan (czyli dane, nazywane najczęściej polami) i zachowanie (czyli procedury, tu: metody). Obiektowy program komputerowy wyrażony jest jako zbiór takich obiektów, komunikujących się pomiędzy sobą w celu wykonywania zadań."

Nie będę już dalej cytował teorii - ale zanim przejdziemy dalej zapoznajcie się z nią używając Wiki lub "wujka Google".
Przejdźmy do naszego nowego przykładu. Od razu mamy w nim obiekty

 

Mamy w nim od razu deklaracje dwóch obiektów. World jest tablicą, a tablice w Processingu są obiektami. Podobnie output ,  który jest narzędziem do zapisywania pliku tekstowego na dysku.

Mamy też tajemnicze RGB - "klasę" czyli "typ",  o którym napisano że dopiero zostanie zdefiniowany.  Tablica i PrintWriter są obiektami predefiniowanymi, które dostajemy w ramach języka Processing. RGB to będzie właśnie nasz własny typ obiektowy służący do przechowywania koloru rozbitego na poszczególne składowe. Zdefiniujemy go za chwilę.

Po czym to poznajemy że jakiś typ jest obiektem? Najłatwiej po sposobie użycia. Żebyśmy mogli używać tablicy World musimy najpierw użyć operatora new. Czyli w przeciwieństwie do int czy float samo zadeklarowanie nie wystarcza.
Podobnie będzie z typem elementu tablicy World. W przeciwieństwie do poprzednich naszych programów, gdzie elementy tablicy były albo typu int albo float, RGB nie jest typem prostym, ale klasą, czyli typem obiektowym. A skoro go sami wymyśliliśmy, to musimy zdefiniować (czyli powiedzieć dokładnie) czym jest:

Składnia takiej deklaracji polega na zamknięciu w strukturę class Name {   }

deklaracji danych i deklaracji funkcji/procedur. Te dane opisują stan obiektu, procedury opisują możliwe sposoby zmiany tego stanu. Dane nazywane są też atrybutami, a funkcje i procedury  metodami.
W naszym przykładzie mamy trzy atrybuty każdego obiektu RGB, trzy liczby typu int nazwane R, G i B. Mogłyby się nazywać też Red, Green, Blue, albo jeszcze inaczej (np. SkladowaCzerwona, SkladowaZielona, SkladowaNiebieska), ale po co komplikować sobie zapis?
Poza atrybutami mamy też kilka metod: Visualise(), isEmpty(), Set() oraz RGB().
Procedura Visualise() ma za zadanie narysowanie graficznej reprezentacji koloru w postaci punktu lub kwadratu w miejscu ustalonym przez współrzędne XY. Współrzędne odnoszą się jednak do położenia w tablicy World , a nie do okna, bo okno może być większe od tablicy dzięki "mnożnikowi" czyli współczynnikowi skali W.
Wyświetlanie odbywa się jednak tylko wtedy gdy obiekt RGB nie jest pusty, co jest sprawdzane specjalną funkcją zwracającą wartość logiczną (boolean).
Funkcja isEmpty() sprawdza czy obiektowi RGB nadano jakiś kolor inny niż czarny - czyli czy któraś ze składowych różni się od 0. 

Do zmiany wszystkich składowych "za jednym zamachem" służy procedura Set().
A do czego służy procedura czy też funkcja o nazwie RGB czyli takiej samej jak cały typ?
To jest tzw. konstruktor. Jest używany przez operator new do zainicjowania obiektu. Powinien zatem nadawać mu jakiś ściśle zdefiniowany stan. Nasz używa instrukcji R=G=B=0; żeby wszystkim składowym nadać stan 0. Konstruktor inkrementuje też zmienną RGB_Counter, dzięki czemu wiemy ile stworzyliśmy już obiektów typu RGB.
Jeszcze jedno. Definiując konstruktor nie podajemy typu jaki zwraca. Taka jest konwencja we wszystkich językach obiektowych jakie znam.  

Skoro typ RGB już mamy zdefiniowany, pora na jego użycie. 

Że RGB jest obiektem, widać już w procedurze setup. Żeby zacząć używać komórki tablicy World nusimy w niej najpierw utworzyć obiekt typu RGB za pomocą operatora new (linia 51).
Potem za pomocą metody Set() nadajemy mu jakiś kolor (linia 54, albo alternatywne 52 i 53 jeśli zdecydujemy się inaczej zakomentować), a za pomocą metody Visualise() rysujemy w środku okna (linia 55), jako że jest to akurat środkowy element tablicy.
Zwróćcie uwagę jak dostajemy się do tych metod.
Indeksowanie tablicy pozwala nam wybrać element tablicy. Na razie bezpiecznie możemy wybrać tylko obiekt w położeniu [Side/2][Side/2] bo tylko on został utworzony. Ale co ma zrobić obiekt? Do tego wyboru służy "kropka" (czyli operator wyboru składowej).
Po kropce możemy w Processingu napisać nazwę dowolnej metody, a także dowolnego atrybutu obiektu (w innych językach obiektowych sprawa jest nieco bardziej skomplikowana).

Nieco podobnie ma się sprawa z typem PrintWriter, który też występuje powyżej. Jest co prawda tworzony nie za pomocą operatora new ale specjalnej funkcji bibliotecznej createWriter() , ale operator new musi zapewne znajdować się wewnątrz tej funkcji.
Natomiast użycie odbywa się podobnie do pisania na konsolę tekstową Processingu, jednak właśnie z użyciem "operatora kropki". Pisząc output.println() wybieramy spośród różnych metod możliwych dla obiektu output typy PrintWriter operację println().

No to teraz pora na deser. Użyjemy sobie naszego typu do jakiejś ładnej symulacji definiując procedurę draw():



Procedura draw() wrzuca najpierw do pliku (linia 70) numer kroku i liczbę istniejących obiektów RGB. Pozwoli nam to wykonać w arkuszu kalkulacyjnym albo naszym własnymgramie wykres wzrostu. Właściwe działanie procedury odbywa się zaś w bloku dla warunku if(!Stop) - czyli dopóki nie nadamy zmiennej Stop wartości true.

W obrębie tego bloku powinniście rozpoznać znany już początek pętli kroku Monte-Carlo (Linie 75-79). Reszta działań odbywa się tylko wtedy gdy wylosowany element tablicy zawiera obiekt typu RGB co sprawdzamy za pomocą warunku World[Y][X] != null . Dla takiej sytuacji losujemy współrzędne Xt i Yt jakiejś sąsiedniej komórki w obrębie sąsiedztwa Moora (linie 82-83) i jeśli jest ona realnie w tablicy i jest pusta (warunek w liniach 84-85) dokonujemy jej zasiedlenia przez nowy obiekt, który będzie "potomkiem" obiektu znajdującego się w komórce Y X.
Potomek nie jest identyczną kopią, ponieważ z dużym prawdopodobieństwem może zajść mutacja jego składowych koloru (w liniach 89-94) dokonana metodą błądzenia losowego i dopiero takie zmutowane składowe są przypisywane na potomka (linia 96).
Na koniec sprawdzamy czy nie zasiedliliśmy właśnie jednego z czterech rogów macierzy, czyli punktów najbardziej odległych od środka (w programie sprawdzamy tylko jeden róg, ale nic nie stoi na przeszkodzie żebyście ten warunek uzupełnili). Jeśli doszliśmy do rogu to dalsze działanie modelu uznajemy za niecelowe i ustawiamy zmienną Stop na true.

A teraz małe wyjaśnienie w ramach post scriptum. Napisałem na początku że z paradygmat u obiektowego korzystaliśmy właściwie cały czas, choć o tym nie wiedzieliśmy. To fakt. Processing jest "nakładką" na język Java, a ten jest w 100% obiektowy. Cała aplikacja w Processingu jest więc obiektem języka Java, którego jedną, dwie lub trzy metody redefiniujemy (czy też  nadpisujemy). Te metody to draw(), setup() i ewentualnie exit(). Rozwiniemy tą kwestię przy okazji "dziedziczenia".

Model wpływu czy model Isinga

Prezentowany poprzednio dynamiczny model wpływu społecznego Nowaka-Latane bardzo przypomina komórkową implementacje bardzo słynnego w fizyce statystycznej modelu Isinga.
Możemy go jeszcze bardziej zbliżyć do wersji fizyków dodając szum (Noise), czyli zupełnie spontaniczne zmiany "poglądu", a także możliwość zamiennego stosowania "reguły większości" ("majority rule") i odwrotnej "reguły mniejszości" ("minority rule"), w której agent zawsze chce mieć pogląd mniejszościowy (np. dlatego, że to wyróżnia z tłumu). Polega to na odwróceniu znaków wyjścia reguły za pomocą zmiennej MajorityRule równej 1 lub -1, pod warunkiem że poglądy są zdefiniowane także jako 1 lub -1 (patrz. setup() )


Ponadto użyjemy też zapisu statystyki symulacji do pliku tekstowego rozdzielanego tabulacjami, co pozwoli wykonać wykres zmian w czasie w Excelu, Calcu lub innym programie czytającym takie dane.
Ponieważ dane do takiego pliku są przez program najpierw wpisywane do buforu w pamięci, a na dysk dopiero jak zgromadzi się ich odpowiednio dużo, plik musi zostać prawidłowo zamknięty w momencie kończenia programu. Dlatego definiujemy też procedurę exit(), wywoływaną automatycznie przez system przy każdym poprawnym zakończeniu w której właśnie znajdują się wywołania sprzątające i zamykające plik wyjściowy (flush() i close() );

Pozostała nam procedura Count() wywoływana w każdym draw() w celu zliczenia liczby "czerwonych" agentów, czyli tych których stan jest równy 1, oraz także wywoływana z draw() właściwa procedura symulacji DoMonteCarloStep().
W porównaniu z wyjściowym programem modelującym tylko wpływ mamy dwie ważne modyfikacje. Po pierwsze sprawdzamy czy "nie wdał nam się szum" (linia 81) losując liczbę z zakresu 0..1 i sprawdzając czy jest mniejsza od zadanego poziomu szumu. Jeśli tak, to po prostu zmieniamy stan agenta na przeciwny (linia 83), a właściwą regułę modelu stosujemy tylko wtedy gdy dany agent w tym kroku nie został trafiony szumem (linie 87-99).
W tej regule po zsumowaniu wpływów ustalamy odpowiedź agenta, ale mnożąc 1 i -1 przez zmienną MajorityRule. Jeśli jest ona równa 1 to wynik się nie zmienia, ale jeśli jest równa -1, o efekt działania reguły zostaje odwrócony.
Poniżej wynik działania programu dla małych kilkuset kroków "reguły większości". Wynik dla reguły mniejszości was zaskoczy :-)