diff utils/harmony/gamma.lisp @ 33:d1010755f507

Large upload of local changes. Many additions, such as harmony and piece-level objects

darcs-hash:20070413100909-f76cc-a8aa8dfc07f438dc0c1a7c45cee7ace2ecc1e6a5.gz

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/utils/harmony/gamma.lisp	Fri Apr 13 11:09:09 2007 +0100
@@ -0,0 +1,233 @@
+;; CSR's gamma functions. I don't really know what's happening here,
+;; but it seems to work... Called from within harmonic_analysis.lisp
+
+(in-package #:amuse-harmony)
+
+(let ((p (make-array 
+          9 :element-type 'double-float
+          :initial-contents ; cf. Lanczos' Approximation on Wikipedia.
+          '(0.99999999999980993d0 676.5203681218851d0 -1259.1392167224028d0 
+            771.32342877765313d0 -176.61502916214059d0 12.507343278686905d0
+            -0.13857109526572012d0 9.9843695780195716d-6 
+            1.5056327351493116d-7)))
+      (c (make-array 8 :element-type 'double-float
+                     :initial-contents 
+                     (mapcar (lambda (x) (float x 1d0))
+                             '(1/12 -1/360 1/1260 -1/1680 1/1188 -691/360360
+                               1/156 -3617/122400)))))
+  (defun gamma/posreal (x)
+    (declare (type (double-float (0d0)) x))
+    (locally 
+        (declare (optimize speed))
+      (labels ((corr (x)
+                 (declare (type double-float x))
+                 (let ((y (/ 1.0 (* x x))))
+                   (do* ((i 7 (1- i))
+                         (r (aref c i) (+ (* r y) (aref c i))))
+                        ((<= i 0) (the (double-float (0d0)) (exp (/ r x))))
+                     (declare (type double-float r)))))
+               (lngamma/lanczos (x)
+                 (declare (type (double-float 0.5d0) x))
+                 (let ((x (1- x)))
+                   (do* ((i 0 (1+ i))
+                         (ag (aref p 0) (+ ag (/ (aref p i) (+ x i)))))
+                        ((>= i 8)
+                         (let ((term1 (* (+ x 0.5d0) (log (/ (+ x 7.5d0) #.(exp 1d0)))))
+                               (term2 (+ #.(log (sqrt (* 2 pi))) (log ag))))
+                           (+ term1 (- term2 7.0d0))))
+                     (declare (type (double-float (0d0)) ag)))))
+               (gamma/xgthalf (x)
+                 (declare (type (double-float 0.5d0) x))
+                 (cond
+                   ((= x 0.5d0) 1.77245385090551602729817d0)
+                   ((< x 5d0) (exp (lngamma/lanczos x)))
+                   ;; the GNU scientific library suggests a third branch
+                   ;; for x < 10, but in fact for our purposes the
+                   ;; errors are under control.
+                   (t
+                    (let* ((p (expt x (* x 0.5)))
+                           (e (exp (- x)))
+                           (q (* (* p e) p))
+                           (pre (* #.(sqrt 2) #.(sqrt pi) q (/ (sqrt x)))))
+                      (* pre (corr x)))))))
+        (if (> x 0.5d0)
+            (gamma/xgthalf x)
+            (let ((g (gamma/xgthalf (- 1d0 x))))
+              (declare (type double-float g))
+              (/ pi g (sin (* pi x)))))))))
+              
+(defun beta (alpha)
+  (let ((sum 0)
+        (product 1))
+    (dotimes (i (length alpha) (/ product (gamma/posreal (float sum 1d0))))
+      (let ((alpha_i (aref alpha i)))
+        (setq product (* product (gamma/posreal (float alpha_i 1d0))))
+        (incf sum alpha_i)))))
+
+;; this is p(d|c_i), where the chord model is represented as a
+;; Dirichlet distribution with parameters alpha_i
+#+nil
+(defun likelihood (observations alpha)
+  ;; observations is a 12d vector summing to 1, alpha summing to 3
+  (assert (< (abs (1- (loop for d across observations sum d)))
+             (* single-float-epsilon 12)))
+  (let ((pairs (loop for alpha_i across alpha
+                     for d_i across observations
+                     if (zerop alpha_i)
+                     do (unless (zerop d_i)
+                          (return-from likelihood 0.0d0))
+                     else 
+                     collect (cons alpha_i d_i))))
+    ;; now all the CARs of pairs have strictly positive alpha_i
+    ;; values, and the d_i are the corresponding observations.
+    (let ((alpha_prime (map 'simple-vector #'car pairs)))
+      (* (/ (beta alpha_prime))
+         (reduce #'* pairs :key (lambda (pair)
+                                  (expt (cdr pair) (1- (car pair)))))))))
+
+
+(defun likelihood (observations alpha &optional (power 1))
+  (assert (< (abs (1- (loop for d across observations sum d)))
+             (* single-float-epsilon 12)))
+  (let* ((alpha_0 (loop for a across alpha sum a))
+	 (alpha_0+1 (1+ alpha_0))
+	 (alpha_00+1 (* alpha_0+1 power))
+	 (alpha_00 (- alpha_00+1 1))
+	 (alpha (map 'simple-vector 
+		     (lambda (x) (* x (/ alpha_00 alpha_0))) 
+		     alpha))
+	 (quantum (expt 0.0005 power))) ; sake of argument
+    (let ((pairs
+           (loop for d_i across observations
+	      for alpha_i across alpha
+	      if (zerop alpha_i)
+	      do (unless (zerop d_i)
+		   (return-from likelihood 0.0d0))
+	      else collect (cons d_i alpha_i))))
+      (assert (= (length pairs) 12))
+      (let ((alpha_prime (map 'simple-vector #'cdr pairs)))
+        (* (/ (beta alpha_prime)) (/ quantum)
+           (reduce #'* pairs :key 
+                   (lambda (pair)
+                     (if (zerop (car pair))
+                         (* (/ (cdr pair)) (expt quantum (cdr pair)))
+                         (* quantum (expt (car pair) (1- (cdr
+							  pair))))))))))))
+
+(defun 4ple-likelihood (pitch-classes chord-probabilities intervals level sum)
+  (let ((observations))
+    (dolist (interval intervals)
+      (push (aref pitch-classes interval) observations)
+      (decf sum (aref pitch-classes interval)))
+    (push sum observations)
+    (likelihood/fourwise (make-array (length observations) :initial-contents (reverse observations))
+			 chord-probabilities level)))
+
+(let ((alpha_0s '((1 . 12) (1/2 . 9) (1/4 . 30))))
+  (defun likelihood/fourwise (observations alpha &optional (power 1))
+    (assert (< (abs (1- (loop for d across observations sum d)))
+               (* single-float-epsilon 4)))
+    (let* ((alpha_0 (loop for a across alpha sum a))
+           (alpha (map 'simple-vector
+                       (lambda (x) (* x (/ (cdr (assoc power alpha_0s))
+                                           alpha_0)))
+                       alpha))
+           (quantum (expt 0.0005 power))) ; sake of argument
+      (let ((pairs
+             (loop for d_i across observations
+                   for alpha_i across alpha
+                   if (zerop alpha_i)
+                   do (unless (zerop d_i)
+                        (return-from likelihood/fourwise 0.0d0))
+                   else collect (cons d_i alpha_i))))
+        (assert (= (length pairs) 4))
+        (let ((alpha_prime (map 'simple-vector #'cdr pairs)))
+          (* (/ (beta alpha_prime)) (/ quantum)
+             (reduce #'* pairs :key 
+                     (lambda (pair)
+                       (if (zerop (car pair))
+                           (* (/ (cdr pair)) (expt quantum (cdr pair)))
+                           (* quantum (expt (car pair) (1- (cdr pair)))))))))))))
+
+(defparameter *alpha-scale* 1)
+(defparameter *beta-scale* 1)
+(defparameter *alpha* #(2.925 1.95 1.625))
+(defparameter *beta* #(0.87 4.46))
+(defparameter *minimal-betas* '((1 . 8.75) (1/2 . 3) (1/4 . 2.5)))
+(defparameter *full-betas-1* '((1 . 14) (1/2 . 6) (1/4 . 8))) ;; gets 144|51
+(defparameter *full-betas* '((1 . 14) (3/4 . 6) (1/2 . 6) (1/4 . 8))) ;; guess...
+(defparameter *betas* *full-betas*)
+
+(defun 3ple-likelihood (pitch-classes chord-probabilities non-chord intervals level sum
+			&optional (alpha *alpha*) (beta *beta*))
+  (declare (ignore chord-probabilities non-chord))
+  (let ((observations))
+    (dolist (interval intervals)
+      (push (aref pitch-classes interval) observations)
+      (decf sum (aref pitch-classes interval)))
+    (push sum observations)
+    (when (= sum 1) ;; So there are no chord notes at all... ?? move to likelihood function
+      (return-from 3ple-likelihood 0))
+    (likelihood/threeplusonewise (make-array (length observations) :initial-contents (reverse observations))
+				 alpha beta
+				 level)
+    #+nil (likelihood/threeplusonewise (make-array (length observations) :initial-contents (reverse observations))
+				 (map 'vector #'(lambda (x) (* x *alpha-scale*))
+				      #(2.925 1.95 1.625))
+				 #+nil (map 'vector #'(lambda (x) (* x *beta-scale*))
+					    #(0.87 4.46))
+				 (make-array 2 :initial-contents (list 0.87 *beta-scale*))
+				 level)))
+
+;;  #(0.87 4.46) worked
+
+;;; observations is a non-relativistic four-vector of proportions:
+;;; #(tonic mediant dominant other).
+;;;
+;;; alpha is a three-vector of Dirichlet parameters for proportions of
+;;; tonic, mediant, dominant of the total chord notes.
+;;;
+;;; beta is a two-vector of Dirichlet parameters for proportions of
+;;; non-chord vs chord notes.
+;;;
+;;; suggested values:
+;;;   alpha: #(2.925 1.95 1.625) ; (tonic, mediant, dominant)
+;;;    beta: #(0.44 1.76) ; (non-chord, chord)
+;;;
+;;; model:
+;;;   p(tmdo|c) = p(o|c)p(tmd|oc)
+;;; where
+;;;   p(o|c) ~ Beta(0.44,1.76) (i.e. p({o,o'}|c) ~ Dir({0.44,1.76})
+;;; and
+;;;   p(tmd|oc) ~ (1-o) Dir({2.925,1.95,1.625})
+(defun likelihood/threeplusonewise 
+    (observations alpha beta &optional (power 1))
+  (assert (< (abs (1- (loop for d across observations sum d)))
+             (* single-float-epsilon 4)))
+  (let* ((quantum (expt 0.00000005 power))) ; sake of argument
+    (let* ((o (aref observations 3))
+           (pairs
+            (loop repeat 3
+                  for d_i across observations
+                  for alpha_i across alpha
+                  if (zerop alpha_i)
+                  do (unless (zerop d_i)
+                       (return-from likelihood/threeplusonewise 0.0d0))
+                  else collect (cons (/ d_i (- 1 o)) alpha_i))))
+      (assert (= (length pairs) 3))
+      (assert (< (abs (1- (loop for (d) in pairs sum d))) (*
+							   single-float-epsilon 3)))
+      (flet ((key (pair)
+               (if (zerop (car pair))
+                   (* (/ (cdr pair)) (expt quantum (cdr pair)))
+                   (* quantum (expt (car pair) (1- (cdr pair)))))))
+        (let ((alpha_prime (map 'simple-vector #'cdr pairs)))
+          (* 
+           ;; p(o|c)
+           (/ (beta beta)) (/ quantum)
+           (reduce #'* (list (cons o (aref beta 0))
+                             (cons (- 1 o) (aref beta 1)))
+                   :key #'key)
+           ;; p(tmd|oc)
+           (/ (beta alpha_prime)) (/ quantum)
+           (reduce #'* pairs :key #'key)))))))