[Gsll-devel] integral-evaluate-spline crashes for xarr containing containing negative numbers

[Mirko Vukovic]

Hello.

To illustrate the problem consider the following routine

(defun spline-integration-example (&optional (lower-limit 0d0))
  "Integral of  pi * sin (pi x) between 0 and 0.5.  Should come close to 1.0"
  (let* ((acc (make-acceleration))
	 (xarr
	  (make-marray
	   'double-float
	   :initial-contents
	   (loop for i from lower-limit 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))))

Doing (spline-integration-example) gives 1.0000108...
Doing (spline-integration-example -0.1d0) causes an error:

Unhandled memory fault at #x0.
   [Condition of type SB-SYS:MEMORY-FAULT-ERROR]

Restarts:
 0: [RETRY] Retry SLIME REPL evaluation request.
 1: [ABORT] Return to SLIME's top level.
 2: [ABORT] Exit debugger, returning to top level.


What is interesting are is the environment prior to the call to
integral-evaluate-spline (#7 below):
Backtrace:
  0: (SB-SYS:MEMORY-FAULT-ERROR)
  1: ("foreign function: #x41A9A2")
  2: ("foreign function: #x41AA80")
  3: ("foreign function: #x2AE6C85E1D1B")
  4: ("foreign function: #x2AE6C85E4683")
  5: ("foreign function: #x2AE6C85E541A")
  6: ("foreign function: #x2AE6C85E5F93")
  7: (INTEGRAL-EVALUATE-SPLINE ..)
      Locals:
        SB-DEBUG::ARG-0 = #<SPLINE {10039177D1}>
        SB-DEBUG::ARG-1 = #<VECTOR-DOUBLE-FLOAT #(-0.1d0 0.1d0
0.30000000000000004d0 0.5d0 0.7d0 ..)>
        SB-DEBUG::ARG-2 = #<VECTOR-DOUBLE-FLOAT
#(-0.3090169943749474d0 0.3090169943749474d0 ..)>
        SB-DEBUG::ARG-3 = 0.5d0
        SB-DEBUG::ARG-4 = 3.395932351202108d-313
        SB-DEBUG::ARG-5 = #<ACCELERATION {1003910A71}>
  8: (SPLINE-INTEGRATION-EXAMPLE -0.1d0)
  9: (SB-INT:SIMPLE-EVAL-IN-LEXENV (SPLINE-INTEGRATION-EXAMPLE -0.1d0)
#<NULL-LEXENV>)
 --more--

Arguments 3 & 4 are incorrect !.  They should be 0d0 and 0.5d0.  I
don't have a clue how that happens.

btw, I am using my definition of integral-evaluate-spline, which I
posted earlier today:

(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.")

Thanks to the properties of integrals, this problem is not a
show-stopper.  Just shifting the integration axis will do the job.
But it is annoying :-)

Hope Liam can help finding the fix.

Mirko