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

Wpływ społeczny - reguła większości

Przyszła pora zaimplementować pierwszy prawdziwy model symulacyjny z psychologii społecznej. Będzie to najprostsza forma modelu dynamicznego wpływu społecznego w postaci probabilistycznego automatu komórkowego 2D.
Komórki automatu (w modelach psychologicznych częściej zwane agentami, ale w tym przypadku byłoby to ewidentnie "na wyrost") posługują się bardzo prostą regułą:
"Dopasowujesz swój stan -1 lub 1 do stanu WIĘKSZOŚCI twoich sąsiadów w sąsiedztwie Moora" co jest implementacją psychologicznego prawa konformizmu:
"Jednostka dopasowuje swoją postawę do postawy dominującej w jej otoczeniu".

Jak widać na powyższym obrazku symulacja ma cyklicznie (w procedurze draw() ) do wykonania trzy zadania: wizualizację, obliczanie statystyk (w Count() ) oraz właściwy model (DoMonteCarloStep()).

• Wizualizacja przegląda wszystkie komórki i te w stanie 1 rysuje w kolorze czerwonym, a inne, zakładamy że w stanie -1 w kolorze białym. Tym razem nie posługujemy się jednak pojedynczymi pikselami, ale prostokątami (rect(x,y,w,h)), a właściwie kwadratami bo oznaczające szerokości i wysokość parametry trzeci i czwarty są równe S. Także pozycjonując komórki posługujemy się tym parametrem, dzięki czemu ich "rogi startowe" oddalone są o wielokrotności S.
• Statystyka jest na razie bardzo prosta - jedynie zliczamy liczbę agentów w stanie 1 czyli czerwonych.
• Wreszcie sam krok modelu polega na wykonaniu tzw. "kroku Monte Carlo" czyli wykonaniu tylu losowań ile jest komórek. Dla każdej komórki sprawdzamy w podwójnej pętli stan jej i jej sąsiedztwa Moora sumując stany komórek. Ponieważ jest to 9 komórek,  to w sumie wychodzi nam albo liczba dodatnia albo ujemna i znak tej liczby staje się nowym stanem wylosowanej komórki centralnej. Warto zwrócić uwagę, że przy obliczaniu prawdziwych współrzędnych komórki zakładamy że świat symulacji zamknięty jest w torus (temu służą operacje modulo %)

Oczywiście program jest jeszcze niekompletny. Potrzebujemy tradycyjnie ;-) deklaracji i inicjalizacji zmiennych w procedurze setup():

Główne zadanie tej części programu polega na utworzeniu tablicy świata (linia 12) i jej zainicjowaniu liczbami -1 i 1 w proporcji zadanej parametrem Ones.
Oczywiście musimy ustalić też rozmiar okna, ale ze względu na sposób wizualizacji nie musi być on równy rozmiarowi tablicy. Wręcz przeciwnie - lepiej żeby był większy. Do wpasowania wizualizacji tablicy świata w okno służy współczynnik S, który obliczamy dzieląc szerokość okna przed długość wiersza tablicy (w linii 18).

Na koniec wprowadzamy jeszcze funkcję exit() obsługująca zakończenie działania programu. Na razie nie robi ona nic istotnego, poza wywołaniem swojego super odpowiednika i powiedzeniem "Do widzenia", ale w dalszych częściach kursu będzie nam się przydawać.

Nie pokaże na razie zrzutu ekranu - sami się przekonajcie co wychodzi :-)
I pomyślcie jak przerobić ten program, żeby dopuszczalny był też stan 0 - "komórka jeszcze nic nie wie, więc nie nie ma zdania".
Więcej o teorii związanej z tym modelem w poblikacjach naukowych profesora Nowaka, ale jest też przystępna wersja polska: WPŁYW SPOŁECZNY JAKO MODEL ROZPRZESTRZENIANIA SIĘ MEMÓW

Processing vs. Java

Processing jest ściśle związany z Javą, podobnie jak C++ jest związane z C.
Można dzięki temu w miarę łagodnie i etapami przejść od programowania w Processingu do programowanie w Javie.
Można najpierw programować Processing w Eclipse - ulubionym IDE starych wyjadaczy Javy, potem w Eclipse programować w Javie, ale z użyciem bibliotek Procesingu, potem mieszając...
Aż wreszcie zostać z samą Javą...
A wtedy ofert pracy nie zabraknie.

http://www.instructables.com/id/Proclipsing-Using-the-Eclipse-IDE-for-Processing-p/

https://youtu.be/0cqSjyvA8EY

Walentynkowy wykres funkcji 2D

Z okazji walentynek zrobiliśmy wariację na temat dwuwymiarowej funkcji ( f(x,y) ).

Korzystamy z ogólnej idei krzywych "sercokształtnych", jednak zamiast rysować krzywą, badamy cały szeroki zakres x, y kolorując wyniki.

Program oparty jest na poprzedniej wizualizacji funkcji 2D więc nie prezentuje go w całości - brakujące elementy są takie jak poprzednio:

Główna różnica polega na rozbudowanej instrukcji wybierającej kolor, która może posługiwać się skalą, albo też tylko kolorować jakościowo różne obszary (większe od zera, zero i mniejsze od zera). Można śmiało powalczyć z kolorystyką, no i z innymi niż zaimplementowana funkcjami.