Szkocka krata i gra w życie - czyli sąsiedztwo Moora

Dotychczasowe przykłady automatu dwuwymiarowego stosowały tzw. sąsiedztwo von Neumana. Alternatywą jest sąsiedztwo Moore, które teraz zaimplementujemy w naszym automacie z regułą modulo:

Różnica polega na tym, że sąsiadami komórki aktualnej, poza sąsiadującymi "przez ścianę", są też komórki sąsiadujące "po rogach". Czyli zamiast 4 sąsiadów komórka ma ich 8.

Reszta programu pozostaje bez zmian, ale efekt mniej już przypomina arabskie mozaiki, a bardziej jakieś tkaniny - np. momentami serwetki, a czasem szkocką kratę...

Mając już taką wersję automatu komórkowego jesteśmy dosłownie o krok od zaimplementowania najsławniejszego z automatów - "Life" czyli "gry w życie" Johna Conwaya.

Po pierwsze deklarujemy trzy zmienne (linie 12-14), które są parametrami modelu - minim to najmniejsza liczba sąsiadów przy których komórka przeżywa, czyli jej stan pozostaje równy 1, maxim to maksymalna liczba sąsiadów przy których komórka przeżywa, a birth  to liczba sąsiadów, przy której komórka "martwa" (o stanie 0) może stać się żywa, czyli posiąść stan 1. Wartości domyślne tych zmiennych odpowiadają parametrom Life Conwaya, ale modyfikując je możecie zbadać też niektóre inne automaty z tej rodziny - czasem są zaskakujące.

Po drugie, ze względów kosmetycznych staramy się żeby szerokość i wysokość świata automatu była wielokrotnością 3 (linia 6) i  wyświetlamy sobie kontrolnie na konsoli realny rozmiar świata i okna, żeby mieć pewność że okno ma wystarczający rozmiar.

Wreszcie zmieniamy inicjalizację automatu w taki sposób żeby był zasiewany wyłącznie jedynkami - bo takie są reguły tego modelu.
Po modyfikacji deklaracji i procedury setup() możemy przejść do implementacji reguły automatu w procedurze draw():

Sama reguła (linie 87-98) zależna jest obowiązkowo od stanu komórki aktualnej. 

• Jeśli jest ona martwa to może ożyć wtedy i tylko wtedy gdy ma birth sąsiadów. 
• Gdy jest już żywa, to może przeżyć jeśli jej liczba sąsiadów mieści się w zakresie od minim do maxim.

W obu przypadkach zapisujemy obie wartości do tablicy WorldNew, bo nie zawiera ona stanu poprzedniego tylko poprzedni od poprzedniego. Można by ten fragment zoptymalizować, ale chyba nie warto...

Można też zredukować liczbę kolorów w pętli wyświetlania co nieco przyśpieszy wykonanie instrukcji switch. Zrobiłem to (linie 54-58), ale nie sprawdzałem jaki jest zysk, bo i tak wąskim gardłem tego programu jest przede wszystkim wyświetlanie.
Poniżej efekt działania po dosyć dużej liczbie kroków, gdzie praktycznie całość świata zajmują już struktury niezmienne "martwe natury" i proste oscylatory...
Ale i tak najciekawsze jest obserwowanie dynamiki - żadem obrazek tego nie odda. Obserwując uważne uda wam się na pewno zaobserwować tzw. szybowce, ale może też inne gatunki "statków kosmicznych"

Przyśpieszanie automatu 2D

Automat prezentowany poprzednio jest wystarczająco szybki, żeby go obserwować, ale gdybyśmy chcieli robić jakieś eksperymenty - np. sprawdzać okresowość, albo gdybyśmy zwiększyli rozmiar świata do 1000x1000 to byłoby prawdopodobnie już zbyt wolno. Przynajmniej na moim komputerze robiło się nużąco ;-)
W jaki sposób możemy ten program przyśpieszyć?
Prawdziwi programiści dokonali by tzw. profilowania kodu i ustalili, w których liniach program "siedzi" najdłużej. Wy na razie musicie mi uwierzyć na słowo - w tym przypadku najkosztowniejsze czasowo jest po prostu rysowanie na ekranie, czyli wywołania procedury point().
Co możemy na to poradzić skoro za każdym razem musimy odrysować stan całego automatu?
Ale czy na pewno?
Zauważcie że jeśli startujemy od pojedynczej komórki (co jak ustaliliśmy jest znacznie ciekawsze niż start z losowego wypełnienia większą ich liczbą) to widzimy wyraźną okresowość działania automatu. Bardzo rzadko faktycznie zmienia się stan wszystkich komórek - zazwyczaj spora cześć, zwłaszcza czarnych pozostaje w tym samym kolorze...
Gdybyśmy wiedzieli czy komórka na ekranie ma już właściwy kolor moglibyśmy rysować tylko to co się zmieniło, zmniejszając liczbę wywołań kosztownej operacji point().
Ale przecież wiemy! Potrzebna informacja jest już w programie - wystarczy porównać komórkę tablicy WorldOld z tablica WorldNew.

Pojawia się jednak problem kosmetyczny. Napis który dotychczas pojawiał się na co krok odświeżanym tle jest teraz na skraju pola, który odświeżany jest stosunkowo rzadko, więc wyniki działania kolejnych wywołań text() nakładają się na siebie.

Musimy więc przenieść napis gdzieś poza obszar wyświetlania świata i samemu zadbać o odświeżanie tła. To właśnie wykonywane jest w liniach 93-96.

Wszystko w porządku?
Nie?
Oczywiście że nie - chyba że sami wpadliście na konieczne modyfikacje setupu.

Po pierwsze okno musi mieć dodatkowe miejsce na wyświetlanie napisu. Będzie więc trochę prostokątne (linia 15).

Ale jest jeszcze jedna, poważniejsza sprawa. Musimy zadbać o to, żeby w pierwszym kroku CAŁA tablica New różniła się od tablicy Old.  Gwarantuje to wypełnienie tablicy New wartościami -1. Tablica Old jest przez system wypełniona zerami z pewną domieszką naszych wylosowanych stanów. Dzięki temu wszystko zadziała...

Czyli udaje nam się przynajmniej okresowo osiągnąć prędkość ok. 70 klatek na sekundę. Na waszych nowych komputerach może być nawet szybciej.

Ale bez przesady z tymi zachwytami ;-) Cudów nie ma i ta optymalizacja nie zawsze zadziała, a wręcz może być niepotrzebnym obciążeniem. Kiedy?
Wtedy gdy wystartujemy z dosyć gęstego zasiewu, np. Dens=0.1
W takiej sytuacji dosyć szybko osiągamy stan, w którym każda komórka zmienia się w każdym kroku i znowu trzeba odrysować każdą. 

Rozbudowa automatu 2D

Teraz rozbudujemy nasz dwuwymiarowy automat komórkowy tak, żeby zabawa była ciekawsza ;-)
Przy okazji porównamy szybkości naszych komputerów.
Na początek zmiany w setupie:

• Po pierwsze zwiększamy tablice, żeby automat miał więcej miejsca na rozwój. Tą liczbę warto lekko tuningować (zwiększając lub zmniejszając małe liczby) ponieważ każde ustawienie automatu ma swoją naturalną szerokość i jeśli obwód torusa (jakim jest świat naszego automatu) pokryje się z naturalną szerokością, to wynik będzie wyglądał tak jak w przestrzeni "nieskończonej".
• Po drugie dodajemy zmienne Div i self, analogicznie do tego co robiliśmy w jednowymiarowej implementacji reguły modulo. 
• Po trzecie poprawiamy losowanie (linia 24), dopasowując je do istnienia zmiennej Div. Musimy użyć rzutowania na (int) bo wynik funkcji random() w Processingu jest typu float co może niemile zaskakiwać osoby przyzwyczajone do C/C++ czy Pascala.
• Możecie też sprawdzić co by było gdyby wykomentować linię 16. U mnie efekt jest fatalny - odblokowanie antyaliasingu bardzo spowalnia program, nie mówiąc już o tym, że dla automatu komórkowego jest to jedynie źródło artefaktów graficznych.

Przejdźmy teraz do zmian w procedurze draw():

Tu też mamy kilka zmian, choć raczej kosmetycznych.

1. W pętli zmiany stanu automatu (linie 55-75) używamy nowych zmiennych Div i self. Div służy oczywiście jako dzielnik operacji modulo (linia 73) a self pozwala zadecydować czy aktualny stan komórki ma/nie ma wpływu na stan następny za pomocą znanej nam już operacji ?:   (linia 66).
2. W pętli rysującej dodaliśmy kolejne możliwe kolory konturów i pikseli (linie 42 i 43)
3. Napis w numerem kroku przesunęliśmy na dół okna zmieniając współrzędną y na wartość WorldSize. Uzupełniliśmy go też o wartość definiowanej przez Processing zmiennej frameRate, która podaje realną wydajność programu w postaci średniej liczby klatek na sekundę.

Ustawiając w procedurze setup() dużą wartość frameRate() wymuszamy na Processingu maksymalny "wysiłek", ale nie mamy żadnej gwarancji że realna liczba klatek będzie choćby zbliżona. W końcu tej roboty jest całkiem sporo - narysowanie każdego punktu to wywołanie procedury point(), a tych wywołań jest ponad ćwierć miliona na klatkę obrazu!
Zresztą sami zobaczcie i podzielcie się informacją w komentarzach na blogu.

Dla Dens=0.1 efekt graficzny nie będzie zbyt ciekawy, za to wynik liczby klatek na sekundę maksymalnie miarodajny.
Dla dalszych prac proponuję jednak wrócić do inicjowania pojedynczą komórką i pomanipulować parametrami modelu self, Div i WorldSize. 

Reguła modulo w DWUWYMIAROWYM automacie komórkowym

Automat dwuwymiarowy od jednowymiarowego technicznie różni się nieznacznie. W zasadzie wystarczy zamiana tablicy jednowymiarowej na dwuwymiarową (linie 5-6). Reszta jest w zasadzie techniczną konsekwencją tej operacji - jednowymiarowe odniesienia zamieniają się na dwuwymiarowe (np. w liniach 18, 22, 34), pojedyncze pętle zamieniają się na podwójne (linie 16-18 i 31-42), no a reguła automatu musi uwzględniać drugi wymiar...

Reguła sięga teraz do lewego i prawego sąsiada każdej komórki i te indeksy (zmienne right i left deklarowane w liniach 47 i 48) obliczane są w pętli zewnętrznej. Uwzględniamy też sąsiadów na górze i na dole, których indeksy obliczamy już w pętli wewnętrznej (zmienne dw i up

 w liniach 52 i 53).  Reszta jest analogiczna...
Dodatkowym elementem jest jeszcze wyprowadzenie na okno tekstu z numerem kroku za pomocą dawno w tych przykładach nie używanej "komendy" text();

Pozostaje nam do prezentacji wynik. Jednak całe działanie automatu 2D nie da się już ująć na pojedynczym obrazku. Potrzebny byłby film (co będzie innym razem). Poniżej więc tylko jedna, dokładnie pięćset szesnasta klatka tego filmu.

Całość zobaczycie, jak zaimplementujecie. Naprawdę warto ;-)

Wariacje na temat "reguły modulo" w automacie jednowymiarowym

Kolejna wersja rozwojowa automatu 1D pozwala nam w prosty sposób decydować o parametrach używanej reguły. Zmienna self (linia 28) jest typu boolean (czyli zmienna logiczna) i decyduje czy aktualny stan komórki wpływa na jej stan następny. Zmienna divider (linia 27) określa przez jaką liczbę dzielimy sumę w regule. Nie jest to jednak ta sama reguła co poprzednio, choć też używa operacji  reszty z dzielenia.

Tym razem zliczamy nie LICZBĘ sąsiadów, ale SUMĘ stanów grupy komórek. Zawsze zawiera ona komórki "skrzydłowe", natomiast aktualny stan komórki środkowej możemy pobrać albo nie - w zależności od wartości zmiennej logicznej self.  Używamy do tego charakterystycznej dla języków C-podobnych konstrukcji warunkowej ?: ( warunek ? wyrażenie dla prawdy : wyrażenie dla nieprawdy ) . 
Ostatecznie sumę stanów dzielimy przez divider i zapisujemy w odpowiedniej komórce tablicy WorldNew (linia 66).

Dla self==true i divider==5 uzyskujemy następujący wzór:

"Reguła modulo" w automacie komórkowym 1D

Przeprowadzamy kolejną modyfikację deterministycznego automatu jednowymiarowego, a ponieważ kod definicji zmiennych oraz procedury setup() jest identyczny jak poprzednio, skupimy się na samej procedurze draw():

Kod draw() zawiera kilka istotnych modyfikacji:

1. Do reguły wchodzi już nie tylko stan bezpośrednich sąsiadów komórki, ale tez sąsiadów odległych o 2 (linie 50-53) stąd poza zmienną right mamy jeszcze morer a zmienna left ma kuzynkę w postaci morel.
2. Może zostać też doliczony stan komórki aktualnej (wykomentowana linia 69)
3. Zamiast wykomentowanej reguły posiadania jednego sąsiada (linia 74) używamy reguły (linia 71) reszty z dzielenia liczby sąsiadów (lub komórek w obszarze o długości 5) przez jakąś stałą. W tym kodzie stała ta wynosi 3, ale warto sprawdzić też inne możliwości.
4. Zmieniona zostaje wizualizacja. Dzięki instrukcji switch (linie 35-42) możemy inaczej pokolorować komórki o stanach 0,1,2 oraz pozostałych.
5. Definiujemy osobne wyświetlanie aktualnego stanu w pasku na dole okna, czyli od y=994 do 999 (w linii 43). Dlatego jednak musimy ograniczyć liczbę wizualizowanych kroków automatu do 994 (linia 31)

Program ma potencjał na różnorodne eksperymenty - można używać lub nie stanu komórki aktualnej, można zmienić dzielnik, no i podobnie jak w poprzednim przypadku można manipulować gęstością początkową.
Oto jeden z wyników...

Deterministyczny automat komórkowy (i synchroniczny)

W poprzednich przykładach jednowymiarowych automatów komórkowych uaktualnienie następowało wg. "reguły Monte Carlo" - losowaliśmy jakąś komórkę, stosowaliśmy do niej reguły zmiany stanu, zmienialiśmy stan i losowaliśmy kolejną. Wykonanie N losowań, czyli takiej liczby losowań i jest komórek umownie nazwaliśmy "krokiem Monte Carlo". Jednak taki sposób uaktualniania nie tylko nie gwarantuje kolejności, ale nawet tego, że w danym kroku każda komórka zostanie uaktualniona. Wręcz przeciwnie - możemy być niemal pewni że w pojedynczym kroku M C będzie wiele komórek pominiętych i wiele wylosowanych dwa razy. A trochę takich wylosowanych trzy albo cztery - tym więcej im większe N.
Taki sposób uaktualniania bliski jest ekologom i naukowcom społecznym - bo wydaje się dobrze odzwierciedlać kolejność zdarzeń (a właściwie jej brak) w systemach zbudowanych ze względnie makroskopowych i obdarzonych "własną wolą" agentów.
Fizycy jednak wolą inną metodę uaktualniania - synchroniczną. W tym wypadku sprawdzamy wszystkie komórki "jednocześnie" - czyli wg. globalnego stanu aktualnego wyliczamy od razu globalny stan następny. Technicznie potrzebne są do tego dwie "tablice świata" - po jednej dla chwili obecnej i kolejnej. Uaktualnienie można zrobić za pomocą pętli po tablicy, lub w jakiś sposób je zrównoleglić (dlatego lubią to fizycy, bo zrównoleglanie algorytmów jest jedną z ich pasji ;-) )
Poniżej kod automatu dla tej samej co poprzednio reguły, ale z uaktualnianiem synchronicznym/równoległym - a zatem też w pełni deterministyczny. Każdy przebieg dla tego samego stanu początkowego będzie identyczny.

• W liniach 6 i 7 tworzymy dwie "tablice świata". W pętli (linie 36-54) w procedurze draw() przeglądamy KOLEJNO komórki tablicy WorldOld, ale wynik stosowania reguły zapisujemy w tablicy WorldNew, więc tablica WorldOld przez cały przebieg pętli pozostaje niezmieniona.
• Dopiero po zakończeniu pętli dokonujemy zamiany zawartości tablic (linie 57-59).
  To akurat w Processingu jest bardzo proste, bo każda tablica jest tu obiektem dostępnym za pomocą uchwytu, więc wystarczy zamienić to co trzymają uchwyty. Tak jakbyśmy zamieniali dwa ciężkie kanistry z wodą przekładając je z ręki do ręki...
  Musimy to jednak zrobić w powietrzu, więc potrzebna jest trzecia ręka do potrzymania na chwilę - to trzecią ręką jest zmienna WorldTmp.
• Inicjacja modelu (linie 13-19) łączy w sobie dwa dotychczasowe sposoby. Jeśli zmienna IDens jest większa od 0 to siejemy komórki losowo z zadaną gęstością, a w przeciwnym wypadku siejemy jedną komórkę w środku tablicy. Oczywiście zasiewamy tablicę WorldOld, bo to ona będzie następnie przetwarzana.
• Wyświetlanie (linie 29-34) też odbywa się z tablicy WorldOld, bo ona trzyma zarówno stan początkowy, jak każdy następny stan aktualny. Tablica WorldNew używana jest tylko w okolicy pętli zmiany stanu.
  
  Wyniki działania programu są zależne od gęstości, ale tak czy inaczej mało przypominają to co widzieliśmy w wersji probabilistycznej. Są znacznie bardziej regularne, co jest dosyć typowe. Fizykom i matematykom się podoba - pewnie dlatego że lubią regularność :-)

.01   .11.50