;;; -*- 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]);





	
|#