home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / lisp.lzh / XLisp-Stat / maximize.lsp < prev    next >
Text File  |  1990-10-03  |  12KB  |  312 lines
  1. (provide "maximize")
  2.  
  3. ;;;;
  4. ;;;; Mode Info Prototype
  5. ;;;;
  6.  
  7. (defproto minfo-proto '(internals))
  8.  
  9. (send minfo-proto :add-method :isnew #'|minfo-isnew|)
  10. (send minfo-proto :add-method :maximize #'|minfo-maximize|)
  11. (send minfo-proto :add-method :loglaplace #'|minfo-loglap|)
  12.  
  13. (defmeth minfo-proto :x () (aref (slot-value 'internals) 3))
  14. (defmeth minfo-proto :scale () (aref (slot-value 'internals) 4))
  15. (defmeth minfo-proto :derivstep () (aref (aref (slot-value 'internals) 9) 1))
  16. (defmeth minfo-proto :tilt () (aref (aref (slot-value 'internals) 9) 6))
  17.  
  18. (defmeth minfo-proto :f (&optional (val nil set))
  19.   (when set
  20.     (send self :set-no-vals-supplied)
  21.     (setf (aref (slot-value 'internals) 0) val))
  22.   (aref (slot-value 'internals) 0))
  23.  
  24. (defmeth minfo-proto :set-no-vals-supplied ()
  25.   (setf (aref (aref (slot-value 'internals) 8) 6) 0))
  26.  
  27. (defmeth minfo-proto :exptilt (&optional (val nil set))
  28.   (if set
  29.       (let ((old (send self :exptilt)))
  30.     (setf (aref (aref (slot-value 'internals) 8) 7) (if val 1 0))
  31.     (if (and (not (or (and old val) (and (not old) (not val))))
  32.          (/= (send self :tilt) 0.0))
  33.         (send self :set-no-vals-supplied))))
  34.   (= 1 (aref (aref (slot-value 'internals) 8) 7)))
  35.  
  36. (defmeth minfo-proto :newtilt (&optional (val nil set))
  37.   (when set
  38.     (setf (aref (aref (slot-value 'internals) 9) 7) (float val))
  39.     (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
  40.   (aref (aref (slot-value 'internals) 9) 7))
  41.  
  42. (defmeth minfo-proto :gfuns (&optional (val nil set))
  43.   (when set
  44.     (if (or (not (consp val))
  45.         (not (every #'functionp val)))
  46.         (error "not all functions"))
  47.     (setf (aref (slot-value 'internals) 1) val)
  48.     (setf (aref (aref (slot-value 'internals) 8) 1) (length val))
  49.     (setf (aref (slot-value 'internals) 10) (repeat 1.0 (length val)))
  50.     (if (/= (send self :tilt) 0.0) (send self :set-no-vals-supplied)))
  51.   (aref (slot-value 'internals) 1))
  52.  
  53. (defmeth minfo-proto :cfuns (&optional (val nil set))
  54.   (when set
  55.     (if (or (not (consp val))
  56.                 (not (every #'functionp val)))
  57.             (error "not all functions"))
  58.     (setf (aref (slot-value 'internals) 2) val)
  59.     (setf (aref (aref (slot-value 'internals) 8) 2) (length val))
  60.     (setf (aref (slot-value 'internals) 7) (repeat 0.0 (length val)))
  61.     (setf (aref (slot-value 'internals) 11) (repeat 0.0 (length val)))
  62.     (send self :set-no-vals-supplied))
  63.   (aref (slot-value 'internals) 2))
  64.  
  65. (defmeth minfo-proto :ctarget (&optional (val nil set))
  66.   (when set
  67.     (if (/= (length val) (length (send self :ctarget)))
  68.         (error "bad target length"))
  69.     (setf (aref (slot-value 'internals) 7) val))
  70.   (aref (slot-value 'internals) 7))
  71.     
  72. (defmeth minfo-proto :fvals ()
  73.   (let* ((fv (aref (slot-value 'internals) 5))
  74.      (n (length (send self :x)))
  75.      (val (select fv 0))
  76.      (grad (select fv (iseq 1 n)))
  77.      (hess (matrix (list n n) (select fv (iseq (+ 1 n) (+ n (* n n)))))))
  78.     (list val grad hess)))
  79.  
  80. (defmeth minfo-proto :copy ()
  81.   (let ((obj (make-object minfo-proto))
  82.     (internals (copy-seq (slot-value 'internals))))
  83.     (dotimes (i (length internals))
  84.          (let ((x (aref internals i)))
  85.            (if (sequencep x)
  86.            (setf (aref internals i) (copy-seq x)))))
  87.     (send obj :add-slot 'internals internals)
  88.     obj))
  89.  
  90. (defmeth minfo-proto :derivscale ()
  91.   (let* ((step (^ machine-epsilon (/ 1 6)))
  92.      (hess (numhess (send self :f) (send self :x) (send self :scale) step))
  93.      (scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess)))))))
  94.     (setf hess (numhess (send self :f) (send self :x) scale step))
  95.     (setf scale (pmax (abs (send self :x)) (sqrt (abs (/ (diagonal hess))))))
  96.     (setf (aref (slot-value 'internals) 4) scale)
  97.     (setf (aref (aref (slot-value 'internals) 9) 1) step)))
  98.  
  99. (defmeth minfo-proto :verbose (&optional (val nil set))
  100.   (when set
  101.     (setf (aref (aref (slot-value 'internals) 8) 5)
  102.           (cond ((integerp val) val)
  103.             ((null val) 0)
  104.             (t 1))))
  105.   (aref (aref (slot-value 'internals) 8) 5))
  106.  
  107. (defmeth minfo-proto :backtrack (&optional (val nil set))
  108.   (if set (setf (aref (aref (slot-value 'internals) 8) 4) (if val 1 0)))
  109.   (aref (aref (slot-value 'internals) 8) 4))
  110.  
  111. (defmeth minfo-proto :maxiter (&optional (val nil set))
  112.     (if set (setf (aref (aref (slot-value 'internals) 8) 3) 
  113.           (if (integerp val) val -1)))
  114.     (aref (aref (slot-value 'internals) 8) 3))
  115.  
  116. (defmeth minfo-proto :tiltscale (&optional (val nil set))
  117.   (when set
  118.     (if (/= (length val) (length (send self :gfuns)))
  119.         (error "wrong size tilt scale sequence"))
  120.     (setf (aref (slot-value 'internals) 10) val))
  121.   (aref (slot-value 'internals) 10))
  122.  
  123. ;;;;
  124. ;;;;
  125. ;;;; Newton's Method with Backtracking
  126. ;;;;
  127. ;;;;
  128.  
  129. (defun newtonmax (f start &key 
  130.                     scale 
  131.                     (derivstep -1.0)
  132.                     (count-limit -1)
  133.                     (verbose 1)
  134.                     return-derivs)
  135. "Args:(f start &key scale derivstep (verbose 1) return-derivs)
  136. Maximizes F starting from START using Newton's method with backtracking.
  137. If RETURN-DERIVS is NIL returns location of maximum; otherwise returns
  138. list of location, unction value, gradient and hessian at maximum.
  139. SCALE should be a list of the typical magnitudes of the parameters.
  140. DERIVSTEP is used in numerical derivatives and VERBOSE controls printing
  141. of iteration information. COUNT-LIMIT limits the number of iterations"
  142.   (let ((verbose (if verbose (if (integerp verbose) verbose 1) 0))
  143.         (minfo (send minfo-proto :new f start 
  144.                      :scale scale :derivstep derivstep)))
  145.     (send minfo :maxiter count-limit)
  146.     (send minfo :derivscale)
  147.     (send minfo :maximize verbose)
  148.     (if return-derivs
  149.         (cons (send minfo :x) (- (send minfo :fvals)))
  150.         (send minfo :x))))
  151.  
  152. ;;;;
  153. ;;;;
  154. ;;;; Nelder-Mead Simplex Method
  155. ;;;;
  156. ;;;;
  157.  
  158. (defun nelmeadmax (f start &key 
  159.                      (size 1)
  160.                      (epsilon (sqrt machine-epsilon)) 
  161.                      (count-limit 500)
  162.                      (verbose t)
  163.                      (alpha 1.0) 
  164.                      (beta 0.5) 
  165.                      (gamma 2.0)
  166.                      (delta 0.5))
  167. "Args: (f start &key (size 1) (epsilon (sqrt machine-epsilon)) 
  168.           (count-limit 500) (verbose t) alpha beta gamma delta)
  169. Maximizes F using the Nelder-Mead simplex method. START can be a
  170. starting simplex - a list of N+1 points, with N=dimension of problem,
  171. or a single point. If start is a single point you should give the
  172. size of the initial simplex as SIZE, a sequence of length N. Default is
  173. all 1's. EPSILON is the convergence tolerance. ALPHA-DELTA can be used to
  174. control the behavior of simplex algorithm."
  175.     (let ((s (send simplex-proto :new f start size)))
  176.       (do ((best (send s :best-point) (send s :best-point))
  177.            (count 0 (+ count 1))
  178.            next)
  179.           ((or (< (send s :relative-range) epsilon) (>= count count-limit))
  180.            (if (and verbose (>= count count-limit))
  181.                (format t "Iteration limit exceeded.~%"))
  182.            (send s :point-location (send s :best-point)))
  183.           (setf next (send s :extrapolate-from-worst (- alpha)))
  184.           (if (send s :is-worse best next)
  185.               (setf next (send s :extrapolate-from-worst gamma))
  186.               (when (send s :is-worse next (send s :second-worst-point))
  187.                     (setf next (send s :extrapolate-from-worst beta))
  188.                     (if (send s :is-worse next (send s :worst-point))
  189.                         (send s :shrink-to-best delta))))
  190.           (if verbose 
  191.               (format t "Value = ~g~%" 
  192.                       (send s :point-value (send s :best-point)))))))
  193.           
  194.  
  195. ;;;
  196. ;;; Simplex Prototype
  197. ;;;
  198.  
  199. (defproto simplex-proto '(f simplex))
  200.  
  201. ;;;
  202. ;;; Simplex Points
  203. ;;;
  204.  
  205. (defmeth simplex-proto :make-point (x)
  206.   (let ((f (send self :f)))
  207.     (if f 
  208.         (let ((val (funcall f x)))
  209.           (cons (if (consp val) (car val) val) x))
  210.         (cons nil x))))
  211.  
  212. (defmeth simplex-proto :point-value (x) (car x))
  213.  
  214. (defmeth simplex-proto :point-location (x) (cdr x))
  215.  
  216. (defmeth simplex-proto :is-worse (x y)
  217.   (< (send self :point-value x) (send self :point-value y)))
  218.  
  219. ;;;
  220. ;;; Making New Simplices
  221. ;;;
  222.  
  223. (defmeth simplex-proto :isnew (f start &optional size)
  224.   (send self :simplex start size)
  225.   (send self :f f))
  226.  
  227. ;;;
  228. ;;; Slot Accessors and Mutators
  229. ;;;
  230.  
  231. (defmeth simplex-proto :simplex (&optional new size)
  232.   (if new
  233.       (let ((simplex 
  234.              (if (and (consp new) (sequencep (car new)))
  235.                  (if (/= (length new) (+ 1 (length (car new))))
  236.                      (error "bad simplex data")
  237.                      (copy-list new))
  238.                  (let* ((n (length new))
  239.                         (size (if size size (repeat 1 n)))
  240.                      ;   (pts (- (* 2 (uniform-rand (repeat n (+ n 1)))) 1)))
  241.                         (diag (* 2 size (- (random (repeat 2 n)) .5)))
  242.                          (pts (cons (repeat 0 n)    
  243.                                    (mapcar #'(lambda (x) (coerce x 'list))
  244.                                            (column-list (diagonal diag))))))
  245.                    (mapcar #'(lambda (x) (+ (* size x) new)) pts)))))
  246.         (setf (slot-value 'simplex) 
  247.               (mapcar #'(lambda (x) (send self :make-point x)) simplex))
  248.         (send self :sort-simplex)))
  249.   (slot-value 'simplex))
  250.  
  251. (defmeth simplex-proto :f (&optional f)
  252.   (when f
  253.         (setf (slot-value 'f) f)
  254.         (let ((simplex 
  255.                (mapcar #'(lambda (x) (send self :point-location x))
  256.                        (send self :simplex))))
  257.           (send self :simplex simplex)))
  258.   (slot-value 'f))
  259.  
  260. (defmeth simplex-proto :sort-simplex ()
  261.   (if (send self :f)
  262.       (setf (slot-value 'simplex) 
  263.             (sort (slot-value 'simplex)
  264.                   #'(lambda (x y) (send self :is-worse x y))))))
  265.  
  266. ;;;
  267. ;;; Other Methods Using List Representation of SImplex
  268. ;;;
  269.  
  270. (defmeth simplex-proto :best-point () (car (last (send self :simplex))))
  271. (defmeth simplex-proto :worst-point () (first (send self :simplex)))
  272. (defmeth simplex-proto :second-worst-point () (second (send self :simplex)))
  273. (defmeth simplex-proto :replace-point (new old)
  274.   (let* ((simplex (send self :simplex))
  275.          (n (position old simplex)))
  276.     (when n 
  277.           (setf (nth n simplex) new)
  278.           (send self :sort-simplex))))
  279. (defmeth simplex-proto :mean-opposite-face (x)
  280.   (let ((face (mapcar #'(lambda (x) (send self :point-location x))
  281.                       (remove x (send self :simplex)))))
  282.     (/ (apply #'+ face) (length face))))
  283.  
  284. ;;;
  285. ;;; Iteration Step Methods
  286. ;;;
  287.  
  288. (defmeth simplex-proto :extrapolate-from-worst (fac)
  289.   (let* ((worst (send self :worst-point))
  290.          (wloc (send self :point-location worst))
  291.          (delta (- (send self :mean-opposite-face worst) wloc))
  292.          (new (send self :make-point (+ wloc (* (- 1 fac) delta)))))
  293.     (if (send self :is-worse worst new) (send self :replace-point new worst))
  294.     new))
  295.  
  296. (defmeth simplex-proto :shrink-to-best (fac)
  297.   (let* ((best (send self :best-point))
  298.          (bloc (send self :point-location best)))
  299.     (dolist (x (copy-list (send self :simplex)))
  300.             (if (not (eq x best))
  301.                 (send self :replace-point 
  302.                       (send self :make-point 
  303.                             (+ bloc 
  304.                                (* fac 
  305.                                   (- (send self :point-location x) bloc))))
  306.                       x)))))
  307.  
  308. (defmeth simplex-proto :relative-range ()
  309.   (let ((best (send self :point-value (send self :best-point)))
  310.         (worst (send self :point-value (send self :worst-point))))
  311.     (* 2 (/ (abs (- best worst)) (+ 1 (abs best) (abs worst))))))
  312.