;;; -*- Mode: LISP; Syntax: Common-Lisp; -*-
;;; Barycentric interpolation representation of real-valued functions.
;;; This is a fast approximate way of computing with a real-valued
;;; function F of a single real value defined on the real interval
;;; [-1,1] by representing F at interpolated points. This is related
;;; to Chebyshev approximation of F, but is only a small piece of the
;;; total suite of programs that might be built up on this notion.
;;; Various generalizations abound, as can be seen by
;;; Trefethen's Chebfun Project, whose programs and documentation
;;; are far more elaborate, and whose papers explain how such notions can
;;; be used in many circumstances. Chebfun uses Matlab, rather than Lisp.
;;; The main function of interest is bcif, which produces an
;;; approximation of a function, f defined in Lisp R->R, of a
;;; single variable. The approximation is really only set up to be
;;; good in the interval [-1,1]. That is, f: [-1,1]->R
;;; Two input items are specified in bcif: the callable function f,
;;; and an integer k, probably less than 1024. That is
;;; (bcif f n x) computes the value of f at location x by
;;; interpolation at n points. If you have a rapidly varying function,
;;; or one that has "edges", you should use more points.
;;; Example.. (setf s (bcif(#'sin 8)), followed by (funcall s 0.5) should give
;;; approx value for sin(0.5).
;;; Here we set up a framework of the interpolation points,
;;; which we store in a global hash table *glpa*, which is filled up as needed.
;;; These key interpolation data points are known as gauss-lobatto points.
;;; Weights for interpolation are also stored.
;;; Oh, Barycentric interpolation is equivalent to Lagrange interpolation, but
;;; the computation is nicer for reasons you can read about elsewhere.
(defvar *glpa* (make-hash-table))
(defvar *wta* (make-hash-table))
(defun glp (n) ;; return a vector of Gauss Lobatto points of length n.
(or
(gethash n *glpa*) ;; if we have already computed them, just return that entry.
(setf (gethash n *glpa*) ;; otherwise compute GL points, store them, return them.
(let((h (/ pi n)))
(declare (double-float h))
(make-array n :element-type 'double-float
:initial-contents
(loop for i from 0 to (1- n) collect
(cos(* h (+ i 0.5d0)))))))))
;; Note that glp[0]=-glp[n].
;; Indeed, if n is even, glp[k]=-glp[n-k] for k=0..n/2-1
;; If n is odd, glp[(n+1)/2]=0, and glp[k]=-glp[n-k] for k=0..(n+1)/2-1
(defun wta (n);; weight array
(or
(gethash n *wta*);; if we have already computed them, just return that entry.
(setf (gethash n *wta*);; otherwise compute weights
(let((h (glp n)))
(make-array
n :element-type 'double-float
:initial-contents
(loop for j fixnum from 0 to (1- n)
collect
(loop for i fixnum from 0 to (1- n) ;; very un-lisp-like code here...
with wj double-float = (aref h j)
with res double-float = 1.0d0
finally (return (/ 1.0d0 res)) unless (= i j)
do (setf res (* res (- wj (aref h i)))))))))))
;; Note that wta[0]=-wta[n].
;; Indeed, if n is even, wta[k]=-wta[n-k] for k=0..n/2-1
;; if n is odd, wta[k]=-wta[n-k] for k=0..(n+1)/2-1
;; (The midpoint weight is NOT zero.)
;; Make an interpolating function out of f:
(defun bcif(f n)
(let ((fa
(make-array
n
:element-type 'double-float
:initial-contents (map 'simple-array f (glp n))))
(g (glp n))
(w (wta n)))
#'(lambda(x)(bcirun fa x g w))))
(defun bcirun(fa xin xa w)
;; bcirun uses an array of function-values of f, called fa.
;; and uses an array of weights. most likely already pre-computed, wta.
;; it evaluates f(xin) using barycentric interpolation from values in fa.
;; as a special trick, when xin=nil, bcirun returns the array fa, rather than a value.
;; this trick allows us to do arithmetic on interpolation arrays. See bciplus.
(if (null xin) fa;; returns fa. Useful info from this function when x=nil
(let*((sn 0d0)(sd 0d0)(term 0d0)(n (length fa))
(x (float xin 1.0d0)) ; make sure it is a double float
(nm1 (1- n)))
(declare (double-float sd sn term x)(fixnum n nm1)
(type (simple-array double-float (*)) fa xa w))
(loop for i fixnum from 0 to nm1 do
(if (= x (aref xa i)) (return-from bcirun (aref fa i)))
(setf term (/ (the double-float (aref w i)) (- x (the double-float (aref xa i)))))
(incf sn (* term (the double-float (aref fa i))))
(incf sd term))
(if (= sd 0.0d0) #+allegro #.excl::*nan-double* #-allegro (error "Division by zero in interpolation")
(/ sn sd)))))
;; that's all that is needed.
;; a totally optional example function that shows we can do functional operations
;; on bci functions. Like adding them. We use a nil argument to extract the data needed
;; from inside the function. We redo this function in a more general setting, as bci+.
;; (see below)
;; We can apply an ordinarily (numerical) lisp function of a single real to real to
;; a bci function this way. We create a new bci fun.
(defun bcifuncall(f s) ;; create h=f.s so that (funcall h x) is like (funcall f(funcall s x)).
(let ((s-fa (funcall s nil)))
(let ((fa (map 'simple-array f s-fa)) ;; compute a new set of points.
(g (glp (length s-fa)))
(w (wta (length s-fa))))
#'(lambda(x)(bcirun fa x g w )))))
;; here is the generalization to a function of 2 arguments, and then
;; we use it to add, multiply, etc
(defun bci2argfun(%%f r s) ;; try computing %%f(r,s)
(let ((r-fa (funcall r nil))
(s-fa (funcall s nil)))
(cond ((= (length r-fa)(length s-fa)) ;; compatible arrays
(let
((fa (map 'simple-array %%f r-fa s-fa))
(g (glp (length r-fa)))
(w (wta (length r-fa))))
#'(lambda(x)(bcirun fa x g w))))
(t(error "lengths do not match so we cannot use ~s on ~s and ~s as BCI forms" %%f r s)))))
;;then
(defun bci+(r s)(bci2argfun #'+ r s))
(defun bci*(r s)(bci2argfun #'* r s))
(defun bci-(r s)(bci2argfun #'- r s))
(defun bci/(r s)(bci2argfun #'/ r s))
(defun bcimax(r s)(bci2argfun #'max r s))
;; etc etc
;; we could do integration, differentiation, finding zeros. See chebfun project web page for
;; ideas.
;; an example
;; (setf s (bcif #'sin 8))
;; make s a function that evaluates sin(x) for x in [-1,1] using
;; 8 point chebyshev / barycentric interpolation. It can be
;; used this way: (funcall s 0.5) computes (sin 0.5), approximately.
;; it gives 0.47942555438789375d0 whereas (sin 0.5d0)
;; gives 0.479425538604203d0.
;; Using (bcif #'sin 16) is a more accurate function: it gives exactly the same for sin(0.5).
;; if you don't like using funcall, you can do this [also]
;; (setf (symbol-function 's) (bcif #'sin 16)) . This
;; means you can write (s 0.5) instead of (funcall s 0.5)
;; (defun h(x) (+ (cos(exp(* 3.0 x)))(sin(exp(* 2.0d0 x))) 12.00d0 (log(+ x 4.0d0))))
;; (setf hfun (bcif #'h 64))
;; (defun g(x)(funcall hfun x))
;; (defun testgh(x)(/ (abs(- (h x)(g x))) (max (abs(h x))(abs (g x)))))
#|
;; a heuristic proof that sin^2 + cos^2 =1
(setf s (bcif #'sin 16))
(setf c (bcif #'cos 16))
(setf g (bci+ (bci* c c ) (bci* s s))) ;; g is the function sin^2 + cos^2
(pprint (loop for i from -1 to 1 by 1/10 collect(funcall g i)))
;; the approximation gets bad outside [-1,1]:
(loop for i from 1 to 5 by 1/4 do (print (1- (funcall g i))))
;; to find the coefficients in the Chebyshev series needed to encode
;; some f, do something like this:
(setf r (bcif #'sin 128)) ;; best to use a power of 2 so we can use FAST discrete cosine transform
(setf rcopy (copy-seq(funcall r nil))) ;; do this because the next step will over-write its argument.
(setf r_transform (fct rcopy 128))
;; fct is fast cosine transform and is defined in file algo749.lisp.
;; The basis for the success of the Chebfun project is that one can
;; judge how many numerical points are needed to represent some
;; function sufficiently accurately by converting the sequence of
;; interpolation points to a Chebyshev series and looking at the
;; coefficients in this transformed space. If, in the Chebyshev space,
;; all coefficients of order > N are consistently at the level of
;; "noise" compared to the larger coefficients, then they may
;; plausibly be dropped entirely. Heuristics based on examining some
;; number of terms can be used to approximate this decision.
;; If the terms do not drop fast enough, we blame the function being
;; approximated: it is ill-behaved in some way. One cure is to split
;; the function into pieces and approximate each piece separately.
|#
;; if maxima..
#+GCL
(defun maxima::$bcifun(f n)
(let ((gs (intern (symbol-name (gensym "$BF")) :maxima))
(ff `(lambda(r)(float (mfuncall (quote ,f) r) 1.0d0))))
(setf (symbol-function gs)(bcif ff n))
gs))
#+GCL (defun maxima::$bcitimes(f g)
(let ((gs (intern (symbol-name (gensym "$BF")) :maxima)))
(setf (symbol-function gs)(bci2argfun #'* f g))
gs)) ;; etc etc
#|
maxima example..
load("bci");
h(x):=sin(10*x)*exp(x);
rr:bcifun(h,30);
plot2d(['rr(x),h(x)],[x,-1.5,1.5]);
plot2d('rr(x)-h(x),[x,-1,1]);
s3:bcifun(lambda([x],x^3),10);
s6: bcitimes(s3,s3);
plot2d(['s6(x),x^6],[x,-1,1]);
|#