#lang racket
(require plot)

;;;
;;; This implements plotting a 2d graph of a one variable function.
;;; Singularities are handled properly (I think).

;;; The points evaluated are chosen using an adaptive strategy:
;;;    - in areas where the graph is smooth, few points are chosen
;;;    - in areas where the graph is oscillating, more points are chosen

;;; See http://goo.gl/Hdi9q   (plot.lisp in the Maxima source)
;;; for code implementing the adaptive strategy.

;;; See chapeter 4.1 in the YACAS book of algorithms for
;;; an explanation of the algorithm:
;;; http://yacas.sourceforge.net/Algochapter4.html#c4s1



; These count functions are used to measure the
; number of evaluations of the function to be drawn.
(define count 0)
(define (reset-count) 
  (set! count 0))
(define (increase-count) 
  (set! count (+ count 1)))

(define (wrap f [value-returned-on-error #f])
  (λ (x)
    (with-handlers ([(λ e #t) (λ x value-returned-on-error)])
      (increase-count)
      (f x))))

(define (hill-or-valley? a b c)
  ; (x1,a) (x2,b) (x3,c) points with a<b<c
  ; Check whether the middle number is below 
  ; or above the two other numbers.
  (or (and (< b a) (< b c))
      (and (> b a) (> b c))))
  
(define (oscilates? a b c d e)
  ; (x1,a) (x2,b) (x3,c) (x4,d) (x5,e) 
  ; points with a<b<c<d<e  
  (or
   ; f undefined somewhere?
   (not (andmap number? (list a b c d e)))
   ; or irregular shape?
   (and (hill-or-valley? a b c)
        (hill-or-valley? b c d)
        (hill-or-valley? c d e))))

(define (smooth? fa fa1 fb fb1 fc eps)
  ; fa fa1 fb fb1 fc are function values
  ; in the points a < a1 < b < b1 < c.
  ; eps is the relative error of the numerical
  ; approximation of the integral below f.
  (define machine-epsilon 2.2204460492503131e-16)
  (and (andmap number? (list fa fa1 fb fb1 fc))
       ; compare two approximations of the area of
       ;    f(x)-min{fa, fa1,fb, fb1,fc}
       ; If the two approximations agree, 
       ; the function is smooth.
       (<= (- (abs (/ (+ fa (* -5. fa1) (* 9. fb) (* -7. fb1) (* 2. fc)) 24.))
              (* eps (- (+ (* 5. fb) (* 8. fb1) (- fc))
                        (min fa fa1 fb fb1 fc))))
           ; this handles rounding problems 
           ; see plot.list source
           (* 150 machine-epsilon))))

(define (adaptive f a b c fa fb fc d eps)
  ; d is maximal depth 
  ; eps is the relative area error (see smooth?)
  ; At most 2^d additional points will be generated  

  ; 1. a < a1 < b < b1 < c
  (define a1 (/ (+ a b) 2.0))
  (define b1 (/ (+ b c) 2.0))
  (define fa1 (f a1))
  (define fb1 (f b1))
  ; 2. 
  (if (or (<= d 0)
          (and (not (oscilates? fa fa1 fb fb1 fc)) ; 3.
               (smooth? fa fa1 fb fb1 fc eps)))    ; 4.
      ; don't refine
      (list (vector a fa) (vector a1 fa1) (vector b fb) 
            (vector b1 fb1) (vector c fc))
      ; refine
      (let* ([eps2 (* 2.0 eps)]
             [left  (adaptive f a a1 b fa fa1 fb (- d 1) eps2)]
             [right (adaptive f b b1 c fb fb1 fc (- d 1) eps2)])
        (append left (rest right)))))


(define (region f a c [d 5] [eps (/ (* 400 400))])  
  ; handles the region between a and c 
  ; a good value for eps is 1/(screen_area)
  (define b (/ (+ a c) 2.))  
  (adaptive f a b c (f a) (f b) (f c) 1 eps))

; split-between : pred xs -> (list (list x ...) ...)
;   The list xs is split into sublists.
;   For all neighbors x1 and x2 in xs, (pred x1 x2) determines whether the list is split.
;   Example:  (split-between = (list 1 1 2 3 3))  =>  '((1 1) (2) (3 3))
(define (split-between pred xs)
  (let loop ([xs xs] [ys '()] [xss '()])
    (match xs
      [(list)                 (reverse (cons (reverse ys) xss))]
      [(list x)               (reverse (cons (reverse (cons x ys)) xss))]
      [(list x1 x2 more ...)  (if (pred x1 x2) 
                                  (loop more (list x2) (cons (reverse (cons x1 ys)) xss))
                                  (loop (cons x2 more) (cons x1 ys) xss))])))

(define (plot2d unwrapped-f [x-min -5] [x-max 5] [y-min -5] [y-max 5])
  ; wrap the function to be drawn, s.t. it 
  ; returns #f in error situations
  (define f (wrap unwrapped-f #f))  
  (define (clipped-to-different-sides? p q)
    ; are the points p and q on different
    ; sides of the y-min / y-max limit ?
    (define py (vector-ref p 1))
    (define qy (vector-ref q 1))
    (or (not (number? py))
        (not (number? qy))
        (and (< py y-min) (> qy y-max))
        (and (> py y-max) (< qy y-min))))
  (define (remove-non-numbers ps)
    (filter (λ (p) (number? (vector-ref p 1))) ps))
  ; 29 is a good value according to plot.lisp
  (define number-of-regions 29) 
  ; region widh
  (define delta (/ (- x-max x-min) number-of-regions))
  ; generate points by dividing the interval
  ; from x-min to x-max into number-of-regions regions,
  ; and calling region, which calls adaptive.
  (define points
    (append*
     (for/list ([i number-of-regions])
       (define a (+ x-min (* delta i)))
       ; avoid rounding error and compute c as:
       (define c (+ x-min (* delta (+ i 1))))    
       (if (= i 0)
           ; keep the first point for the first region
           (region f a c) 
           ; otherwise remove the first point (which is 
           ; present as the end point if the preceding region)
           (rest (region f a c))))))
  ; Split the point list into groups. Split between two points 
  ; if they are on different sides of y-min and y-max.
  ; All connected points will be drawn with (lines ...).
  (define connected-points
    (split-between clipped-to-different-sides? points))
  ; Display both the adaptive plot and the original plot
  ; for visual comparision.
  (displayln
   (list
    (begin0
      (plot (list (map lines 
                       (map remove-non-numbers 
                            connected-points)))
            #:x-min x-min #:x-max x-max 
            #:y-min y-min #:y-max y-max)
      (displayln (format "adaptive number of evaluations: ~a" count))
      (reset-count))
    (plot (function (wrap unwrapped-f +inf.0) x-min x-max) 
          #:y-min y-min  #:y-max y-max)))
  (displayln (format "original number of evaluations: ~a\n\n" count))
  (reset-count))
   
; Examples
(plot2d /)     ; check division by zero and singularity at x=0
(plot2d tan)   ; check multiple singularities
(plot2d (λ (x) (- (/ x))))   ; check symmetry
(plot2d (λ (x) (/ (abs x)))) ; check connected components of same sign
(plot2d (λ (x) (* 1e6 x)))   ; check that lines with high slopes are drawn
(plot2d (λ (x) (* 1e50 x)))   ; check that lines with high slopes are drawn