;****************************************************************************
;;; MICRO CONNECTIONIST QUANTIZER
;;; (C)1990, Peter Desain & Henkjan Honing
;;; in Common Lisp (uses loop macro)
;;; see https://www.mcg.uva.nl/abstracts/dh-91-b.html

;;; utilities

(define-modify-macro multf (factor) *)
(define-modify-macro divf (factor) /)
(define-modify-macro zerof () (lambda(x) (declare (ignore x)) 0))

(defun print-state (time intervals)
  "Print elements of interval vector"
  (loop initially (format t "~%~2D: " time)
        for index below (length intervals) do
        (format t "~2,1,5$ " (aref intervals index))))

(defmacro with-adjacent-intervals 
 (vector (a-begin a-end a-sum b-begin b-end b-sum) &body body)
  "Setup environment for each interaction of (sum-)intervals"
  `(loop with length = (length ,vector)
         for ,a-begin below (1- length)
         do (loop for ,a-end from ,a-begin below (1- length)
                  sum (aref ,vector ,a-end) into ,a-sum 
                  do (loop with ,b-begin = (1+ ,a-end)
                           for ,b-end from ,b-begin below length
                           sum (aref ,vector ,b-end) into ,b-sum
                           do ,@body))))

;;; interaction function

(defun delta (a b minimum peak decay)
  "Return change for two time intervals"
  (let* ((inverted? (<= a b))
         (ratio (if inverted? (/ b a)(/ a b)))
         (delta-ratio (interaction ratio peak decay))
         (proportion (/ delta-ratio (+ 1 ratio delta-ratio))))
    (* minimum (if inverted? (- proportion) proportion))))
  
(defun interaction (ratio peak decay)
  "Return change of time interval ratio"
  (* (- (round ratio) ratio)
     (expt (abs (* 2 (- ratio (floor ratio) 0.5))) peak)
     (expt (round ratio) decay)))

;;; quantization procedures

(defun quantize! 
 (intervals &key (iterations 20) (peak 5) (decay -1))
  "Quantize data of inter-onset intervals"
  (let* ((length (length intervals))
         (changes (make-array length :initial-element 0.0))
         (minimum (loop for index below length 
                        minimize (aref intervals index))))
    (loop for count to iterations 
          do (print-state count intervals)
          (update! intervals minimum changes peak decay))))

(defun update! (intervals minimum changes peak decay)
  "Update all intervals synchronously"
  (with-adjacent-intervals intervals 
    (a-begin a-end a-sum b-begin b-end b-sum)
    (let ((delta (delta a-sum b-sum minimum peak decay)))
      (propagate! changes a-begin a-end (/ delta a-sum))
      (propagate! changes b-begin b-end (- (/ delta b-sum)))))
  (enforce! changes intervals))

(defun propagate! (changes begin end change)
  "Derive changes of basic-intervals from sum-interval change"
  (loop for index from begin to end 
        do (incf (aref changes index) change)))

(defun enforce! (changes intervals)
  "Effectuate changes to intervals"
  (loop for index below (length intervals)
        do (multf (aref intervals index) 
                  (1+ (aref changes index)))
           (zerof (aref changes index))))

#|
;;; examples
(quantize! (vector 1.1 2.0 2.9))

(quantize! (vector 1.177 0.592 0.288 0.337 0.436 0.337 0.387 0.600 
                  0.634 0.296 0.280 0.296 0.346 1.193) :iterations 10)

|#