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