2009-10-26
printpaper.org
Ponad miesiąc temu uruchomiłem nowy projekt: printpaper.org.
Główna funkcja już jest i działa. Mam jeszcze sporo pomysłów do zrealizowania (w tym wersje językowe), ale muszę się przy tym bardzo powstrzymywać, bo staram się utrzymać minimalny charakter tej usługi. Informacje o nowościach będą się pojawiać na Twitterze na kanale printpaperorg.
Wesołego drukowania!
2009-10-03
helper.h
Sporo piszę w C (kod poniżej dotyczy C99). W międzyczasie urodziłem kilka jednolinijkowych makr mających przede wszystkim zwiększyć komfort pisania i przejrzystość kodu; służą głównie do alokacji pamięci.
Najczęściej używam NEW(type, n) dzięki któremu zamiast pisać:
wystarczy, że napiszę:
Są też NEW_CLEAR(t, n) i RESIZE(t, n), które odpowiadają calloc(..) i realloc(..). Niby nic wielkiego, ale zauważyłem, że często tego używam i dużo przyjemniej się pracuje nad kodem - właśnie dlatego postanowiłem się tym podzielić.
Do ściągnięcia: helper.tar.gz
Nie jestem przekonany co do tej nazwy - jest zbyt ogólna; ale to jest takie małe i ma w sobie prawie tylko proste, pomocnicze makra, że tylko taka mi pasowała.
Dodatkowo przy wstępnym portowaniu kodu z języków, w których tablice mają indeksy startujące od 1, a nie od 0 zrobiłem NEW1 i FREE1; Np:
Potem oczywiście lepiej zmodyfikować kod do normalnych indeksów, ale przydają się po to żeby na początku uniknąć błędów przy modyfikacjach w locie i skupić się na innych możliwych błędach.
Niedawno rozszerzyłem to o makra alokujące tablice, które pamiętają ilość elementów. Dzięki temu można napisać:
To jest świeże i kod dotyczący tych tablic nie jest jeszcze dobrze przetestowany - czekam na bugi i pomysły. Zależało mi na tym, żeby to było możliwie proste, oszczędne i żeby tablice te bez problemu działały z normalnymi funkcjami, które oczekują na wejściu wskaźnika. Główną motywacją było to, że nie znoszę przechowywać wielkości tablicy w osobnej zmiennej i jeszcze pamiętać jak się ona nazywa (oraz pisać funkcji, które obok argumentu typu wskaźnikowego mają kolejny argument z ilością elementów). Ogólnie można by całą ideę rozszerzyć i zapisywać trochę więcej informacji oprócz ilości elementów i zaimplementować niezły Vector, ale to wydaje się już za bardzo kłócić z zakładaną prostotą.
Działa to tak: ARRAY_NEW zwraca normalny wskaźnik zadeklarowanego typu. Dokładniej ARRAY_NEW alokuje normalną tablicę+miejsce na jeden size_t, zapisuje długość tablicy w size_t na początku i zwraca wskaźnik już do samych danych (czyli +sizeof(size_t)). Trzeba tylko pamiętać, żeby taką tablicę zwolnić za pomocą ARRAY_FREE.
Pętlę for taką jak kawałek wyżej można uprościć za pomocą makra FOREACH:
Wśród makr jest też ARRAY_RESIZE(array, new_size). Służy do zmiany rozmiaru tablicy alokowanej przy pomocy ARRAY_NEW. Przykład:
Ponieważ są to normalne wskaźniki można ich normalnie używać jak każdych innych i napisać np:
lub trochę krócej:
Dziś dodałem też makro ACCESS (ARRAY_ACCESS), które kłóci się trochę z resztą kodu, bo psuje czytelność, ale ma taką zaletę, że wplata assert(..) (z assert.h) dbający o dostęp tylko w ramach zaalokowanej pamięci na elementy.
Żeby tego użyć zamiast:
Należy napisać:
Jako jedno z zastosowań widzę implementacje złożonych algorytmów, które trudno się debuguje.
Należy jednak pamiętać, że ACCESS to makro i w tym wypadku nie powinno się podawać do niego argumentów, które mają efekty uboczne, bo zostaną wywołane kilka razy.
ACCESS korzysta z assert z assert.h, dlatego jeżeli doda się na początku programu #define NDEBUG (przed pierwszym #include lub jako parametr dla kompilatora), to w miejscach użycia ACCESS wygeneruje się kod tak samo szybki jak przy normalnych odwołaniach do tablicy.
Makra RESIZE i ARRAY_RESIZE w obecnej formie korzystają z niestandardowego operatora __typeof__ i mogą nie działać pod czymś innym niż gcc. Reszta kodu wymaga tylko malloc, calloc, realloc (oraz memset dla makr ZERO i ARRAY_ZERO, które służą do czyszczenia tablic).
Najczęściej używam NEW(type, n) dzięki któremu zamiast pisać:
SomeType *table = (SomeType *)malloc(sizeof(SomeType)*1000);
wystarczy, że napiszę:
SomeType *table = NEW(SomeType, 1000);
Są też NEW_CLEAR(t, n) i RESIZE(t, n), które odpowiadają calloc(..) i realloc(..). Niby nic wielkiego, ale zauważyłem, że często tego używam i dużo przyjemniej się pracuje nad kodem - właśnie dlatego postanowiłem się tym podzielić.
Do ściągnięcia: helper.tar.gz
Nie jestem przekonany co do tej nazwy - jest zbyt ogólna; ale to jest takie małe i ma w sobie prawie tylko proste, pomocnicze makra, że tylko taka mi pasowała.
Dodatkowo przy wstępnym portowaniu kodu z języków, w których tablice mają indeksy startujące od 1, a nie od 0 zrobiłem NEW1 i FREE1; Np:
double *vec = NEW1(double, 2000);
vec[1] = 12323; // pierwszy element
vec[2000] = 33; // ostatni element
// vec[0] = 232; - błąd
vec[1] = 12323; // pierwszy element
vec[2000] = 33; // ostatni element
// vec[0] = 232; - błąd
Potem oczywiście lepiej zmodyfikować kod do normalnych indeksów, ale przydają się po to żeby na początku uniknąć błędów przy modyfikacjach w locie i skupić się na innych możliwych błędach.
Niedawno rozszerzyłem to o makra alokujące tablice, które pamiętają ilość elementów. Dzięki temu można napisać:
double *v = ARRAY_NEW(double, 1000);
//...
for (int i = 0; i < LEN(v); ++i) {
v[i] = sin(2.*M_PI*i/LEN(v));
}
//...
ARRAY_FREE(v);
//...
for (int i = 0; i < LEN(v); ++i) {
v[i] = sin(2.*M_PI*i/LEN(v));
}
//...
ARRAY_FREE(v);
To jest świeże i kod dotyczący tych tablic nie jest jeszcze dobrze przetestowany - czekam na bugi i pomysły. Zależało mi na tym, żeby to było możliwie proste, oszczędne i żeby tablice te bez problemu działały z normalnymi funkcjami, które oczekują na wejściu wskaźnika. Główną motywacją było to, że nie znoszę przechowywać wielkości tablicy w osobnej zmiennej i jeszcze pamiętać jak się ona nazywa (oraz pisać funkcji, które obok argumentu typu wskaźnikowego mają kolejny argument z ilością elementów). Ogólnie można by całą ideę rozszerzyć i zapisywać trochę więcej informacji oprócz ilości elementów i zaimplementować niezły Vector, ale to wydaje się już za bardzo kłócić z zakładaną prostotą.
Działa to tak: ARRAY_NEW zwraca normalny wskaźnik zadeklarowanego typu. Dokładniej ARRAY_NEW alokuje normalną tablicę+miejsce na jeden size_t, zapisuje długość tablicy w size_t na początku i zwraca wskaźnik już do samych danych (czyli +sizeof(size_t)). Trzeba tylko pamiętać, żeby taką tablicę zwolnić za pomocą ARRAY_FREE.
Pętlę for taką jak kawałek wyżej można uprościć za pomocą makra FOREACH:
FOREACH(i, v) {
v[i] = sin(2.*M_PI*i/LEN(v));
}
v[i] = sin(2.*M_PI*i/LEN(v));
}
Wśród makr jest też ARRAY_RESIZE(array, new_size). Służy do zmiany rozmiaru tablicy alokowanej przy pomocy ARRAY_NEW. Przykład:
double *tab = ARRAY_NEW(double, 100);
//...
printf("%d\n", LEN(tab)); // drukuje "100"
ARRAY_RESIZE(tab, 1024);
tab[1023] = 123; // ok
printf("%d\n", LEN(tab)); // drukuje "1024"
//...
printf("%d\n", LEN(tab)); // drukuje "100"
ARRAY_RESIZE(tab, 1024);
tab[1023] = 123; // ok
printf("%d\n", LEN(tab)); // drukuje "1024"
Ponieważ są to normalne wskaźniki można ich normalnie używać jak każdych innych i napisać np:
write(fd, tab, sizeof(*tab)*LEN(tab));
lub trochę krócej:
write(fd, tab, SIZE(tab));
Dziś dodałem też makro ACCESS (ARRAY_ACCESS), które kłóci się trochę z resztą kodu, bo psuje czytelność, ale ma taką zaletę, że wplata assert(..) (z assert.h) dbający o dostęp tylko w ramach zaalokowanej pamięci na elementy.
Żeby tego użyć zamiast:
tab[i] = tab[i - 2] + tab[i - 1];
Należy napisać:
ACCESS(tab, i) = ACCESS(tab, i - 2) + ACCESS(tab, i - 1);
Jako jedno z zastosowań widzę implementacje złożonych algorytmów, które trudno się debuguje.
Należy jednak pamiętać, że ACCESS to makro i w tym wypadku nie powinno się podawać do niego argumentów, które mają efekty uboczne, bo zostaną wywołane kilka razy.
ACCESS korzysta z assert z assert.h, dlatego jeżeli doda się na początku programu #define NDEBUG (przed pierwszym #include lub jako parametr dla kompilatora), to w miejscach użycia ACCESS wygeneruje się kod tak samo szybki jak przy normalnych odwołaniach do tablicy.
Makra RESIZE i ARRAY_RESIZE w obecnej formie korzystają z niestandardowego operatora __typeof__ i mogą nie działać pod czymś innym niż gcc. Reszta kodu wymaga tylko malloc, calloc, realloc (oraz memset dla makr ZERO i ARRAY_ZERO, które służą do czyszczenia tablic).
2009-09-06
Czwarty wymiar piątaka
Przypadkiem to wyszło; właściwie nic takiego, ale wizualnie po prostu fajne. Efekt występuje szczególnie w aparatach z sensorem CMOS - tutaj widok z mojej kamery internetowej. Chciałem wychwycić monetę tak jakby stała w miejscu, żeby lepiej zaobserwować co się z nią dzieje kiedy szybko się obraca - przy normalnych ustawieniach była rozmazana tak jak dla oka, tyle że tym razem ustawiłem bardzo krótki czas naświetlania (rzędu tysięcznych sekundy; tak jakbym używał stroboskopu) i wtedy pojawiło wygięcie. O tym efekcie/iluzji było jakiś czas temu w internecie głośno szczególnie przez wykonaną iphonem fotografię obracającego się śmigła. W tanich chińskich aparatach (szczególnie w tanich komórkach) ten sam efekt występuje w normalnych warunkach już przy bardzo powolnym ruchu. Podobny efekt uzyskuje się też przez skanowanie obracających się obiektów (kto nie wkładał głowy albo ręki do skanera ten nie zrozumie ;) ):
Te konkretne skany są ze stycznia 2003. Było tego dużo więcej - efekt burzy mózgów ze Snują. Wtedy jeszcze tego nie znaliśmy - wyszło naturalnie, ale chyba każdy kto trochę dłużej eksperymentuje ze skanerem lub czymś podobnym, w końcu na to natrafia.
Nie wiem kto to odkrył jako pierwszy, ale już w 1988 opanował do perfekcji i wykorzystał w swoim filmie "The Fourth Dimension" Zbigniew Rybczyński. To, że o tym dziś piszę zupełnie przypadkowo (naprawdę) zbiega się z tym, że wczoraj Zbig we własnej osobie wygłosił wykład w toruńskim CSW z okazji festiwalu animacji Pole Widzenia 2. Poniżej jego szkice dotyczące tego filmu:
Efekt przesunięcia fragmentów obrazu w czasie bardzo fajnie popchnęli dalej ludzie pracujący nad deformowalnym ekranem dotykowym. Fizyczna instalacja nazywa się Khronos Projector (2005). Dostępna jest też uproszczona wersja on-line w postaci apletów Java. Tak to wygląda w realu:
Tutaj pełniejsza prezentacja:
Pewno można by to dużo dalej posunąć z dzisiejszymi ekranami dotykowymi i szybkimi kartami graficznymi. Ciekawie by się grało w FPS-a, w którym gra odbywałaby się w różnym czasie w hipotetycznej przestrzeni, w której światło porusza się bardzo wolno i w zależności od np. natężenia pola grawitacyjnego można by lokalnie zaglądać w przeszłość (i nawet modyfikować jej niektóre aspekty)... jak znam życie pewno już coś podobnego jest; nie wiem, nie gram za często ;)
2009-08-16
Blip #komentarze z poziomu bloga
Od kilku dobrych tygodniu rozwijam pewną ideę - głównie w chwilach kiedy mam dosyć tego co muszę robić, ale jednak coś bym porobił. Jest to pewna usługa webowa... ale nie będę na razie tłumaczył o co dokładnie chodzi, sory. (zbyt wiele zmiennych dotyczących pomysłu jeszcze muszę ustalić, a nie chciałbym być źle zrozumiany)
Piszę o tym, bo jako jeden z pierwszych testów wykorzystałem to do umożliwienia umieszczania komentarzy do moich blipnięć (czyli mini postów po lewej) z poziomu tego bloga. Wcześniej blipnięcia można było komentować tylko będąc użytkownikiem samego Blipa - i tylko wtedy można je było oglądać.
Żeby przejść na stronę komentarzy odpowiedniego posta wystarczy kliknąć mały dymek obok daty publikacji.
Jeżeli jakiś post już został skomentowany to dymek przy nim ma kolor zielony.
Ta aplikacja do komentarzy to oddzielna strona - niezależna od Bloggera czy Blipa - i nie ma z nimi nic wspólnego. To co udostępniam to tylko jedno z bardzo wielu zastosowań mojej aplikacji.
Pod komentarzami można się podpisywać, ale nie trzeba. Nie ma żadnej rejestracji - chociaż pewną namiastkę chyba wprowadzę. Taką małą, dosyć fajną cechą jest to, że komentarze pojawiają się na żywo - jak w chacie. Tzn. jeżeli ktoś wstawi komentarz to inni ludzie, którzy oglądają aktualnie tę samą stronę komentarzy zobaczą nowy komentarz prawie od razu - bez przeładowania strony.
Na razie jest to wersja bardzo robocza - nic jeszcze nie jest ustalone - wszystko może się zmienić (łącznie z nazwą); ale uwagi, pomysły i informacje o wszelkich wykrytych błędach w działaniu będą dla mnie bardzo cenne (np w komentarzach do tego posta).
Kilka informacji o komentarzach, bardzo w skrócie:
Linki zamieniają się automatycznie na linki, nie ma żadnego formatowania. Dodatkowo można robić linki wewnętrzne przez haszowanie tekstu (Np. #tekst albo #[tekst ze spacją]) oraz komentarze do komentarzy przez skróty do urli komentarzy (np. #c:11). Strony komentarzy skrótów pojawiających się w komentarzach dostają komentarz z informacją wsteczną gdzie się pojawiły.
Na razie mam pod to bardzo tani i raczej słaby hosting - osobiście mi wystarcza, ale dla hashbasha zaraz będzie za mały. Już planuję przerzutkę na coś poważniejszego. Dodatkowo ponieważ jest to bardzo preview to prosiłbym o nie siłowanie się z serwerem, bo nie zaimplementowałem jeszcze kilku ważnych zabezpieczeń przed wandalami (np. można floodować bez ograniczeń; szczególnie jeżeli idzie o (niekompletne) api; work bardzo in progress).
Piszę o tym, bo jako jeden z pierwszych testów wykorzystałem to do umożliwienia umieszczania komentarzy do moich blipnięć (czyli mini postów po lewej) z poziomu tego bloga. Wcześniej blipnięcia można było komentować tylko będąc użytkownikiem samego Blipa - i tylko wtedy można je było oglądać.
Żeby przejść na stronę komentarzy odpowiedniego posta wystarczy kliknąć mały dymek obok daty publikacji.
Jeżeli jakiś post już został skomentowany to dymek przy nim ma kolor zielony.
Ta aplikacja do komentarzy to oddzielna strona - niezależna od Bloggera czy Blipa - i nie ma z nimi nic wspólnego. To co udostępniam to tylko jedno z bardzo wielu zastosowań mojej aplikacji.
Pod komentarzami można się podpisywać, ale nie trzeba. Nie ma żadnej rejestracji - chociaż pewną namiastkę chyba wprowadzę. Taką małą, dosyć fajną cechą jest to, że komentarze pojawiają się na żywo - jak w chacie. Tzn. jeżeli ktoś wstawi komentarz to inni ludzie, którzy oglądają aktualnie tę samą stronę komentarzy zobaczą nowy komentarz prawie od razu - bez przeładowania strony.
Na razie jest to wersja bardzo robocza - nic jeszcze nie jest ustalone - wszystko może się zmienić (łącznie z nazwą); ale uwagi, pomysły i informacje o wszelkich wykrytych błędach w działaniu będą dla mnie bardzo cenne (np w komentarzach do tego posta).
Kilka informacji o komentarzach, bardzo w skrócie:
Linki zamieniają się automatycznie na linki, nie ma żadnego formatowania. Dodatkowo można robić linki wewnętrzne przez haszowanie tekstu (Np. #tekst albo #[tekst ze spacją]) oraz komentarze do komentarzy przez skróty do urli komentarzy (np. #c:11). Strony komentarzy skrótów pojawiających się w komentarzach dostają komentarz z informacją wsteczną gdzie się pojawiły.
Na razie mam pod to bardzo tani i raczej słaby hosting - osobiście mi wystarcza, ale dla hashbasha zaraz będzie za mały. Już planuję przerzutkę na coś poważniejszego. Dodatkowo ponieważ jest to bardzo preview to prosiłbym o nie siłowanie się z serwerem, bo nie zaimplementowałem jeszcze kilku ważnych zabezpieczeń przed wandalami (np. można floodować bez ograniczeń; szczególnie jeżeli idzie o (niekompletne) api; work bardzo in progress).
2009-07-29
Wsteczna propagacja błędu - sinusy na wejściu
W lutym pisałem o zabawie z implementacją zadania na sieci neuronowe. Od tego czasu praktycznie nic przy tym nie grzebałem. Dziś (spontanicznie) postanowiłem zobaczyć co by się stało gdybym na wejścia dorzucił jeszcze sinusy zależne od x i osobne zależne od y. Generalnie pomysł użycia sinusów wziąłem od Fouriera; ale nie myślałem, że będzie taki fajny rezultat. W porównaniu z sieciami bez tych wejść zbieganie do obrazu źródłowego jest bardzo szybkie (przynajmniej dla małych sieci) - dla sieci o dwóch warstwach ukrytych po 32 neurony, po kilkuset tysiącach przykładowych pikseli obraz jest już rozpoznawalny. A dodaję tylko 16 neuronów na wejściu, czyli jak mam 32 w pierwszej ukrytej to przybywa 16*32 wag, czyli tragedii nie ma.
Niestety nie mam teraz czasu na sklejenie filmu z większą ilością przykładów, ale skleiłem poniższe porównanie dwóch konfiguracji nowej sieci z sinusami i sieci starej bez sinusów - pokazuję iterację nr 2mln i 10mln (tzn 10mln losowych przykładów i wstecznych propagacji).
Obraz jeszcze trochę zbiega w kolejnych milionach iteracjach, ale dobrze by było zmniejszać stałą uczenia, bo zachowuje się trochę zbyt chaotycznie. Niestety na razie brak czasu na eksperymenty.
Dorzuciłem jeszcze wejście z szumem (dokładnie rozkład jednostajny 0.1-0.9). Chyba Keyer mi podrzucił ten pomysł jakiś miesiąc temu - nie jestem pewien czy dobrze to zapamiętałem - miało to chyba poprawić zbieganie. Dalej chyba Gruby to chyba sugerował. Być może wszystko pokręciłem i źle to robię. Na oko nie zauważyłem dużej zmiany (suma błędów na pikselach po 5-10mln iteracji) - może trochę szybciej, nie wiem na ile się zasugerowałem... Ale nawet po dłuższym uczeniu widzę, że wagi przy tym wejściu utrzymuja się na poziomie do +/-0.2 - na razie nie mam zdania co z tym zrobić.
A co do metody to muszę jeszcze poeksperymentować z różnymi funkcjami na wejściu, zmienna stała uczenia i przeróżnych konfiguracjach. Jestem też ciekaw jak obraz generowany przez sieć wyglada w powiększeniu oraz jak algorytm poradzi sobie przy bardziej jednolitych teksturach (skały, kamienie, drewno, czy jakieś powierzchnie organiczne). Może by to miało jakieś zastosowanie przy generowaniu tekstur proceduralnych? Ciekawe też jak bardzo można przyciać dokładność wag by nie zakłócić za bardzo całego obrazu.
---
UPDATE (2009-11-10 01:27): Filmik z porównaniem procesów uczenia się zwykłych sieci (pierwszy wiersz) z sieciami z sinusami na wejściu takimi jak wyżej (w drugim wierszu):
Jakiś czas temu Maja podpowiedziała mi, że dorzucanie tak zmodyfikowanych wejść przypomina tzw. metody kernelowe - czyli ma to jakąś nazwę. Ogólnie wizualnie ciekawe wyniki dostawałem też przy innych funkcjach okresowych (np. sign(sin(.)) czy arcsin(sin(.)) ), ale przy zdjęciach uczenie najlepiej przebiega przy sinusach... Nie próbowałem jeszcze całej bazy razem z cosinusami.
Zrezygnowałem na razie z wejścia losowego.
Niestety nie mam teraz czasu na sklejenie filmu z większą ilością przykładów, ale skleiłem poniższe porównanie dwóch konfiguracji nowej sieci z sinusami i sieci starej bez sinusów - pokazuję iterację nr 2mln i 10mln (tzn 10mln losowych przykładów i wstecznych propagacji).
Obraz jeszcze trochę zbiega w kolejnych milionach iteracjach, ale dobrze by było zmniejszać stałą uczenia, bo zachowuje się trochę zbyt chaotycznie. Niestety na razie brak czasu na eksperymenty.
Dorzuciłem jeszcze wejście z szumem (dokładnie rozkład jednostajny 0.1-0.9). Chyba Keyer mi podrzucił ten pomysł jakiś miesiąc temu - nie jestem pewien czy dobrze to zapamiętałem - miało to chyba poprawić zbieganie. Dalej chyba Gruby to chyba sugerował. Być może wszystko pokręciłem i źle to robię. Na oko nie zauważyłem dużej zmiany (suma błędów na pikselach po 5-10mln iteracji) - może trochę szybciej, nie wiem na ile się zasugerowałem... Ale nawet po dłuższym uczeniu widzę, że wagi przy tym wejściu utrzymuja się na poziomie do +/-0.2 - na razie nie mam zdania co z tym zrobić.
A co do metody to muszę jeszcze poeksperymentować z różnymi funkcjami na wejściu, zmienna stała uczenia i przeróżnych konfiguracjach. Jestem też ciekaw jak obraz generowany przez sieć wyglada w powiększeniu oraz jak algorytm poradzi sobie przy bardziej jednolitych teksturach (skały, kamienie, drewno, czy jakieś powierzchnie organiczne). Może by to miało jakieś zastosowanie przy generowaniu tekstur proceduralnych? Ciekawe też jak bardzo można przyciać dokładność wag by nie zakłócić za bardzo całego obrazu.
---
UPDATE (2009-11-10 01:27): Filmik z porównaniem procesów uczenia się zwykłych sieci (pierwszy wiersz) z sieciami z sinusami na wejściu takimi jak wyżej (w drugim wierszu):
Jakiś czas temu Maja podpowiedziała mi, że dorzucanie tak zmodyfikowanych wejść przypomina tzw. metody kernelowe - czyli ma to jakąś nazwę. Ogólnie wizualnie ciekawe wyniki dostawałem też przy innych funkcjach okresowych (np. sign(sin(.)) czy arcsin(sin(.)) ), ale przy zdjęciach uczenie najlepiej przebiega przy sinusach... Nie próbowałem jeszcze całej bazy razem z cosinusami.
Zrezygnowałem na razie z wejścia losowego.
Labels: neuron
2009-07-14
Pan Magoo
Na specjalne zamówienie tym razem coś bez pierwiastków i wykresów.
Niedawno znalazłem świetnego bloga obracającego się wokół stylu kreskówek z lat pięćdziesiątych: Cartoon Modern - niestety już nie rozwijany, ale archiwa bardzo bogate (jest też kontynuacja w postaci mini bloga na tumblr). Można tam znaleźć szczególnie dużo informacji o studiu animacji UPA (filmografia na Wikipedii) - to oni są odpowiedzialni z krótki film Rooty Toot Toot - nomanacja do Oskara 1952 (ostatnio na blipie przyczepiłem do tego link). I to właśnie oni wymyślili Pana Magoo. A teraz do rzeczy:
O tym jak Magoo latał (Oskar 1955)
Pan Magoo w centrum lotów kosmicznych
Pan Magoo na siłowni
A na deser:
Magoo w reklamie piwa
Niedawno znalazłem świetnego bloga obracającego się wokół stylu kreskówek z lat pięćdziesiątych: Cartoon Modern - niestety już nie rozwijany, ale archiwa bardzo bogate (jest też kontynuacja w postaci mini bloga na tumblr). Można tam znaleźć szczególnie dużo informacji o studiu animacji UPA (filmografia na Wikipedii) - to oni są odpowiedzialni z krótki film Rooty Toot Toot - nomanacja do Oskara 1952 (ostatnio na blipie przyczepiłem do tego link). I to właśnie oni wymyślili Pana Magoo. A teraz do rzeczy:
O tym jak Magoo latał (Oskar 1955)
Pan Magoo w centrum lotów kosmicznych
Pan Magoo na siłowni
A na deser:
Magoo w reklamie piwa
2009-06-16
1.f/sqrtf(float): FPU vs Q3 hack vs SSE
Na wczorajszym seminarium (piątkowym w poniedziałek) pojawiło się pytanie: o ile szybszy jest szybki odwrotny pierwiastek kwadratowy z Quake 3 od normalnego 1.f/sqrtf(float)? Generalnie hack jest bardzo szczegółowo opisany, przedyskutowany i świetnie zbadany - warto przejrzeć papiery w przypisach do artykułu na wikipedii (szczególnie 4 i 5).
Ale nie wiedziałem jakie są konkretne liczby, dlatego napisałem prosty test mający to zmierzyć. Oprócz porównania ze zwykłym liczeniem na FPU postanowiłem to jeszcze porównać z taką samą szybką metodą, ale z dodatkowym przybliżeniem oraz z liczeniem na SSE (instrukcja rsqrtps - Packed Single-Precision Floating-Point Square Root Reciprocal). Wszystkie te metody są bardzo przybliżone - na FPU jest najdokładniej; mimo iż liczby są na floatach to na FPU wszystko jest liczone na rejestrach 80-bitowych. Za to na SSE wszystko jest liczone na liczbach 32-bitowych - gdzieś czytałem, że jest to realizowane podobną metodą co właśnie szybka odwrotny pierwiastek kwadratowy. Konkretne rezultaty poniżej.
Kilka uwag:
Jeżeli chodzi o metodę z Quake 3 to biorę minimalnie poprawioną wersję z magiczną stałą 0x5f375a86 wykorzystując wynik z Fast Inverse Square Root (pdf) - w tekście dużo detali oraz odpowiednia stała dla wersji na double. (W Understanding and improving fast inverse square root hack uzyskano prawie taką samą stałą: 0x5f375a84).
Niżej metodę na FPU nazwałem "Normal", a metodę szybką "Fast". Oprócz SSE jest jeszcze metoda szybka z dodatkową iteracją przybliżenia pod koniec, która nazywam "FastExtra" - na podstawie: Magical Square Root Implementation In Quake III .
Uwaga: flaga -O2 w gcc psuje kod szybkiej metody, dlatego poniższe testy kompilowałem z -O1.
TEST: 1.f/sqrtf(float)
Test wygląda tak: Biorę 1048576 pseudolosowych liczb z zakresu (0.0, 1000.0) (zawsze te same), które zapisuję w tablicy i 200 razy wykonuję obliczenie odwrotności pierwiastka kwadratowego na każdej z kolejnych liczb - czyli wykonuję 209715200 takich operacji. Wynik każdej operacji zapisuję w innej tablicy niż tablica wejściowa. Testuję kolejno każdą z metod. Czas jest liczony w sekundach. Błąd jest liczony w porównaniu z metodą na FPU.
Wyniki na laptopie z Pentium 4-M 2GHz, DDR-266:
Wyniki na desktopie z AMD Athlon X2 BE-2300, DDR2-800:
Wersja SSE jest wyraźnie najszybsza, ale największy błąd wyniósł nawet 3e-3. Warto zwrócić uwagę na to jak różny jest średni błąd na różnych maszynach, ale to chyba normalne przy SSE - będzie to wyraźniejsze w drugim teście.
TEST: Normalize(float3)
Ponieważ odwrotność pierwiastka wykorzystuję przede wszystkim do normalizowania trójwymiarowych wektorów, postanowiłem jeszcze przetestować jak każda z tych metod się do tego nadaje. Tutaj obliczając błąd sprawdzam jak daleka jest długość wyjściowych wektorów od pożądanej jedynki.
Wyniki na Pentium 4-M 2GHz, DDR-266:
Wyniki na AMD Athlon X2 BE-2300, z DDR2-800:
Przyspieszenie w przypadku SSE jest zadowalające - na szybszej maszynie wręcz bardzo. Na laptopie (Pentium 4M) wersja Fast jest całkiem blisko SSE, chociaż znowu ma wyraźnie gorszą dokładność. Jednak tam zyski z SSE najprawdopodobniej giną przesłonięte przez opóźnienia wywołane fatalnie wolną pamięcią.
Jak to zwykle przy SSE, ważne jest odpowiednie ułożenie danych tak by dało się obliczać cztery wektory na raz, co w niektórych warunkach może być kłopotliwe, jeżeli nie niemożliwe. Dla wszystkich metod ułożyłem dane wektor po wektorze - tzn. xyzxyzxy..., za to dla SSE dane są wcześniej przekształcone w ciąg: xxxxyyyyzzzzxxx... . We wszystkich przypadkach pamięć jest wyrównana do granicy 16-bajtowej. Jeżeli idzie o implementację dla SSE, to wygląda następująco:
Mięso w pętli jest tłumaczone przez gcc na:
Na koniec: 1.f/sqrtf(0.f)
Przy okazji normalizacji interesującym było jeszcze obejrzeć jak się te metody zachowują dla wektorów zerowych. FPU i SSE jako wynik 1/sqrt(0) zwracają inf (FPU bardzo powoli - można to jakoś przełączyć?), co dalej przemnożone przez 0.f daje nan. Za to w przypadku metod Fast i FastExtra dostałem bardzo duże wartości (Fast: 1.98e+19, FastExtra: 2.97e+19), które śmiało można przemnożyć przez 0. Wynik będzie taki, że normalizacja wektora na FPU i SSE daje wektor [NaN, NaN, NaN], za to pozostałe przybliżone metody dadzą [0.f, 0.f, 0.f]. Zależy od zastosowania, który wynik jest bardziej pożądany, chociaż pewno w większości przypadków lepiej jest standardowo dorzucić ifa po obliczeniu sumy kwadratów współrzędnych, ale może to skomplikować obliczenia na SSE.
Błądy ujemny i dodatnie:
Sprawdziłem jeszcze (na Pentium 4M), dla normalizacji wektora, jak dokładnie wyglądają najgorsze błędy dodatnie i ujemne:
Jak widać błąd ujemny jest minimalny - nawet mniejszy niż przy FPU. Błąd ujemny wskazuje na to o ile długość znormalizowanego wektora przekracza jedynkę. Czyli w przypadku metod Fast i FastExtra większość błędów powoduje wyprodukowanie wektorów trochę krótszych niż wektory normalne. Przy tak dużym błędzie absolutnym jest to fantastyczna własność - przynajmniej w moich zastosowaniach.
Ale nie wiedziałem jakie są konkretne liczby, dlatego napisałem prosty test mający to zmierzyć. Oprócz porównania ze zwykłym liczeniem na FPU postanowiłem to jeszcze porównać z taką samą szybką metodą, ale z dodatkowym przybliżeniem oraz z liczeniem na SSE (instrukcja rsqrtps - Packed Single-Precision Floating-Point Square Root Reciprocal). Wszystkie te metody są bardzo przybliżone - na FPU jest najdokładniej; mimo iż liczby są na floatach to na FPU wszystko jest liczone na rejestrach 80-bitowych. Za to na SSE wszystko jest liczone na liczbach 32-bitowych - gdzieś czytałem, że jest to realizowane podobną metodą co właśnie szybka odwrotny pierwiastek kwadratowy. Konkretne rezultaty poniżej.
Kilka uwag:
Jeżeli chodzi o metodę z Quake 3 to biorę minimalnie poprawioną wersję z magiczną stałą 0x5f375a86 wykorzystując wynik z Fast Inverse Square Root (pdf) - w tekście dużo detali oraz odpowiednia stała dla wersji na double. (W Understanding and improving fast inverse square root hack uzyskano prawie taką samą stałą: 0x5f375a84).
Niżej metodę na FPU nazwałem "Normal", a metodę szybką "Fast". Oprócz SSE jest jeszcze metoda szybka z dodatkową iteracją przybliżenia pod koniec, która nazywam "FastExtra" - na podstawie: Magical Square Root Implementation In Quake III .
Uwaga: flaga -O2 w gcc psuje kod szybkiej metody, dlatego poniższe testy kompilowałem z -O1.
TEST: 1.f/sqrtf(float)
Test wygląda tak: Biorę 1048576 pseudolosowych liczb z zakresu (0.0, 1000.0) (zawsze te same), które zapisuję w tablicy i 200 razy wykonuję obliczenie odwrotności pierwiastka kwadratowego na każdej z kolejnych liczb - czyli wykonuję 209715200 takich operacji. Wynik każdej operacji zapisuję w innej tablicy niż tablica wejściowa. Testuję kolejno każdą z metod. Czas jest liczony w sekundach. Błąd jest liczony w porównaniu z metodą na FPU.
Wyniki na laptopie z Pentium 4-M 2GHz, DDR-266:
Input: 200x1048576 positive floats (max: 999.998291)
Normal time : 9.333797
Fast time : 2.626831 (speedup: 3.553254)
FastExtra time: 4.312687 (speedup: 2.164265)
SSE time : 1.755530 (speedup: 5.316797)
NORMAL/FAST:
mean error: -0.0000602813
error variance: 0.0000000182
max error: -0.0563874133
NORMAL/FASTEXTRA:
mean error: -0.0000001185
error variance: 0.0000000000
max error: -0.0001132605
NORMAL/SSE:
mean error: 0.0000000009
error variance: 0.0000000002
max error: -0.0032330328
Wyniki na desktopie z AMD Athlon X2 BE-2300, DDR2-800:
Input: 200x1048576 positive floats (max: 999.998291)
Normal time : 6.598592
Fast time : 1.681536 (speedup: 3.924146)
FastExtra time: 2.198866 (speedup: 3.000907)
SSE time : 0.680382 (speedup: 9.698363)
NORMAL/FAST:
mean error: -0.0000602813
error variance: 0.0000000182
max error: -0.0563874133
NORMAL/FASTEXTRA:
mean error: -0.0000001185
error variance: 0.0000000000
max error: -0.0001132605
NORMAL/SSE:
mean error: 0.0000014102
error variance: 0.0000000001
max error: -0.0032330328
Wersja SSE jest wyraźnie najszybsza, ale największy błąd wyniósł nawet 3e-3. Warto zwrócić uwagę na to jak różny jest średni błąd na różnych maszynach, ale to chyba normalne przy SSE - będzie to wyraźniejsze w drugim teście.
TEST: Normalize(float3)
Ponieważ odwrotność pierwiastka wykorzystuję przede wszystkim do normalizowania trójwymiarowych wektorów, postanowiłem jeszcze przetestować jak każda z tych metod się do tego nadaje. Tutaj obliczając błąd sprawdzam jak daleka jest długość wyjściowych wektorów od pożądanej jedynki.
Wyniki na Pentium 4-M 2GHz, DDR-266:
Input: 200x1048576 3d non zero vectors
Normal time : 10.012473
Fast time : 6.510379 (speedup: 1.537925)
FastExtra time: 7.751822 (speedup: 1.291628)
SSE time : 5.167708 (speedup: 1.937508)
NORMAL:
mean error: 0.0000000000
error variance: 0.0000000000
max error: -0.0000000503
FAST:
mean error: 0.0009001758
error variance: 0.0000003486
max error: 0.0017512076
FASTEXTRA:
mean error: 0.0000017375
error variance: 0.0000000000
max error: 0.0000046443
SSE:
mean error: -0.0000003470
error variance: 0.0000000133
max error: -0.0003259508
Wyniki na AMD Athlon X2 BE-2300, z DDR2-800:
Input: 200x1048576 3d non zero vectors
Normal time : 7.214732
Fast time : 4.266968 (speedup: 1.690833)
FastExtra time: 5.219265 (speedup: 1.382327)
SSE time : 2.117634 (speedup: 3.406978)
NORMAL:
mean error: 0.0000000000
error variance: 0.0000000000
max error: -0.0000000503
FAST:
mean error: 0.0009001758
error variance: 0.0000003486
max error: 0.0017512076
FASTEXTRA:
mean error: 0.0000017375
error variance: 0.0000000000
max error: 0.0000046443
SSE:
mean error: -0.0000209675
error variance: 0.0000000057
max error: -0.0002579081
Przyspieszenie w przypadku SSE jest zadowalające - na szybszej maszynie wręcz bardzo. Na laptopie (Pentium 4M) wersja Fast jest całkiem blisko SSE, chociaż znowu ma wyraźnie gorszą dokładność. Jednak tam zyski z SSE najprawdopodobniej giną przesłonięte przez opóźnienia wywołane fatalnie wolną pamięcią.
Jak to zwykle przy SSE, ważne jest odpowiednie ułożenie danych tak by dało się obliczać cztery wektory na raz, co w niektórych warunkach może być kłopotliwe, jeżeli nie niemożliwe. Dla wszystkich metod ułożyłem dane wektor po wektorze - tzn. xyzxyzxy..., za to dla SSE dane są wcześniej przekształcone w ciąg: xxxxyyyyzzzzxxx... . We wszystkich przypadkach pamięć jest wyrównana do granicy 16-bajtowej. Jeżeli idzie o implementację dla SSE, to wygląda następująco:
// SIZE - ilość wektorów
float xyz4[3*SIZE]; // wektory wejściowe (xxxxyyyyzzzzxxx...) (16-byte aligned)
float xyzs[3*SIZE]; // wektory wyjściowe - znormalizowane (16-byte aligned)
...
for (int i = 0; i < 3*SIZE; i += 12) {
__m128 r0 = _mm_load_ps(xyz4 + i);
__m128 r1 = _mm_load_ps(xyz4 + i + 4);
__m128 r2 = _mm_load_ps(xyz4 + i + 8);
__m128 r3 = r0;
__m128 r4 = r1;
__m128 r5 = r2;
r0 = _mm_mul_ps(r0, r0);
r1 = _mm_mul_ps(r1, r1);
r2 = _mm_mul_ps(r2, r2);
r0 = _mm_add_ps(r0, r1);
r0 = _mm_add_ps(r0, r2);
r0 = _mm_rsqrt_ps(r0);
r3 = _mm_mul_ps(r3, r0);
r4 = _mm_mul_ps(r4, r0);
r5 = _mm_mul_ps(r5, r0);
_mm_store_ps(xyzs + i, r3);
_mm_store_ps(xyzs + i + 4, r4);
_mm_store_ps(xyzs + i + 8, r5);
}
Mięso w pętli jest tłumaczone przez gcc na:
movl -24(%ebp), %eax
movaps (%eax,%ebx), %xmm3
leal 16(%ebx), %edx
movaps (%eax,%edx), %xmm4
leal 32(%ebx), %ecx
movaps (%eax,%ecx), %xmm5
movaps %xmm3, %xmm0
mulps %xmm3, %xmm0
movaps %xmm4, %xmm1
mulps %xmm4, %xmm1
movaps %xmm5, %xmm2
mulps %xmm5, %xmm2
addps %xmm1, %xmm0
addps %xmm2, %xmm0
rsqrtps %xmm0, %xmm0
mulps %xmm0, %xmm3
mulps %xmm0, %xmm4
mulps %xmm0, %xmm5
movl -40(%ebp), %eax
movaps %xmm3, (%eax,%ebx)
movl -40(%ebp), %eax
movaps %xmm4, (%eax,%edx)
movl -40(%ebp), %eax
movaps %xmm5, (%eax,%ecx)
Na koniec: 1.f/sqrtf(0.f)
Przy okazji normalizacji interesującym było jeszcze obejrzeć jak się te metody zachowują dla wektorów zerowych. FPU i SSE jako wynik 1/sqrt(0) zwracają inf (FPU bardzo powoli - można to jakoś przełączyć?), co dalej przemnożone przez 0.f daje nan. Za to w przypadku metod Fast i FastExtra dostałem bardzo duże wartości (Fast: 1.98e+19, FastExtra: 2.97e+19), które śmiało można przemnożyć przez 0. Wynik będzie taki, że normalizacja wektora na FPU i SSE daje wektor [NaN, NaN, NaN], za to pozostałe przybliżone metody dadzą [0.f, 0.f, 0.f]. Zależy od zastosowania, który wynik jest bardziej pożądany, chociaż pewno w większości przypadków lepiej jest standardowo dorzucić ifa po obliczeniu sumy kwadratów współrzędnych, ale może to skomplikować obliczenia na SSE.
Błądy ujemny i dodatnie:
Sprawdziłem jeszcze (na Pentium 4M), dla normalizacji wektora, jak dokładnie wyglądają najgorsze błędy dodatnie i ujemne:
NORMAL:
max +error: 0.0000000502
max -error: -0.0000000503
FAST:
max +error: 0.0017512076
max -error: -0.0000000384
FASTEXTRA:
max +error: 0.0000046443
max -error: -0.0000000468
SSE:
max +error: 0.0003250295
max -error: -0.0003259508
Jak widać błąd ujemny jest minimalny - nawet mniejszy niż przy FPU. Błąd ujemny wskazuje na to o ile długość znormalizowanego wektora przekracza jedynkę. Czyli w przypadku metod Fast i FastExtra większość błędów powoduje wyprodukowanie wektorów trochę krótszych niż wektory normalne. Przy tak dużym błędzie absolutnym jest to fantastyczna własność - przynajmniej w moich zastosowaniach.
Starsze posty dostępne w archiwum
Ładuję...
