wtorek, 31 stycznia 2017

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 ;-)

niedziela, 29 stycznia 2017

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:

czwartek, 26 stycznia 2017

"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:

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.
Może zostać też doliczony stan komórki aktualnej (wykomentowana linia 69)
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.
Zmieniona zostaje wizualizacja. Dzięki instrukcji switch (linie 35-42) możemy inaczej pokolorować komórki o stanach 0,1,2 oraz pozostałych.
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...

poniedziałek, 23 stycznia 2017

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

niedziela, 22 stycznia 2017

"Zasiewanie" z zadaną gęstością (CA)

Wracamy dziś do jednowymiarowego automatu komórkowego (ang. Cellular Automaton, stąd skrót CA). Program opiera się niemal w całości na poprzednim przykładzie z jedną różnicą - zamiast pojedynczego ziarna automatu umieszczonego w środku "tablicy świata" wprowadzamy możliwość początkowego zasiania tablicy losowo z zadaną gęstością.

Sprowadza się to do jednego dodatkowego parametru modelu (IDens) oraz jednej pętli w procedurze setup(). Reszta kodu pozostaje niemal bez zmian. Czym się różni, zobaczycie, jeśli zaimplementujecie dokładnie ten powyższy kod.

wtorek, 17 stycznia 2017

Obliczanie Pi - przyśpieszone

A teraz obiecana przyśpieszona wersja programu obliczającego liczbę Pi metodą Monte Carlo:

Pierwotny program na moim laptopie wykonywał ok 1500 losowań na sekundę. Ten jest 1000 razy szybszy! Wykonuje około 1500000 losowań na sekundę. Składa się na to kilka modyfikacji:

W obrębie procedury draw() wstawiona jest pętla wykonująca pewną liczbę losowań (zmienna nst) w każdym wywołaniu. Dzięki temu zmniejszamy czas obsługi okna odbywającej się w związku z każdym wywołaniem draw(). Zamiast 1500 takich operacji mamy trochę ponad 70, przy nst == 20000. Czyli program losuje punkty i wrzuca je na wewnętrzną reprezentację grafiki okna, a okno odświeża nie 1500, ale jedynie 75 razy na sekundę. To daje gigantyczną oszczędność!
W procedurze setup() wyłączamy wygładzanie grafiki za pomocą noSmooth() .
Skoro rysujemy pojedyncze punkty to jest ono zbędne, a nawet niepożądane. Ta modyfikacja daje jeszcze około pięciokrotne przyśpieszenie
Zamiast losować składową koloru i odcień szarości dla wizualizowanych punktów używamy reszty z dzielenia (operator %). Jest to operacja dużo szybsza niż uzyskanie kolejnej liczby pseudolosowej, a efekt wizualny jest nieodróżnialny - i tak losujemy pozycje punktów, więc to że kolory mają sekwencyjne odcienie jest nie do zauważenia. Przy okazji oszczędzamy też liczbę wywołań generatora pseudolosowego co przy miliardach nawrotów pętli może mieć znaczenie (ale to temat na osobny post, a nawet kilka wpisów)
Tak znaczące przyśpieszenie programu pozwala osiągnąć w rozsądnym czasie kilka miliardów losowań. I tu pojawia się problem, bo w zakresie liczby typu int (32 bity ze znakiem) mieści się zaledwie trochę ponad 2 miliardy. A potem, czyli przy kolejnym powiększeniu o 1, liczba ta się przewija i staje UJEMNA! Nie ma sensu dalej liczyć. Stąd sprawdzamy w programie kiedy nastąpi ten moment i wyłączamy dalsze wywołania pętli draw() za pomocą noLoop();
Ale przedtem obliczamy uzyskaną wartość Pi i wyświetlamy w oknie graficznym czcionką o rozmiarze 18 pkt.

Oto wynik działania:

Tak oto uzyskaliśmy bardzo szybki program, aczkolwiek jego dokładność nie jest oszałamiająca. Wynika to z ograniczeń zastosowanych typów - float i int. Można by dokładność zwiększyć używając pojemniejszych typów - na początek także standardowe double i long. Ale to by rodziło kolejne problemy, bo są to typy języka Java nie do końca wspierane przez Processing.
Tak więc o ich stosowaniu będzie innym razem.

poniedziałek, 16 stycznia 2017

Obliczanie liczby Pi metodę Monte Carlo

"Metoda Monte Carlo (MC) jest stosowana do modelowania matematycznego procesów zbyt złożonych, aby można było przewidzieć ich wyniki za pomocą podejścia analitycznego.
Istotną rolę w metodzie MC odgrywa losowanie (wybór przypadkowy) wielkości charakteryzujących proces, przy czym losowanie dokonywane jest zgodnie z rozkładem, który musi być znany."

Oryginalny przykład w Wikipedii jest napisany w C++:

#include <iostream>
#include <cmath>
#include <ctime>
#include <cstdlib>
using namespace std;
int main()
{
     srand(time(NULL)); //zainicjalizowanie maszyny generujacej liczby pseudolosowe
     int n;
     int nk = 0;
     double x,y;
     float s;
     
     cout << "Podaj liczbe losowanych pkt:" << endl;
     cin >> n;

     for(int i = 1; i <= n; i++)
     {           
         x = ((double)rand() / (RAND_MAX)) * 2 - 1;
         y = ((double)rand() / (RAND_MAX)) * 2 - 1;
         if(x*x + y*y <= 1)
         {
             nk++;
         }
     }
     
     cout << "Liczba pkt. w kole wynosi: " << nk << endl;
     cout << "Liczba pkt. w kwadracie wynosi: " << n << endl;
     s = 4. * nk / n;
     cout << "Liczba pi wynosi: " << s;
}

Jednak przerobienie go na Processing trwa chwilę i gdyby nie dodanie wyświetlania to kod byłby nawet prostszy.

Kilka elementów jest tu warte uwagi.

Po pierwsze zniknęła główna pętla, bo została zastąpiona przez niejawną pętle wywołań funkcji draw() (od linii 24). Liczbę wywołań tej funkcji życzymy sobie mieć jak największą, zatem w funkcji setup() deklarujemy frameRate na 10000. To jednak oczywiście tylko "pobożne życzenie" - nie sądzę żeby udało się to na komputerze jaki macie do dyspozycji
Ponadto znacznie uprościły się instrukcje losowania x i y (linie 26 i 27). Sam język daje nam do dyspozycji funkcję random(float), która może zwracać liczbę w zakresie 0..1, co w C++ uzyskiwaliśmy przez dzielenie rand()/RAND_MAX;
A obliczenia wartości Pi wykonujemy nie na końcu programu, lecz za każdym razem gdy operator naciśnie jakiś klawisz na klawiaturze. Służy temu funkcja obsługi zdarzenia keyPressed(). Wyniki jej działania pojawiają się nie w oknie programu lecz na konsoli tekstowej. W szczególności dzięki predefiniowanej w Processingu zmiennej frameRate możemy dowiedzieć się ile naprawdę punktów na sekundę oblicza nasz program.
Wreszcie korzystając z tego, że w przeciwieństwie do C++ mamy od razu do dyspozycji okno graficzne rysujemy też wylosowane punkty, a właściwie 1/4 z nich - te które mają obie współrzędne dodatnie. Dzięki użyciu liczby losowej w kolorach cały czas będzie widać postęp działania programu, także wtedy, gdy już każdy punkt okna zostanie co najmniej raz zamalowany.

Uruchomienie tego programu zajmie wam chwilę. Niestety uzyskanie sensownego przybliżenia liczby Pi będzie trwało znacznie dłużej. Ten program nie grzeszy szybkością. Głównym ograniczeniem jest liczba wywołań draw(). Na moim, dosyć starym komputerze nie udaje uzyskać więcej niż ok. 1500 wywołań draw() w ciągu sekundy.
Ale każde wywołanie draw() wiąże się z obsługą ekranu - co jest czasochłonne.
To daje nam furtkę do znaczącego przyśpieszenia obliczeń w tym programie.

SAMI WYMYŚLCIE JAK! ROZWIĄZANIE W NASTĘPNYM POŚCIE
POWODZENIA!

Reference - czyli co jest nam dane

W przeciwieństwie do C,C++ czy Javy język Processing jest zorientowany na grafikę, stąd większość jego elementów właśnie do grafiki się odnosi. Natomiast z punktu widzenia struktur czysto syntaktycznych jest w porównaniu z językami ogólnymi mocno zubożony - w szczególności w obszarze programowania obiektowego. Jednak zawsze można użyć klasy wykonanej bezpośrednio w języku Java i to zadziała.
Aktualna składnia jezyka (już w wersji 3) i zawartość jego biblioteki standardowej jest dostępna z tej listy:
https://processing.org/reference/