Re: proposal for new m-file function

[Ben Abbott]

On Mar 22, 2012, at 11:23 AM, Olaf Till wrote:

> On Thu, Mar 22, 2012 at 08:37:09AM -0400, Ben Abbott wrote:
>> 
>> 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
> 
> Should be fixed, new patches attached.
> 
> Olaf

I ran some tests, and pushed your changeset.

        http://hg.savannah.gnu.org/hgweb/octave/

Thanks

Ben