splines need work

[Paul Price]

psSpline1D has not yet crystallised. There are many questions marked by
"XXX" in the code that should have been raised a long time ago. This is not
really a bug against the psLib implementation (because the SDRS needs work as
well) but it serves as a central place where we can sort things out.

A large part of the problem came from my past misunderstanding of splines --- I
thought we could fit a spline with n knots to N points, where N > n. This is
not so, but we need N == n. I don't see much of a use for linear splines, so I
suggest we concentrate solely on cubic splines. If we want linear splines in
the future, we can make a new type, or use an enum. I think we also want to
remove the generality of the knot locations --- it's not done in general, and I
don't see at first glance how it might be done, since one needs the y values at
those points which is what we're trying to generate in the first place! These
decisions will help simplify psSpline1D.

We want to fit a cubic spline to a set of points, (x,y). Let's dispense with
the errors (yErr) because I haven't seen any implementation of splines that deal
with errors (because it's meant for interpolation, not fitting). We shouldn't
be specifying the location of knots when we allocate the spline, because they're
set according to the particular data that's used to set the spline.

We need to decide on the type of spline (natural, clamped, not-a-knot). (I am
disappointed that the current implementation does not mention what type of
spline it assumes.) Unless there's a reasoned objection, let's go with natural
splines (I figure they look smoother than not-a-knot, but not overdamped like
clamped).

Now, the current definition of psSpline1D is sufficient for our needs.

We need to change the allocator. All the allocator need do is allocate the
psSpline1D, because we dispensed with the order, and the number of spline pieces
is set when we do a psVectorFitSpline1D.

psSpline1D *psSpline1DAlloc(void);

We get rid of psSpline1DAllocGeneric --- the knots are set by psVectorFitSpline1D.

psSpline1DEval and psSpline1DEvalVector are fine, but we need to change
psVectorFitSpline1D:

psSpline1D psVectorFitSpline1D(const psVector *x, const psVector *y);

We have dispensed with yErr. There's no need to feed in a spline, because
there's nothing to set; the function will call psSpline1DAlloc, set the knots
(from the x vector, or the indices if !x), and set the polynomials.

I think that's it. Please let me know if you have concerns, questions or
recommendations. I don't think the above changes need a great deal of more
work, but it's mainly a reorganisation and simplification of what's already there.

Some useful resources:
www.amath.washington.edu/~bloss/amath352_lectures/chapter4.pdf
www.amath.washington.edu/courses/352-winter-2002/spline_note.pdf

[gusciora@…]

Everything in this bug report sounds good to me. I was never completely
comfortable with the way the spline functionality was spec'ed. I had meant to
discuss all the XXX remarks during a previous trip to UH-Manoa, but I guess
the spline stuff was passed over.

On the generality of the knot locations: as long as their is an associated
value (y-value) for every knot location (x-value), it doesn't require any
extra work to achieve full generality of knot locations. The code I wrote
does not assume the knots are evenly spaced. So, as long as we require that
the x- and y-vectors have the same length, we should be okay.

Misc comment: I did specify, though not by name (natural, clamped, not-a-knot)
what type of spline I was implementing. In psSpline.c:

The first and second derivatives at the endpoints, undefined in the SDR,
are here defined to be 0.0. They can be modified via ypo and yp1.

and

XXX: The second derivatives at the endpoints, undefined in the SDR,
are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV.

Anyway, this corresponds to a "natural spline" so it won't require any extra
work.

Additional comment: It might be worth considering whether or not to change the
definition of the psSpline typedef to hold the second derivatives of the
splines, and not the 1-D polynomials of each spline. The reason is that my
understanding of splines is derived from Numerical Recipes; their algorithm
produces only the 2nd derivatives (not the polynomials) when fitting the
spline and uses only the 2nd derivatives when evaluating the spline. I don't
know which was is better. Their code appears much more concise. However,
while the NR code which fits the splines requires less operations (since they
don't calculate the polynomials), the NR code which evaluates the splines
requires more operations.

[gusciora@…]

I think all these problems were fixed with the previous release.

[gusciora@…]

Done.

[Paul Price]

Bug has been resolved.... closing.

[Paul Price]

Bugs have been fixed... closing.