Optymalizowanie prostego programu próbnika gibbs Common Lisp

Question

Jako ćwiczenie przepisałem przykładowy program w poście na blogu Gibbs sampler in various languages (revisited) autorstwa Darrena Wilkinsona.

Kod pojawi się poniżej. Ten kod działa na moim (5-letni) maszynie w około 53 sekund, przy użyciu SBCL 1.0.56, tworząc rdzeń obrazu za pomocą buildapp, a następnie uruchomić go z

time ./gibbs > gibbs.dat 

Ponieważ było to, w jaki sposób obliczono czasy dla inne języki w poście, myślałem, że zrobię coś porównywalnego Kod C w poście działa w około 25 sekund. Chciałbym spróbować przyspieszyć kod Lisp jeśli to możliwe.

############################## 
gibbs.lisp 
############################## 
(eval-when (:compile-toplevel :load-toplevel :execute) 
    (require :cl-rmath) (setf *read-default-float-format* 'double-float)) 

(defun gibbs (N thin) 
    (declare (fixnum N thin)) 
    (declare (optimize (speed 3) (safety 1))) 
    (let ((x 0.0) (y 0.0)) 
    (declare (double-float x y)) 
    (print "Iter x y") 
    (dotimes (i N) 
     (dotimes (j thin) 
    (declare (fixnum i j)) 
    (setf x (cl-rmath::rgamma 3.0 (/ 1.0 (+ (* y y) 4)))) 
    (setf y (cl-rmath::rnorm (/ 1.0 (+ x 1.0)) (/ 1.0 (sqrt (+ (* 2 x) 2)))))) 
     (format t "~a ~a ~a~%" i x y)))) 

(defun main (argv) 
    (declare (ignore argv)) 
    (gibbs 50000 1000)) 

Potem zbudował wykonywalny gibbs z wywołaniem sh gibbs.sh z gibbs.sh jak

################## 
gibbs.sh 
################## 
buildapp --output gibbs --asdf-tree /usr/share/common-lisp/source/ --asdf-tree /usr/local/share/common-lisp/source/ --load-system cl-rmath --load gibbs.lisp --entry main 

mam 6 kompilatora notatek podczas kompilacji z SBCL 1.0.56, które odtwarzają poniżej. Nie jestem pewien, co z nimi zrobić, ale byłbym wdzięczny za wszelkie wskazówki.

; compiling file "/home/faheem/lisp/gibbs.lisp" (written 30 MAY 2012 02:00:55 PM): 

; file: /home/faheem/lisp/gibbs.lisp 
; in: DEFUN GIBBS 
;  (SQRT (+ (* 2 X) 2)) 
; 
; note: unable to 
; optimize 
; due to type uncertainty: 
; The result is a (VALUES (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT)) 
;       &OPTIONAL), not a (VALUES FLOAT &REST T). 

;  (/ 1.0d0 (SQRT (+ (* 2 X) 2))) 
; 
; note: unable to 
; optimize 
; due to type uncertainty: 
; The second argument is a (OR (DOUBLE-FLOAT 0.0) 
;        (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX 
;                DOUBLE-FLOAT). 
; 
; note: forced to do static-fun Two-arg-/ (cost 53) 
;  unable to do inline float arithmetic (cost 12) because: 
;  The second argument is a (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT. 
;  The result is a (VALUES (OR (COMPLEX DOUBLE-FLOAT) (DOUBLE-FLOAT 0.0)) 
;        &OPTIONAL), not a (VALUES DOUBLE-FLOAT &REST T). 

;  (CL-RMATH:RGAMMA 3.0d0 (/ 1.0d0 (+ (* Y Y) 4))) 
; 
; note: doing float to pointer coercion (cost 13) 

;  (SQRT (+ (* 2 X) 2)) 
; 
; note: doing float to pointer coercion (cost 13) 

;  (CL-RMATH:RNORM (/ 1.0d0 (+ X 1.0d0)) (/ 1.0d0 (SQRT (+ (* 2 X) 2)))) 
; 
; note: doing float to pointer coercion (cost 13) 
; 
; compilation unit finished 
; printed 6 notes 

; /home/faheem/lisp/gibbs.fasl written 
; compilation finished in 0:00:00.073 

UPDATE 1: Rainer Joswig's answer zwrócił uwagę, że argument SQRT może być ujemna, co był źródłem niejasnych notatek kompilatora widziałem, mianowicie

The result is a (VALUES (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT)) 
    ;       &OPTIONAL), not a (VALUES FLOAT &REST T). 

Kompilator narzekał, że ponieważ nie wiedział, czy wartość argumentu była dodatnia, wynik może być liczbą zespoloną. Ponieważ w tym przykładzie wartość x jest próbką różniącą się od rozkładu gamma, jest zawsze większa niż 0. Została wskazana przez Stephana na liście dyskusyjnej użytkowników SBCL, (patrz druga wiadomość w wątku "Optimizing a simple Common Lisp Gibbs sampler program", że ten może być rozwiązany poprzez uznanie X będzie większa lub zerowy następująco,

(declare (type (double-float 0.0 *) x)) 

Zobacz Common Lisp Hyperspec do odpowiedniej dokumentacji o FLOAT types i Interval Designators.

Thi Wydaje się, że trochę przyspiesza kod. Obecnie jest niezawodnie poniżej 52 sekund, ale nadal niewiele zyskuje. Pozostawia również notatki o

note: doing float to pointer coercion (cost 13) 

Jeśli nie jest to naprawić jakiegoś powodu, chciałbym wiedzieć, dlaczego. Poza tym wyjaśnienie znaczenia notatki byłoby interesujące. Co konkretnie oznacza słowo pointer? Czy jest to związane z faktem wywoływania funkcji C? Również koszt 13 nie wydaje się bardzo przydatny. Co mierzy się?

Rainer zasugerował także, że może istnieć możliwość zmniejszenia zużycia, co może zmniejszyć czas uruchamiania. Nie wiem, czy redukcja redukcji jest możliwa, czy zmniejszyłaby czas wykonywania, , ale byłbym zainteresowany opiniami i podejściami. Ogólnie rzecz biorąc, wydaje się, że niewiele można zrobić, aby poprawić wydajność tej funkcji. Może jest za mały i prosty.

Answer 1

Należy pamiętać, że Common Lisp ma operatora specjalnego THE. Pozwala zadeklarować typy dla wyników wyrażenia. To na przykład pozwala zawęzić typy, jeśli to możliwe.

Na przykład, jaki jest wynik (SQRT somefloat)? Może to być wartość zmiennoprzecinkowa, ale może to być liczba zespolona, ​​jeśli somefloat jest ujemna. Jeśli wiesz, że jakiś pływak jest zawsze dodatni (i tylko wtedy), możesz napisać: (the double-float (sqrt somefloat)). Kompilator może wtedy być w stanie wygenerować bardziej wydajny kod.

Należy również pamiętać, że Common Lisp ma deklaracje OPTIMIZE. Jeśli potrzebujesz najszybszego kodu, upewnij się, że odpowiednio je ustawiłeś. Możliwe tylko dla poszczególnych funkcji. Zwykle lepiej niż globalnie zmieniać optymalizację, aby była bardzo agresywna.

Common Lisp ma funkcję DISASSEMBLE, która pozwala spojrzeć na skompilowany kod.

Następnie jest makro TIME. Interesujące informacje, które z niego otrzymasz, obejmują to, ile to robi. Z arytmetyką podwójnego pływania jest prawdopodobnie duża suma. Przydałaby się pomoc na liście mailingowej SBCL. Może ktoś może ci powiedzieć, jak tego uniknąć.

Answer 2

Działa to dla mnie:

(sqrt (the (double-float 0d0) (+ (* 2d0 x) 2d0)))