Re: Fixing octave-forge spline

[Joe Koski]

Update on this problem:

Apparently I just got lucky in my initial test cases and never encountered a
three point spline fit. My suggested fix isn't correct. Looking at csape.m,
it appears to my Fortran trained eyes, that, when you have only three data
points (n = 3) to fit with the 'variational' case, at line 106 the routine
is trying to use trisolve to solve a one-by-one matrix problem. Perhaps a
simple inverse of a single coefficient is all that it needs to work properly
with three points. The routine csape.m carefully checks for the four point
case in the fit, but never three. Apparently the MATLAB interp1 routine does
this check, because it works for these cases.

Sorry for the false alarm about a potential fix. I suspect the same type of
problem occurs with the 'not-a-knot' case. I'll see if I can kluge a better
fix. All my Octave code ends up looking like Fortran.

Joe Koski

on 1/30/04 1:16 PM, Joe Koski at address@hidden wrote:

> A while back I submitted a "bug report" regarding the use of the
> octave-forge spline routine. It apparently repeated an already known
> problem. For completeness, here is my e-mail sent Dec. 29.
> 
> ================
> In substituting the octave-forge spline.m for the Matlab interp1 routine, I
> discovered that, while interp1 apparently returns a fit, csape.m crashes
> when presented with three points in the spline fit. Example:
> 
> octave:1> x=rand(1:3)
> x =
> 
> 0.45231  0.39556  0.48169
> 
> octave:2> y=rand(1:3)
> y =
> 
> 0.73379  0.42766  0.69522
> 
> octave:3> pp=spline(x,y)
> error: `ldg' undefined near line 169 column 29
> error: evaluating argument list element number 1
> error: evaluating assignment expression near line 169, column 18
> error: evaluating if command near line 59, column 3
> error: called from `csape' in file
> `/sw/share/octave/2.1.46/site/m/octave-forge/splines/csape.m'
> error: evaluating assignment expression near line 33, column 7
> error: called from `spline' in file
> `/sw/share/octave/2.1.46/site/m/octave-forge/splines/spline.m'
> error: evaluating assignment expression near line 3, column 3
> octave:3> 
> 
> When given four points, spline returns a value. For three points, shouldn't
> spline/csape return either an error message or a value for pp? Apparently
> interp1 does return pp for three points.
> 
> Joe Koski
> ==============
> 
> In order to get my data analysis routine to work, I modified spline.m as
> follows.
> 
> Original spline.m:
> =================
> ## Author:  Kai Habel <address@hidden>
> ## Date: 3. dec 2000
> ## 2001-02-08 Paul Kienzle
> ##   * copied from csapi.m
> 
> function ret = spline (x, y, xi)
> 
> ret = csape (x, y, 'not-a-knot');
> 
> if (nargin == 3)
>   ret = ppval (ret, xi);
> endif
> 
> endfunction
> =================
> 
> My replacement:
> =================
> function ret = spline (x, y, xi)
> 
> lx=length(x);
> if(lx ==3)
>   ret = csape (x, y, 'variational');
> else
>   ret = csape (x, y, 'not-a-knot');
> end
> 
> if (nargin == 3)
>   ret = ppval (ret, xi);
> endif
> 
> Endfunction
> =================
> This change got the data analysis routine to work with seemingly reasonable
> answers. My assumption (correct or incorrect) was that csape was running out
> of spline boundary conditions, and 'not-a-knot' was not a correct choice for
> a case with only three data points. I chose 'variational' because it gives a
> "natural" spline, one of the two options ("clamped" and "free") given in my
> 1985 edition of Burden and Faires. I'm most certainly not an expert on
> spline routines.
> 
> Is this the correct modification? If it is, then csape.m should trap the
> three-point 'nat-a-knot' case and replace it with 'variational' (or some
> other boundary condition) and a warning message.
> 
> Thanks for any input or suggestions.
> 
> Joe Koski
> 
> 
> 
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------
> 



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------