#|

    class6.lsp       More Programming in Lisp: searching and local variables


    - Generic fixed-point function.
    
    - Improved fixed-point with a helper.
    
    - Use of a local variable in the fixed-point function.
    
    - Functors:  functions that build other functions (num. deriv)
    
    - Newton's method as a fixed point with numerical derivative.
    
    - Finding extrema without derivatives:  Golden section search.
    
|#


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Newton's method and fixed-point recursions          ;;
;;                      (application of functional prog)   ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Find square root of number, x_1 = x_0 - f(x_0)/f'(x_1)

(defun FIXED-POINT (f x)
  "Args (f x): fixed-point iteration starting from x."
  ; This version evaluates the function f twice!
  (if (< (abs (- x (funcall f x))) .001)
      x
      (fixed-point f (funcall f x))
      ))

(fixed-point    
 #'(lambda (x) (- x (/ (- (^ x 2) 3)    ;  f(x) = x^2 - 3
                       (* 2 x))))    
 3
 )


;--- Avoid the double evaluation with a helper function

(defun FP-HELPER (f x0 x1)
  (if (< (abs (- x0 x1)) .001)
      x1
      (fp-helper f x1 (funcall f x1))
      ))

(defun FIXED-POINT-2 (f x)
  "Args (f x): fixed-point iteration starting from x."
  ; This version evaluates f once, returning last function result.
  (fp-helper f x (funcall f x) )
  )


(fixed-point-2
 #'(lambda (x) (- x (/ (- (^ x 2) 3)    ;  f(x) = x^2 - 3
                       (* 2 x))))    
 3
 )

;--- Why not just store the value?

(defun FIXED-POINT-3 (f x)
  "Args (x map): fixed-point iteration starting from x."
  ; This version uses the let form to hold local value.
  (let ((x1 (funcall f x))  )
    (if (< (abs (- x x1)) .001)
        x1
        (fixed-point f x1)
        )))

(fixed-point-3
 #'(lambda (x) (- x (/ (- (^ x 2) 3)    ;  f(x) = x^2 - 3
                       (* 2 x))))    
 3
 )


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Newton's method with numerical derivatives.         ;;
;;                            (function arguments, output) ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


; --- Set the perturbation by local interval variable

(defun NUM-DERIV (f &key (eps .001))
  "Args (f ): Returns numerical derivative function."
  (let ((eps 0.001)  )
    (lambda (x) (/ (- (funcall f (+ x eps))
                        (funcall f (- x eps)))
                     (* 2 eps)
                     ))))
                   

; --- Nicer alternative: make this an optional, keyed argument.

(defun NUM-DERIV (f &key (eps .001))
  "Args (f :eps): Returns numerical derivative function, perturbed eps."
  (lambda (x) (/ (- (funcall f (+ x eps))
                    (funcall f (- x eps)))
                 (* 2 eps)
                 )))

(plot-function (num-deriv #'(lambda (x) (^ x 2))) 0 1)

(plot-function (num-deriv #'sqrt) 0 4)



; --- New Newton's method, using numerical derivatives

(defun NEWTON (f x0)
  "Args (f x0): Newton's method for finding zero of f near x0."
  (let ((fp (num-deriv f))  )
    (fixed-point #'(lambda (x) (- x (/ (funcall f  x)
                                       (funcall fp x)
                                       )))
                 x0)
    ))

(newton #'(lambda (x) (- (^ x 2) 3)) 3)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                               ;;
;;     Newton's method for extrema, with numerical derivatives.  ;;
;;                                                               ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(defun NEWTON-EXT (f x0)
  "Args (f x0): Newton's method for finding min/max of f near x0."
  ; Needs some test to check whether its a max or a min.
  (let* ((fp  (num-deriv f))
         (fpp (num-deriv fp))   )
    (fixed-point #'(lambda (x) (- x (/ (funcall fp  x)
                                       (funcall fpp x)
                                       )))
                 x0)
    ))

(newton-ext #'(lambda (x) (^ (- x 4) 2)) 8)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                               ;;
;;    Golden section search for the extrema without derivatives  ;;
;;                                   (local functions)           ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(defun GOLDEN-SECTION (f a b)
  "Args (f a b): Golden section search for max of f in [a,b]."
  (let ((gr (/ (- (sqrt 5) 1) 2))   )
    (labels ((new-x (a b)
                    (+ a (* gr gr (- b a))))
             (new-y (a b)
                    (+ a (* gr (- b a))))
             (reduce (f a x y b fx fy)
                     (if (< (- b a) .001) (* .5 (+ a b))
                         (if (> fx fy)
                             (let ((new (new-x a y))  )
                               (reduce f a new x y (funcall f new) fx))
                             (let ((new (new-y x b))  )
                               (reduce f x y new b fy (funcall f new)))
                             )))
             )
      (let ((x (new-x  a b))
            (y (new-y  a b))  )
        (reduce f a x y b (funcall f x) (funcall f y))
        ))))
        
(golden-section #'(lambda (x) (- (^ (- x 4) 2))) -5  10)