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