Re: proposal for new m-file function

[Ben Abbott]

On Mar 22, 2012, at 11:16 AM, Ben Abbott wrote:

> On Mar 22, 2012, at 10:37 AM, Ben Abbott wrote:
> 
>> On Mar 22, 2012, at 10:12 AM, Martin Helm wrote:
>> 
>>> Am 22.03.2012 13:17, schrieb Ben Abbott:
>>>> On Mar 22, 2012, at 6:31 AM, Martin Helm wrote:
>>>> 
>>>>> Am 22.03.2012 10:43, schrieb Carlo de Falco:
>>>>>> 
>>>>>> 2012/3/22 <address@hidden
>>>>>> <mailto:address@hidden>>
>>>>>> 
>>>>>> Message: 2
>>>>>> Date: Wed, 21 Mar 2012 20:33:01 -0400
>>>>>> From: Ben Abbott <address@hidden <mailto:address@hidden>>
>>>>>> To: octave maintainers mailing list <address@hidden
>>>>>> <mailto:address@hidden>>
>>>>>> Subject: proposal for new m-file function
>>>>>> Message-ID: <address@hidden
>>>>>> <mailto:address@hidden>>
>>>>>> Content-Type: text/plain; charset="us-ascii"
>>>>>> 
>>>>>> A user's question about fixed points piecewise-linear fitting led
>>>>>> to the development of a function looks to be a good fit for
>>>>>> Octave's core.
>>>>>> 
>>>>>> 
>>>>>>  
>>>>>> https://mailman.cae.wisc.edu/pipermail/help-octave/2012-March/050900.html
>>>>>> 
>>>>>> The idea was to fit a piece-wise polynomial to a set of data.
>>>>>> After some discussion between myself and Martin Helm, the attached
>>>>>> ppfit.m was produced. The name was chosen to match the existing
>>>>>> ppval().
>>>>>> 
>>>>>> The function provides a least-squares fit of a 1D interpolation
>>>>>> with specified break positions to a set of data.
>>>>>> 
>>>>>> Demos and tests are included.
>>>>>> 
>>>>>> Any concern about adding this to Octave's core ?
>>>>>> 
>>>>>> Ben
>>>>>> 
>>>>>> googling for the name 'ppfit' I found this function:
>>>>>> 
>>>>>> https://www.assembla.com/code/zaxxon_scripts/subversion/nodes/trunk/project/scripts/ppfit.m
>>>>>> 
>>>>>> which seems to do the same as yours with the additional option of
>>>>>> returning a spline with N continuous derivatives
>>>>>> Would it be possible to add that option to your code?
>>>>>> The license of the linked function looks like BSD so it should be no
>>>>>> problem to get the code from there.
>>>>>> c.
>>>>> Wouldn't it be much cleaner to add that additional fitting options
>>>>> (quadratics and higher order splines) to interp1 (without breaking
>>>>> matlab compatibility of course)?
>>>>> I think that is the place where such functionality should naturally
>>>>> live. I could look at it over the weekend and propose a patch for it.
>>>> Ok. I'll wait on your changeset.
>>>> 
>>>> Is it possible for interp1 to return the available methods and whether the 
>>>> underlying interpolants are linear functions ?
>>>> 
>>>> Ben
>>> 
>>> What do you have in mind, is something like interp1("methods") which
>>> returns some cellarray for example with the method names and a flag
>>> indicating linear/nonlinear an option?
>> 
>> Something like below?
>> 
>>      [methods, linear] = interp1 ("methods")
> 
> I just realized I can easily check if the interpolant is linear by ...
> 
>       yi = interp1 (xi, eye (numel (xi)), x, "method");
>       s = sum (yi, 2);
>       s = s / mean (s);
>       if (sqrt (mean (abs (s).^2)) < sqrt (eps))
>               islinear = true;
>       endif
> 
> So, I'll only need to have access to the different methods.
> 
>       methods =  interp1 ("methods");
> 
> Ben

Martin,

I've been giving this more thought ... I'm not sure what you had in mind, but 
after test driving the ppfit.m that Carlo linked to, I'm inclined to;

(1) Rename Jonas Lundgren's ppfit.m to __ppfit__.m and then call it from 
interp1.m and ppfit.m to obtain the piecewise-polynomials for a specified order.

(2) Modify interp1.m to allow the method to be specified by "linear", 
"nearest", "spline", "cubic", "pchip", or numeric value which would be 
interpreted as the order of polynomial used as the interpolant. Thus, y1 and y2 
below would be equivalent.

        y1 = interp1(xb, yb, x, n);

        pp = __pfit__ (xb, yb, xb, n);
        y2 = ppval (pp, x);

(3) Modify my ppfit.m to accept the same methods as interp1.m (trivial). I'd 
also like to support __ppfit__() as it is. Thus, pp1 and pp2 below would be 
equivalent.

        pp1 = __ppfit __(x, y, xb, n);

        pp2 = ppfit (x, y, xb, method, weights, "global");

Thoughts ?

Ben