Re: proposal for new m-file function

[Olaf Till]

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

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

ppfit.patch
Description: Text Data

interp1.changeset
Description: Text document

signature.asc
Description: Digital signature