#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 a) (> b c))))
(define (oscilates? a b c d e)
; (x1,a) (x2,b) (x3,c) (x4,d) (x5,e)
; points with a**** (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
**