Re: proposal for new m-file function

[Ben Abbott]

On Mar 23, 2012, at 9:23 PM, Martin Helm wrote:

> Am 24.03.2012 01:25, schrieb 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
> 
> I found today some time to test the ppfit function from JonasLundgren myself 
> and came more and more to the conclusion that I cannot beat it - not even the 
> way it can handle the pure interpolation which is included in his function as 
> a subset of its features, while working myself through some algorithms for 
> higher order splines.
> The way his implementation handles and constructs the basis for the splines 
> is very clever.
> I fully agree that what you propose is the best solution.
> What I had myself in mind does not even come near to it.

I agree. Jonas Lundgren's algorithm is impressive.

I also spent some time with Lundgren's script. I had hoped to find a way to do 
the equivalent of the command below in an efficient way.

        a = interp1 (xb, eye (numel (xb)), x, method, "extrap");

Rather than modifying his script, perhaps it is best just to do something like 
...

        yb = eye (numel (xb));
        yb = mat2cell (yb, size (yb, 1), ones (1, size (yb, 2)));
        func = @(yb) ppval (ppfit (xb, yb, xb, order), x)(:);
        a = horzcat (cellfun (func, yb, "UniformOutput", false){:});

In any event, I'll email Jonas and let him know we plan to use his script.

Ben