#| #################################################################### 
kent v0.1 : the metropilis simulator annealler + back select
copyright (c) 2007, Tim Menzies (tim@menzies.us)

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation;  version 2.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA.
#################################################################### |#

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;; lib
(defmacro o (x) 
  `(progn 
    (format t "[~a]=[~a]~%" (quote ,x) ,x) ,x))

(defmacro oo (&rest l)
  `(progn 
    ,@(mapcar #'(lambda(x) `(o ,x)) l)))

(defmacro show (&rest tests)
  `(dolist (one (reverse ',tests) t)
    (print  one) (princ " => " ) (princ (eval one)) ))  

(defun say (n) 
  (let ((r 0.001))
    (* (round (/ n r))
       r)))

(defun no-op1 (x) x)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;; maths utils

(defun between (min max)
  (+ min (* (random 1.0) (- max min))))

(defun gaussian (m s)
  (+ m (+ s (box-muller))))

(defun box-muller ()
  "returns a number 0..1"
  (let (x1 x2  w)
    (loop
      (setf x1 (-  (random 2.0) 1)
            x2 (-  (random 2.0) 1)
            w  (+  (* x1 x1) (* x2 x2)))
      (if (< w 1)
        (return (* x1 (sqrt (/ (* -2.0 (log w)) w))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;;  ranges

(defun nominal (range &key  nudge  value (validp nil))  
  (cond (validp                    (and (member value range) value))
        ((and nudge (zerop nudge)) value)
        (t                         (elt range  (random (length range))))))

(defun ration-int (min max &key  nudge value (validp nil))
  (ration  min  max 
           :nudge nudge :value value :validp validp :filter 'round))

(defun ration (min max &key  nudge value validp  (filter 'no-op1))
 (cond ((< max min) (ration  max  min 
                             :nudge nudge
                             :value value :nudge nudge 
                             :validp validp :filter filter))
         (validp     (and (<= min value max) value))
         (nudge      (let
                       ((delta (* (- max min) nudge 0.5)))
                       (setf min   (- value delta)
                             max   (+ value delta))
                       (funcall filter (between min max))))
        (t           (funcall filter (between min max)))))

(defun normal-int (m s &key  nudge value validp (filter 'round))
  (normal m s :nudge nudge :value value 
          :validp validp :filter filter))

(defun normal (m s &key nudge value validp  (filter 'no-op1))
  (let ((min (- m (* 4 s)))
        (max (+ m (* 4 s))))
    (cond (validp  (and (<= min value max) value))
          (nudge   (ration min max :nudge nudge))
          (t       (funcall filter (gaussian m s))))))

(defun ordinal (range &key  nudge value (validp nil))
  (let ((dim (- (length range) 1)))
    (cond (validp (and (member value range) value))
          (nudge  (let*
                    ((pos0  (position value range))
                     (pos1  (ration-int 0 dim 
                                        :value pos0 :nudge nudge)))
                    (elt range pos1)))
          (t     (let
                   ((pos1 (ration-int  0  dim)))
                   (elt range pos1))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;; things

(defstruct thing name functor args) 

(defun question-mark-reader (stream char)
  (declare (ignore char))
  (let ((thing (read stream t nil)))
    `(access ',thing cache things)))

(set-macro-character #\? #'question-mark-reader)

(defmacro def (name in defs of functor has &rest args)
  (declare (ignore in of has))
  `(setf  ,defs (def1 ',name ',functor ',args ,defs)))

(defun def1 (name functor args defs) 
  (let* ((new (make-thing :name name  :functor functor :args args))
         (cell (cons name new)))
    (cond ((assoc name defs)  (rplacd (assoc name defs) new))
          (defs               (setf defs (cons cell defs)))
          (t                  (setf defs (list cell))))
    defs))

(defun access1 (x)
  (print 3)
  (print x)
  (let ((f (thing-functor x))
        (a (thing-args    x)))
     (apply f a)))

(defmacro access (key)
  `(mutiple-value-setq (value cache things)
                       (access1 ,key cache things)))

(defun access1 (key cache things)
  (let (value
         (tmp (assoc key cache)))
    (cond (tmp  (setf value (cdr tmp)))
          (t    (let ((maker  (cdr (assoc  key things))))
                  (setf value  (access1 maker)
                        cache  (cons (cons key  value) cache)))))
    (values value cache things)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;;  simulated annealer
(defun metropolis (&key (steps 10000)
                        (cost 'cost)
                        (next 'mutate)
                        (first 'start))
  (let* (best
          cache
          (time 1) 
          (bestS (expt 10 32))
          new
          newS
          (current (funcall first cache))
          (currentS (funcall cost cache)))
    (loop
      (if (<= time 0 ) 
          (return (values best bestS)))
      (setf new   (funcall next current)
            newS  (funcall cost new)) 
      (if (< newS bestS)
          (progn  (format t " !~a ~a~%"  (say time) (say (- bestS newS)))
                  (setf best   (copy-list new)
                        bestS  newS)))
      (if (accept (- newS currentS) time)
          (setf current  (copy-list new)
                currentS newS))
      (setf time 
            (- time (/ 1.0 steps)))
      )))

(defun accept(d n)
  (let ((e 2.781281828))
  (cond ((< d 0)                                    t)
        ((< (random 1.0) (expt e (* -1.0 (/ d n)))) t)
        (t                                          nil))))

(defun mutate (wme)
  (mapcar #'(lambda (part)
              (let* ((old (cdr part))
                     (new (+ old (between -0.5 0.5))))
                (cons (car part) new)))
          wme))

(defun defs () 
 (let (cache things)
   (def a in things of ration     has 0 12)
   (def b in things of ordinal    has 1 2 3)
   (def c in things of ration-int has 0 12)
   (def d in things of normal     has 20 2)
   (def e in things of nominal    has name age shoesize)
   (def f in things of normal     has 200 20) 
   (format t "[~a] [~a]~%" ?a ?a)
   t
   ))
  
;(defmodel cost
;  (abs (+ (* 3 ?a ) (* -1 ?b) 2)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;;  demos

(defun ranges ()
  (show   (dotimes (i 50) (format t " ~a"  (normal-int  10 2)))
          (dotimes (i 50) (format t " ~a"  (normal  10 2)))
          (dotimes (i 10) 
            (princ  (ordinal '(1 2 3  4 5 6 7 8 9 10) 
                             :value 6 :nudge 1)))
          (dotimes (i 10) 
            (princ  (ordinal '(1 2 3  4 5 6 7 8 9 10) 
                             :value 6 :nudge 0.01)))
          (ordinal '(1 2 3) :value 3 :validp t)
          (ordinal '(1 2 3) :value 9 :validp t)
          (dotimes (i 10) 
            (princ  (ordinal '(1 2 3))))
          (ration-int 0 12 :value 6 :nudge 1)
          (ration-int 0 12 :value 6 :nudge 0.01)
          (ration-int 0 12 :value 18 :validp t)
          (ration 0 12 :value 6 :nudge 1)
          (ration 0 12 :value 6 :nudge 0.01)
          (ration 0 12 :value 18 :validp t)
          (ration 0 12 :value 8 :validp t)
          (ration 0 12)
          (nominal  '(a b c))
          (nominal  '(a b c) :nudge 1 :value 'b)
          (nominal  '(a b c) :validp t :value 'b)
          (nominal  '(a b c) :validp t :value 'z)
          ))