Re: proposal for new m-file function

[Ben Abbott]

On Mar 22, 2012, at 6:16 AM, Olaf Till wrote:

> On Thu, Mar 22, 2012 at 10:06:52AM +0100, Olaf Till wrote:
>> On Wed, Mar 21, 2012 at 08:33:01PM -0400, Ben Abbott wrote:
>>> 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
>> 
>> The loop could be vectorized:
>> 
>> --- original_ppfit.m    2012-03-22 09:18:22.000000000 +0100
>> +++ ppfit.m     2012-03-22 09:32:04.000000000 +0100
>> @@ -113,12 +113,8 @@
>>            "Input data value pairs must have the same size.")
>>   endif
>>   ## Use linear least squares for the initial guess
>> -  a = zeros (numel (x), numel (xi));
>> -  for ii = 1:numel (xi)
>> -    yi = zeros (numel (xi), 1);
>> -    yi(ii) = 1;
>> -    a(:, ii) = interp1 (xi, yi, x, method, "extrap");
>> -  endfor
>> +  yi = eye (numel (xi));
>> +  a = interp1 (xi, yi, x(:), method, "extrap");
>>   ## y = (q*r) * yi + errors
>>   ws = warning ("off", "all"); 
>>   [q, r, p] = qr (a .* w, 0);
>> 
>> but for a bug in interp1 which transposes the returned matrix if both
>> "spline" and "extrap" are given.
>> 
>> Olaf
> 
> Attached is the bugfix for interp1.m.
> 
> Olaf

Olaf, with your changeset for interp1.m applied and with or without the change 
you proposed to ppfit.m, "demo ppfit" does not work for "cubic".

Also, perhaps bsxfun() should be used to avoid the broadcasting warnings 
(unless they'll be turned off in the next release ?)

warning: operator -: automatic broadcasting operation applied
warning: product: automatic broadcasting operation applied
warning: operator -: automatic broadcasting operation applied
warning: product: automatic broadcasting operation applied
warning: operator -: automatic broadcasting operation applied
warning: product: automatic broadcasting operation applied
warning: operator -: automatic broadcasting operation applied
warning: product: automatic broadcasting operation applied
warning: operator -: automatic broadcasting operation applied
warning: product: automatic broadcasting operation applied
warning: operator -: automatic broadcasting operation applied
warning: product: automatic broadcasting operation applied

Can you take a look ?

Ben