Skip to content
Browse files

Added Chebyshev approximation.

  • Loading branch information...
1 parent d43e34a commit cade7cb7730b8eed997495d68ef5d541484dd42d @tpapp committed Apr 13, 2012
Showing with 188 additions and 11 deletions.
  1. +4 −2 cl-num-utils.asd
  2. +124 −0 src/chebyshev.lisp
  3. +19 −8 src/interval.lisp
  4. +8 −1 src/package.lisp
  5. +33 −0 tests/chebyshev.lisp
View
6 cl-num-utils.asd
@@ -39,7 +39,8 @@
(:file "optimization")
(:file "differentiation")
(:file "rootfinding")
- (:file "quadrature"))))
+ (:file "quadrature")
+ (:file "chebyshev"))))
:depends-on (alexandria iterate let-plus anaphora))
(defsystem :cl-num-utils-tests
@@ -75,6 +76,7 @@
;; (:file "interactions")
(:file "differentiation")
(:file "rootfinding")
- (:file "quadrature"))))
+ (:file "quadrature")
+ (:file "chebyshev"))))
:depends-on
(iterate metabang-bind anaphora lift alexandria cl-num-utils))
View
124 src/chebyshev.lisp
@@ -0,0 +1,124 @@
+;;; -*- Mode:Lisp; Syntax:ANSI-Common-Lisp; Coding:utf-8 -*-
+
+(in-package #:cl-num-utils)
+
+(deftype double-float-vector (&optional (length '*))
+ `(simple-array double-float (,length)))
+
+(declaim (inline chebyshev-recursion))
+(defun chebyshev-recursion (x value previous-value)
+ "Chebyshev polynomial recursion formula."
+ (- (* 2 x value) previous-value))
+
+(declaim (inline chebyshev-root))
+(defun chebyshev-root (m i)
+ "Return the iTH root of the Mth Chebyshev polynomial as double-float."
+ (assert (within? 0 i m))
+ (- (cos (/ (* (+ i 1/2) (float pi 1d0)) m))))
+
+(defun chebyshev-roots (m)
+ "Return the roots of the Mth Chebyshev polynomial as a vector of
+double-floats."
+ (aprog1 (make-array m :element-type 'double-float)
+ (dotimes (i m)
+ (setf (aref it i) (chebyshev-root m i)))))
+
+(defun chebyshev-regression (f n-polynomials
+ &optional (n-points n-polynomials))
+ "Chebyshev polynomial regression using the given number of polynomials and
+points (zeroes of the corresponding Chebyshev polynomial)."
+ (check-types (n-polynomials n-points) positive-fixnum)
+ (assert (<= n-polynomials n-points) ()
+ "Can't identify ~A coefficients with only ~A points."
+ n-polynomials n-points)
+ (locally (declare (optimize speed)
+ (type positive-fixnum n-polynomials n-points))
+ (let+ ((z (the double-float-vector (chebyshev-roots n-points)))
+ (f-at-z (map 'double-float-vector
+ (lambda (z) (coerce (funcall f z) 'double-float)) z))
+ (coefficients (make-array n-points :element-type 'double-float))
+ (values z)
+ previous-values
+ ((&flet weighted-sum (values)
+ (/ (loop for v across values
+ for f across f-at-z
+ summing (* f v))
+ (/ n-points 2)))))
+ (declare (type double-float-vector
+ z f-at-z values previous-values))
+ (loop for j from 0 below n-polynomials
+ do (setf (aref coefficients j)
+ (if (zerop j)
+ (/ (reduce #'+ f-at-z) n-points)
+ (progn
+ (cond
+ ((= j 1) (weighted-sum z))
+ ((= j 2) (setf values
+ (map 'double-float-vector
+ (lambda (z)
+ (chebyshev-recursion z z 1d0))
+ z)))
+ ((= j 3)
+ (setf previous-values values
+ values (map 'double-float-vector
+ #'chebyshev-recursion
+ z previous-values z)))
+ (t (map-into previous-values
+ #'chebyshev-recursion z values previous-values)
+ (rotatef values previous-values)))
+ (weighted-sum values)))))
+ coefficients)))
+
+(defun chebyshev-evaluate (coefficients x)
+ "Return the sum of Chebyshev polynomials, weighted by COEFFICIENTS, at X."
+ (let ((value (coerce x 'double-float))
+ (previous-value 1d0)
+ (sum 0d0))
+ (dotimes (index (length coefficients))
+ (incf sum (* (aref coefficients index)
+ (cond
+ ((= index 0) 1d0)
+ ((= index 1) x)
+ (t (setf previous-value (chebyshev-recursion x value previous-value))
+ (rotatef value previous-value)
+ value)))))
+ sum))
+
+
+
+(declaim (inline ab-to-cinf cinf-to-ab))
+
+(defun cinf-to-ab (x a b c)
+ "Map x in [c,plus-infinity) to z in [a,b] using x -> (x-c)/(1+x-c)+(b-a)+a."
+ (let ((xc (- x c)))
+ (assert (<= 0 xc) () "Value outside domain.")
+ (+ (* (/ xc (1+ xc)) (- b a)) a)))
+
+(defun ab-to-cinf (z a b c)
+ "Inverse of cinf-to-ab."
+ (let ((z-norm (/ (- z a) (- b a))))
+ (assert (within? 0 z-norm 1) () "Value outside domain.")
+ (+ c (/ z-norm (- 1 z-norm)))))
+
+(defun chebyshev-approximate (f interval n-polynomials
+ &key (n-points n-polynomials) closed-left? closed-right?)
+ "Return a closure approximating F on the given INTERVAL (may be infinite on
+either end) using the given number of Chebyshev polynomials."
+ (chebyshev-approximate-implementation f interval n-polynomials n-points
+ closed-left? closed-right?))
+
+(defgeneric chebyshev-approximate-implementation (f interval n-polynomials n-points
+ closed-left? closed-right?)
+ (:documentation "Implementation of CHEBYSHEV-APPROXIMATE.")
+ (:method (f (interval plusinf-interval) n-polynomials n-points
+ closed-left? closed-right?)
+ (let+ ((a (if closed-left?
+ (chebyshev-root n-points 0)
+ -1d0))
+ (left (coerce (plusinf-interval-left interval) 'double-float))
+ (coefficients
+ (chebyshev-regression (lambda (z)
+ (funcall f (ab-to-cinf z a 1d0 left)))
+ n-polynomials n-points)))
+ (lambda (x)
+ (chebyshev-evaluate coefficients (cinf-to-ab x a 1d0 left))))))
View
27 src/interval.lisp
@@ -2,26 +2,37 @@
(in-package #:cl-num-utils)
-;;;; An interval is an ordered pair of real numbers. It is not
-;;;; necessarily decreasing, as there can be negative intervals (eg
-;;;; for reverse plots), but some functions (eg interval-containing
-;;;; and interval-intersection) return positive intervals by
-;;;; construction.
+;;; TODO: rewrite interface
+;;; TODO: open/closed, general accessors LEFT, RIGHT, CLOSED-LEFT? CLOSED-RIGHT?
(defstruct interval
"A pair of numbers designating an interval on the real line. Using the
constructor INTERVAL, LEFT <= RIGHT is enforced."
(left 0 :type real :read-only t)
(right 0 :type real :read-only t))
+(defstruct plusinf-interval
+ "Interval [left,∞)."
+ (left nil :type real :read-only t))
+
+(defstruct minusinf-interval
+ "Interval (-∞,right]."
+ (right nil :type real :read-only t))
+
(define-structure-let+ (interval) left right)
(declaim (inline interval))
-
(defun interval (left right)
"Create an INTERVAL."
- (assert (<= left right) ())
- (make-interval :left left :right right))
+ (cond
+ ((and (typep left 'real) (typep right 'real))
+ (assert (<= left right) ())
+ (make-interval :left left :right right))
+ ((and (typep left 'real) (eq right t))
+ (make-plusinf-interval :left left))
+ ((and (typep right 'real) (eq left nil))
+ (make-minusinf-interval :right right))
+ (t (error "not implemented / under construction"))))
(defun interval-length (interval)
"Difference between left and right."
View
9 src/package.lisp
@@ -295,4 +295,11 @@
#:row
#:row-with-type
#:normalize1
- #:abs-diff))
+ #:abs-diff
+ #:chebyshev-recursion
+ #:chebyshev-root
+ #:chebyshev-roots
+ #:chebyshev-regression
+ #:chebyshev-evaluate
+ #:chebyshev-approximate
+ #:chebyshev-approximate-implementation))
View
33 tests/chebyshev.lisp
@@ -0,0 +1,33 @@
+;;; -*- Mode:Lisp; Syntax:ANSI-Common-Lisp; Coding:utf-8 -*-
+
+(in-package #:cl-num-utils-tests)
+
+(deftestsuite chebyshev-tests (cl-num-utils-tests)
+ ()
+ (:equality-test #'==))
+
+(defun maximum-on-grid (f interval &optional (n-grid 1000))
+ (loop for index below n-grid
+ maximizing (funcall f
+ (interval-midpoint interval
+ (/ index (1- n-grid))))))
+
+(defun approximation-error (f f-approx interval &optional (n-grid 1000))
+ (maximum-on-grid (lambda (x)
+ (abs-diff (funcall f x) (funcall f-approx x)))
+ interval n-grid))
+
+(defun test-chebyshev-approximate (f interval n-polynomials test-interval
+ &rest rest)
+ (let ((f-approx (apply #'chebyshev-approximate f interval n-polynomials rest)))
+ (approximation-error f f-approx test-interval)))
+
+(addtest (chebyshev-tests)
+ chebyshev-open-inf
+ (ensure (<= (test-chebyshev-approximate (lambda (x) (/ x (+ 4 x)))
+ (interval 2 t) 15 (interval 2 102))
+ 1e-5))
+ (ensure (<= (test-chebyshev-approximate (lambda (x) (exp (- x)))
+ (interval 0 t) 15 (interval 0 10)
+ :n-points 30)
+ 1e-4)))

0 comments on commit cade7cb

Please sign in to comment.
Something went wrong with that request. Please try again.