view utils/harmony/gamma.lisp @ 61:c911d65ae94d

pointy-clicky misspelling correction

darcs-hash:20070627102839-dc3a5-16cb5156dafee4b4a1aa3442274fab6c7260ed9c.gz

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