[Gsll-devel] Correction for integral-evaluate-spline

[Mirko Vukovic]

On Tue, Jan 27, 2009 at 10:46 PM, Liam Healy <lhealy at common-lisp.net> wrote:
> GSL's  gsl_spline_ functions do exactly what you are looking for:
> "remember" xarr and yarr so you don't have to keep sending them.
> That's why they're missing from the arglist of
> integral-evaluate-spline in GSLL.  The  gsl_interp_ functions require
> the arrays, and that's why the GSLL equivalents take them as
> arguments.
>
> Please see and compare:
> http://www.gnu.org/software/gsl/manual/html_node/Evaluation-of-Interpolating-Functions.html
> double gsl_interp_eval_integ (const gsl_interp * interp, const double
> xa[], const double ya[], double a, double b, gsl_interp_accel * acc)
> which becomes #'intergral-evaluate-interpolation in GSLL
> and
> http://www.gnu.org/software/gsl/manual/html_node/Higher_002dlevel-Interface.html
> Function: double gsl_spline_eval_integ (const gsl_spline * spline,
> double a, double b, gsl_interp_accel * acc)
> which becomes #'integral-evaluate-spline in GSLL.
> They are two different functions with two different sets of arguments,
> though they do the same thing.

OK, but is there an error in evaluation.lisp integral-evaluate-spline?

It has as arguments x, low, high, acceleration.  What is the `x' doing
there?  The call for gsl_spline_eval_integ is:
gsl_spline_eval_integ (const gsl_spline * spline, double a, double b,
gsl_interp_accel * acc).  If so, should the definition be:

(defmfun integral-evaluate-spline (spline low high acceleration)
  "gsl_spline_eval_integ"
  (((mpointer spline) :pointer) (low :double) (high :double)
   ((mpointer acceleration) :pointer))
  :c-return :double
  :documentation			; FDL
  "Find the numerical integral of an interpolated function over the
   range [low, high], using the spline object spline,
   data arrays xa and ya and the accelerator acceleration.")

(That may have been the original stumbling block that lead me to
wonder down the incorrect path)

>
> Liam
>
> On Tue, Jan 27, 2009 at 11:29 AM, Mirko Vukovic <mirko.vukovic at gmail.com> wrote:
>> Hello,
>>
>> In order to use integral-evaluate-spline, I had to modify its
>> definition.  The one below works, but I am not sure if there is a
>> better way.  Here is the definition:
>> (defmfun integral-evaluate-spline (spline xa ya low high acceleration)
>>  "gsl_spline_eval_integ"
>>  (((mpointer spline) :pointer) ((c-pointer xa) :pointer) ((c-pointer
>> ya) :pointer) (low :double) (high :double)
>>   ((mpointer acceleration) :pointer))
>>  :c-return :double
>>  :documentation                        ; FDL
>>  "Find the numerical integral of an interpolated function over the
>>   range [low, high], using the spline object spline,
>>   data arrays xa and ya and the accelerator acceleration.")
>>
>> I removed a reference to a double-float `x' and added the x&y array
>> pointers into the arguments, as they are necessary for gsl's
>> `gsl_interp_eval_integ'.
>
> They are needed for this function and they are there, in
> integral-evaluate-interpolation.  You copied
> and pasted the definition of integral-evaluate-spline, then put in
> arguments that don't belong there.

You are right.  I messed up there.

>
>> (I am really not sure what is the purpose of the ..._e functions).
>>
>> Here is a bit of test code that exercises the defintion by integrating
>> pi*sin(pi x) between 0 & 0.5 which should equal 1.0:
>>
>> (let* ((acc (make-acceleration))
>>       (xarr
>>        (make-marray
>>         'double-float
>>         :initial-contents
>>         (loop for i from 0.0d0 below 1d0 by 0.2d0
>>            collect i)))
>>       (yarr
>>        (make-marray
>>         'double-float
>>         :initial-contents
>>         (loop for x across (cl-array xarr)
>>            collect (sin (* x pi)))))
>>       (spline (make-spline *cubic-spline-interpolation* xarr yarr)))
>>  (* pi (integral-evaluate-spline spline xarr yarr 0d0 0.5d0 acc)))
>>
>> I wonder if there is a way for `make-spline' to "remember" xarr and
>> yarr so they don't have to be passed to `integral-evaluate-spline'.
>>
>> As always, thank for the heavy lifting :-)
>>
>> Mirko
>>
>> _______________________________________________
>> Gsll-devel mailing list
>> Gsll-devel at common-lisp.net
>> http://common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel
>>
>