Mercurial > hg > octave-terminal
changeset 9768:31900e17b5f5
improve Matlab compatibility & performance of ppval/mkpp and some associated funcs
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 02 Nov 2009 14:04:06 +0100 |
parents | 0df32e0b2074 |
children | 9a1c4fe44af8 |
files | scripts/ChangeLog scripts/general/__splinen__.m scripts/polynomial/mkpp.m scripts/polynomial/pchip.m scripts/polynomial/ppval.m scripts/polynomial/spline.m scripts/polynomial/unmkpp.m |
diffstat | 7 files changed, 82 insertions(+), 43 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,12 @@ +2009-11-02 Jaroslav Hajek <highegg@gmail.com> + + * polynomial/mkpp.m: Improve Matlab compatibility. + * polynomial/ppval.m: Ditto. + * polynomial/unmkpp: Update. + * polynomial/pchip.m: Update and optimize. + * polynomial/spline.m: Update. + * general/__splinen__.m: Update. + 2009-10-23 Jaroslav Hajek <highegg@gmail.com> * general/tril.m, general/triu.m: Remove sources.
--- a/scripts/general/__splinen__.m +++ b/scripts/general/__splinen__.m @@ -37,7 +37,7 @@ endif yi = y; for i = length(x):-1:1 - yi = permute (spline (x{i}, yi, xi{i}), [length(x),1:length(x)-1]); + yi = permute (spline (x{i}, yi, xi{i}(:)), [length(x),1:length(x)-1]); endfor [xi{:}] = ndgrid (cellfun (@(x) x(:), xi, "UniformOutput", false){:});
--- a/scripts/polynomial/mkpp.m +++ b/scripts/polynomial/mkpp.m @@ -27,14 +27,10 @@ ## lowest. There must be one row for each interval in @var{x}, so ## @code{rows (@var{p}) == length (@var{x}) - 1}. ## -## You can concatenate multiple polynomials of the same order over the -## same set of intervals using @code{@var{p} = [ @var{p1}; @var{p2}; -## @dots{}; @var{pd} ]}. In this case, @code{rows (@var{p}) == @var{d} -## * (length (@var{x}) - 1)}. -## -## @var{d} specifies the shape of the matrix @var{p} for all except the -## last dimension. If @var{d} is not specified it will be computed as -## @code{round (rows (@var{p}) / (length (@var{x}) - 1))} instead. +## @var{p} may also be a multi-dimensional array, specifying a vector-valued +## or array-valued polynomial. The shape is determined by @var{d}. If @var{d} is +## not given, the default is @code{size (p)(1:end-2)}. If @var{d} is given, the +## leading dimensions of @var{p} are reshaped to conform to @var{d}. ## ## @seealso{unmkpp, ppval, spline} ## @end deftypefn @@ -44,14 +40,24 @@ print_usage (); endif pp.x = x(:); - pp.P = P; - pp.n = length (x) - 1; - pp.k = columns (P); + n = length (x) - 1; + if (n < 1) + error ("mkpp: at least one interval is needed"); + endif + nd = ndims (P); + k = size (P, nd); if (nargin < 3) - d = round (rows (P) / pp.n); + if (nd == 2) + d = 1; + else + d = prod (size (P)(1:nd-1)); + endif endif pp.d = d; - if (pp.n * prod (d) != rows (P)) + pp.P = P = reshape (P, prod (d), [], k); + pp.orient = 0; + + if (size (P, 2) != n) error ("mkpp: num intervals in x doesn't match num polynomials in P"); endif endfunction
--- a/scripts/polynomial/pchip.m +++ b/scripts/polynomial/pchip.m @@ -115,7 +115,7 @@ c2 = -c3 - del1; c3 = c3 / h; - coeffs = [c3.'(:), c2.'(:), d1.'(:), f1.'(:)]; + coeffs = cat (3, c3, c2, d1, f1); pp = mkpp (x, coeffs, szy(1:end-1)); if (nargin == 2)
--- a/scripts/polynomial/ppval.m +++ b/scripts/polynomial/ppval.m @@ -19,9 +19,11 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi}) ## Evaluate piece-wise polynomial @var{pp} at the points @var{xi}. -## If @code{@var{pp}.d} is a scalar greater than 1, or an array, -## then the returned value @var{yi} will be an array that is -## @code{d1, d1, @dots{}, dk, length (@var{xi})]}. +## If @var{pp} is scalar-valued, the result is an array of the same shape as @var{xi}. +## Otherwise, the size of the result is @code{[pp.d, length(@var{xi})]} if +## @var{xi} is a vector, or @code{[pp.d, size(@var{xi})]} if it is a multi-dimensional +## array. If pp.orient is 1, the dimensions are permuted as in interp1, to +## @code{[pp.d, length(@var{xi})]} and @code{[pp.d, size(@var{xi})]} respectively. ## @seealso{mkpp, unmkpp, spline} ## @end deftypefn @@ -33,23 +35,45 @@ if (! isstruct (pp)) error ("ppval: expects a pp structure"); endif - if (isempty (xi)) - yi = []; - else - transposed = (columns (xi) == 1); - xi = xi(:); - xn = length (xi); - idx = lookup (pp.x, xi, "lr"); - dx = (xi - pp.x(idx)).'; - dx = reshape (dx(ones(1,prod(pp.d)),:),[pp.d,xn]); - c = reshape (pp.P(:,1), pp.n, prod (pp.d)); - yi = reshape (c(idx,:).', [pp.d, xn]); - for i = 2 : pp.k; - c = reshape (pp.P(:,i), pp.n, prod (pp.d)); - yi = yi .* dx + reshape (c(idx,:).', [pp.d, xn]); - endfor - if (transposed && isscalar (pp.d) && pp.d == 1) - yi = yi.'; - endif + + ## Extract info. + x = pp.x; + P = pp.P; + d = pp.d; + k = size (P, 3); + nd = size (P, 1); + + ## Determine resulting shape. + if (d == 1) # scalar case + yisz = size (xi); + elseif (isvector (xi)) # this is special + yisz = [d, length(xi)]; + else # general + yisz = [d, size(xi)]; endif + + ## Determine intervals. + xi = xi(:); + xn = numel (xi); + + idx = lookup (x, xi, "lr"); + + ## Offsets. + dx = (xi - x(idx)).'; + dx = dx(ones (1, nd), :); # spread (do nothing in 1D) + + ## Use Horner scheme. + yi = P(:,idx,1); + for i = 2:k; + yi .*= dx; + yi += P(:,idx,i); + endfor + + ## Adjust shape. + yi = reshape (yi, yisz); + if (d != 1 && pp.orient == 1) + ## Switch dimensions to match interp1 order. + yi = shiftdim (yi, length (d)); + endif + endfunction
--- a/scripts/polynomial/spline.m +++ b/scripts/polynomial/spline.m @@ -186,7 +186,7 @@ c = c(1:n-1,:); b = b(1:n-1,:); a = a(1:n-1,:); - coeffs = [d(:), c(:), b(:), a(:)]; + coeffs = cat (3, d.', c.', b.', a.'); ret = mkpp (x, coeffs, szy(1:end-1)); if (nargin == 3)
--- a/scripts/polynomial/unmkpp.m +++ b/scripts/polynomial/unmkpp.m @@ -30,10 +30,7 @@ ## (@var{i}, :)} contains the coefficients for the polynomial over ## interval @var{i} ordered from highest to lowest. If @code{@var{d} > ## 1}, @code{@var{p} (@var{r}, @var{i}, :)} contains the coefficients for -## the r-th polynomial defined on interval @var{i}. However, this is -## stored as a 2-D array such that @code{@var{c} = reshape (@var{p} (:, -## @var{j}), @var{n}, @var{d})} gives @code{@var{c} (@var{i}, @var{r})} -## is the j-th coefficient of the r-th polynomial over the i-th interval. +## the r-th polynomial defined on interval @var{i}. ## @item @var{n} ## Number of polynomial pieces. ## @item @var{k} @@ -54,7 +51,10 @@ endif x = pp.x; P = pp.P; - n = pp.n; - k = pp.k; + n = size (P, 2); + k = size (P, 3); d = pp.d; + if (d == 1) + P = reshape (P, n, k); + endif endfunction