;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; whole model is short
; whole model is presented with paper
; in an open source language
; stored in a version control system that anyone can access
; comes with functions called eg1 eg2, etc
; random number generator's seed can be reset
; comes with a decision maker that return minimal descretized theories
;   and those theories come with a "brittleness estimate"
; no globals (there will be multiple what-ifs)

; OODA loop: observe orient decide act
; lots of support for observe, little for orient decide act

;; Kolb's model therefore works on two levels - a four-stage cycle:

;;    1. Concrete Experience - (CE)
;;    2. Reflective Observation - (RO)
;;    3. Abstract Conceptualization - (AC)
;;    4. Active Experimentation - (AE)

;; and a four-type definition of learning styles, (each representing the combination of two preferred styles, rather like a two-by-two matrix of the four-stage cycle styles, as illustrated below), for which Kolb used the terms:

;;    1. Diverging (CE/RO)
;;    2. Assimilating (AC/RO)
;;    3. Converging (AC/AE)
;;    4. Accommodating (CE/AE)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;unit tests From Peter Seiblel's excellent text: http://gigamonkeys.com/book

(defparameter *test-name* nil)
(defparameter *tests* nil)

(defmacro deftest (name parameters &body body)
  `(progn
     (push ',name *tests*) 
     (defun ,name ,parameters
       (let ((*test-name* (append *test-name* (list ',name))))
         ,@body))))

(defmacro check (&body forms)
  `(combine-results
    ,@(loop for f in forms collect `(report-result ,f ',f))))

(defmacro with-gensyms ((&rest names) &body body)
  `(let ,(loop for n in names collect
              `(,n (make-symbol ,(string n))))
     ,@body))

(defmacro combine-results (&body forms)
  (with-gensyms (result)
    `(let ((,result t))
       ,@(loop for f in forms collect
              `(unless ,f (setf ,result nil)))
       ,result)))

(let ((passes 0) (fails 0))

  (defun tests-reset()
    (setf passes 0)
    (setf fails 0))

  (defun tests-report ()
    (format t "~%PASSES: ~a (~a %)~%FAILS : ~a~%"
            passes (* 100 (/ passes (+ passes fails)))
            fails))

  (defun report-result (result form)
    (if result
        (and (incf passes)
             (format t "."))
        (and (incf fails)
             (format t "~%fail ... ~a: ~a~%" *test-name* form)))
    result)
  )

(defmacro dotests (&body body)
  `(progn
     (tests-reset)
     ,@body
     (tests-report)))

(defun tests ()
  (dotests
    (mapc #'funcall (reverse *tests*))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; random number generator

(defparameter *seed0* 10013)
(defparameter *seed*  *seed0*)

(defun reset-seed () (setf *seed* *seed0*))

(defun park-miller-randomizer ()
  "The Park-Miller multiplicative congruential randomizer
  (CACM, October 88, Page 1195). Creates pseudo random floating
  point numbers in the range 0.0 < x <= 1.0."
  (let ((multiplier 
	 16807.0d0);16807 is (expt 7 5)
	(modulus 
	 2147483647.0d0)) ;2147483647 is (- (expt 2 31) 1)
    (let ((temp (* multiplier *seed*)))
      (setf *seed* (mod temp modulus))
      (/ *seed* modulus))))

(defun my-random (n)
  "Returns a pseudo random floating-point number in range 0.0 <= number < n"
  (let ((random-number (park-miller-randomizer)))
    (* n (- 1.0d0 random-number))))

(defun my-randone ()
  "returns  a float 0..1"
  (let ((big most-positive-fixnum ))
    (/ (my-random big) big)))

(defun my-random-int (n)
  "Returns a pseudo-random integer in the range 0 <= n-1."
  (let ((random-number (/ (my-random 1000.0) 1000)))
    (floor (* n random-number))))

(defun within (range)
  (let ((min (first range))
        (max (rest range)))
  (+ min (* (- max min) (my-randone)))))

(defun within-int (range)
  (let ((min (first range))
        (max (rest  range)))
    (+ min (* (my-random-int (1+ (- max  min)))))))

(defun random-demo ()
  (let (counts out)
    (labels 
	((sorter (x y) (< (car x) (car y)))
         (zap    ()    (setf out nil)
	               (reset-seed) 
	               (setf counts (make-hash-table)))
	 (inc    (n)   (setf (gethash n counts) 
			     (1+  (gethash n counts 0)))) 
	 (cache  (k v) (push (list k v) out)))
      (zap)
      (dotimes (i 10000) 
        (inc  (my-random-int 5)))
      (maphash #'cache counts)
      (sort out #'sorter))))

(deftest test-random ()
  (check
    (equalp (random-demo) (random-demo))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; sample from a normal distribution

(defun normal (m s) 
  (let ((deltam (* s (box_muller))))
    (if (>= (my-randone) 0.5)
        (+ m deltam)
        (- m deltam))))

(defun box_muller ()
  (multiple-value-bind (w x) (w)
    (* x  (sqrt (/ (* -2.0 (log w)) w)))))

(defun w ()
  (let* ((x1 (* 2.0 (my-randone)))
         (x2 (* 2.0 (my-randone)))
         (w  (+ (* x1 x1) (* x2 x2))))
    (if (> w 1)
        (w)
        (values w x1))))

(defun sum (l &key (what #'identity))
  "Sum a list of numbers."
  (let ((n 0))
    (dolist (x l n)
      (incf n (funcall what x)))))  

(defun mean (l &key (value #'identity))
  "Return the mean and sum of list of numbers."
  (let ((sum 0)
	(n   0))
    (dolist (one l (/ sum n))
      (incf n)
      (incf sum (funcall value one)))))

(defun stdev (n sum sumSq) 
  "Compute the mean and standard deviation."
  (sqrt (/ (- sumSq(/ (* sum sum) n)) (- n 1))))

(defun list2stats (l &key (value #'identity))
  "Return the mean and standard deviation of a list of numbers."
  (let ((n 0) (sum 0) (sumSq 0))
    (dolist (x l (values (/ sum n ) (stdev n sum sumSq)))
      (incf n)
      (let ((one (funcall value x)))
        (incf sum one)
        (incf sumSq (* one one))))))

(deftest test-normal ()
  (let ((mean0  10)
        (stdev0 2)
        (epsilon 0.2)
        out)
    (reset-seed)
    (dotimes (i 500) 
      (push (normal mean0 stdev0) out))
    (multiple-value-bind (mean1 stdev1)
        (list2stats out)
      (check
        (and (<= (- mean0  epsilon) mean1  (+ mean0  epsilon))
             (<= (- stdev0 epsilon) stdev1 (+ stdev0 epsilon)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun poisson (lamba)
  (let ((l (exp (* -1 lamba)))
        (k 0)
        (p 1))
    (poisson1 l k p)))

(defun poisson1 (l k p)
  (if (< p l)
      (- k 1)
      (poisson1 l (1+ k) (* p (my-randone)))))

(deftest test-poisson ()
  (let ((mean0  1)
        (epsilon 0.05)
        out)
    (reset-seed)
    (dotimes (i 5000) 
      (push (poisson mean0) out))
    (check (<= (- mean0  epsilon) (mean out)  (+ mean0  epsilon)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defstruct opt
  (value       '(30 . 500)) ; value of requirements
  (cost         '(1 . 100)) ; cost of requirements
  (versions     '(2 .   6)) ; how many versions
  (implemented     0.8)        ; requirements mplemented at each version
  (requirements '(1 .  25)) ; number of requimrents
  (b0-%-mean       40)
  (b0-%-sigma       5)
  (lamba             2) ; rate of new requirement discovery
  (sigma             10)) ; variance in requirement value
)

(defstruct requirement
  value 
  cost )

(defun new-requirement (o)
  (make-requirement :value (within (opt-value o))
                    :cost  (within (opt-cost  o))))

(defstruct wme
  versions 
  b0  
  requirements 
  heap
  plan
  done
  revalue-after)

(defun wme-new (o)
  (let* ((w            (make-wme))
         (versions     (within-int (opt-versions o)))
         (requirements (within-int (opt-requirements o)))
         (b0           (round (* requirements
                                 (/ (normal (opt-b0-%-mean o) 
                                            (opt-b0-%-sigma o))
                                    100)))) 
         (heap        (todo requirements o))
         plan)
    (dotimes (i b0)
      (push (pop heap) plan))
    (setf (wme-heap w)          heap
          (wme-b0   w)          b0
          (wme-todo w)          plan
          (wme-versions w)      versions
          (wme=requirements w)  requirements
          (wme-revalue-after w) (wme-revalue-after0 w))
    p))

(defun wme-stop (w) 
  (> (my-randone) (/ 1 (expt (wme-versions w) 0.33))))
    
(defun wme-revalue-after0 (w)
  (let* ((b0 (wme-b0 w))
         (base-requirements (subseq (wme-heap w) 0 (1- b0)))
         ($base-requirements
          (sum base-requirements :what #'requirement-cost)))
    (/ $base-requirements (wme-versions w))))

;;;;;;;;;;;;;
(defun todo (n o)
  (let (out)
    (dotimes (i n out)
      (push (new-requirement o) out)))) 

(defun run (&optional (sorter #'sort-on-value) (pruner #'indentity))
  (let* ((o (make-opt))
         (w (wme-new o)))
    (run1 w o sorter pruner)

(defun step (w o > /)
 (if (wme-plan w)
     (unless (wme-stop w)
       (sort-and-prune           w > /)
       (adjust-some-values       w)
       (retire-some-requirements w)
       (add-new-requirements     w)
       (step                     w o > /))))

(defun sort-and-prine (w > /)
       (funcall / (funcall > w)))

(defun retire-some-requirements (w o)
  (let ((n      (length (wme-plan w)))
        (enough (round (* n (opt-implemented o)))))
    (dotimes (i enough)
      (push (pop (wme-plan w)) (wme-done w)))))

(defun adjust-some-requirements (w)
  (labels ((move () (push (pop (wme-plan w) (wme-done ))))
           (stay () t))
    (let ((enough (wme-revalue-after w)))
      (adjust-some-requirements1 w enough #'wme-plan #'wme-heap #'move #'stay))))

(defun adjust-some-requirements1 (w enough plan heap move stay)
  (if (<= enough 0)
      (adjust-some-requirements1 w enough stay stay)
)

(defun add-some-requirements (w o)
  (dotimes (i (poisson (opt-lamba o)))
    (push (new-requirement o) (wme-plan w))))

(defun sort-on-value (w)
  (setf (wme-heap w) 
        (sort (wme-heap w) #'> :key #'wme-value)))