;; function interpolation
;; Yann LeCun, June 2006.

(libload "libnum-config")

#? **** Function Interpolation
;; This defines a class for approximating scalar functions
;; by interpolating between a number of control points.
;; This provides a variety of interpolation methods, including,
;; linear, polynomial, cubinc spline, and Akima splines.
;; Methods are available for computing integrals and 
;; first and second derivatives.
;; The class implementation uses GSL's interpolation functions.
;; Example of use:
;; {<code>
;; ;; create cubic spline interpolator
;; (defparameter z (new interpolator 2 [0 1 2 3 4] [0 5 2 8 1]))
;; = ::interpolator:82b15b8
;; ? (==> z eval 0.5)
;; = 3.6317
;; ? (==> z eval 1)
;; = 5
;; ? (==> z derivative 1)
;; = -1.0357
;; = 8.0179
;; ? (==> z derivative2 1)
;; = -18.1071
;; ? (==> z integral 0 1)
;; = 3.2545
;; ? (delete z)
;; </code>}


#? interpolator
;; The interpolator class.
(defclass interpolator object
  ((-idx1- (-double-)) xarray)
  ((-idx1- (-double-)) yarray)
  ((-gptr-) xp)
  ((-gptr-) yp)
  ((-gptr-) gslinterp)
  ((-gptr-) accel)
  )

#? (new interpolator <type> <xarray> <yarray>)
;; Create a new interpolator object.
;; Argument <type> determines the type of interpolation
;; to be performed: 0: linear, 1: polynomial, 2: cubic spline, 
;; 3: cubic spline with periodic boundary conditions,
;; 4: akima spline, 5: akima spline with periodic BC.
;; <xarray> and <yarray> are idx1 of doubles that contain
;; the abscissa and ordinate of each control point.
;; The values in <xarray> must sorted by increasing value
;; (no check is performed).
;; Different types of interpolation require a different
;; minimum number of points to function properly (e.g. Akima
;; requires 5).
(defmethod interpolator interpolator (type xa ya)
  ((-int-) type)
  ((-idx1- (-double-)) xa)
  ((-idx1- (-double-)) ya)
  (cpheader "#include <gsl/gsl_interp.h>")
  (when (<> (idx-dim xa 0) (idx-dim ya 0))
    (error "x and y arrays have different dimensions"))
  (let ((n (idx-dim xa 0)))
    ((-int-) n)
    (setq xarray (double-matrix n))
    (setq yarray (double-matrix n))
    (idx-copy xa xarray)
    (idx-copy ya yarray)
    #{{ $xp = IDX_PTR($xarray, double);
        $yp = IDX_PTR($yarray, double);
        $accel = (void *)gsl_interp_accel_alloc();
        switch($type) {
         case 0: 
           $gslinterp = (void *)gsl_interp_alloc(gsl_interp_linear, $n);
           break;
         case 1: 
           $gslinterp = (void *)gsl_interp_alloc(gsl_interp_polynomial, $n);
           break;
         case 2: 
           $gslinterp = (void *)gsl_interp_alloc(gsl_interp_cspline, $n);
           break;
         case 3: 
           $gslinterp = (void *)gsl_interp_alloc(gsl_interp_cspline_periodic, $n);
           break;
         case 4: 
           $gslinterp = (void *)gsl_interp_alloc(gsl_interp_akima, $n);
           break;
         case 5: 
           $gslinterp = (void *)gsl_interp_alloc(gsl_interp_akima_periodic, $n);
           break;
         default: run_time_error("illegal interpolation type (must be 0-5)"); 
        }
        if ($n < gsl_interp_min_size((gsl_interp *)$gslinterp)) {
          gsl_interp_accel_free( (gsl_interp_accel *)$accel);
          gsl_interp_free( (gsl_interp *)$gslinterp);
          run_time_error("not enough points for this type of interpolator");
        }
        gsl_interp_init( (gsl_interp *)$gslinterp, $xp, $yp, $n);
      }#}
    ()))

(defmethod interpolator -destructor ()
  (when (<> (to-gptr ()) gslinterp) 
    #{ gsl_interp_free( (gsl_interp *)$gslinterp) #})
  (when (<> (to-gptr ()) accel) 
    #{ gsl_interp_accel_free( (gsl_interp_accel *)$accel) #}) ())

#? (==> <interpolator> eval <x>)
;; evaluate interpolated function at abscissa <x>
(defmethod interpolator eval (x)
  ((-double-) x)
  (to-double #{ gsl_interp_eval((gsl_interp *)$gslinterp,$xp,$yp,$x,(gsl_interp_accel *)$accel) #}))

#? (==> <interpolator> derivative <x>)
;; evaluate first derivative of interpolated function at abscissa <x>
(defmethod interpolator derivative (x)
  ((-double-) x)
  (to-double #{ gsl_interp_eval_deriv((gsl_interp *)$gslinterp,$xp,$yp,$x,(gsl_interp_accel *)$accel) #}))

#? (==> <interpolator> derivative2 <x>)
;; evaluate 2nd derivative of interpolated function at abscissa <x>
(defmethod interpolator derivative2 (x)
  ((-double-) x)
  (to-double #{ gsl_interp_eval_deriv2((gsl_interp *)$gslinterp,$xp,$yp,$x,(gsl_interp_accel *)$accel) #}))

#? (==> <interpolator> integral <a> <b>)
;; evaluate integral of interpolated function between
;; abscissas <a> and <b>.
(defmethod interpolator integral (a b)
  ((-double-) a b)
  (to-double #{ gsl_interp_eval_integ((gsl_interp *)$gslinterp,$xp,$yp,$a,$b,(gsl_interp_accel *)$accel) #}))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(dhc-make-with-libs
 ()
 (list libnum-libgslcblas libnum-libgsl )
 (interpolator interpolator -destructor eval derivative derivative2 integral))


    