Mercurial > octave
view scripts/polynomial/ppval.m @ 28912:0de38a6ef693
maint: Use Octave convention of space after function name in scripts dir.
* cplxpair.m, gradient.m, integral.m, integral3.m, interp2.m, interpft.m,
num2str.m, pol2cart.m, quad2d.m, quadgk.m, quadl.m, repelem.m, sph2cart.m,
convhull.m, delaunay.m, delaunayn.m, movegui.m, print_usage.m,
__strip_html_tags__.m, colormap.m, gray2ind.m, imformats.m, importdata.m,
javachk.m, condest.m, isbanded.m, krylov.m, lscov.m, rref.m, inputParser.m,
publish.m, symvar.m, validateattributes.m, ode15i.m, ode15s.m, ode23.m,
ode23s.m, ode45.m, fminbnd.m, fminsearch.m, glpk.m, qp.m, get_forge_pkg.m,
installed_packages.m, annotation.m, axis.m, camorbit.m, campos.m, camva.m,
camzoom.m, daspect.m, legend.m, pbaspect.m, __gnuplot_legend__.m, light.m,
patch.m, plotyy.m, __pie__.m, reducepatch.m, ribbon.m, streamline.m, trisurf.m,
figure.m, ndgrid.m, __gnuplot_draw_axes__.m, __opengl_print__.m, padecoef.m,
polygcd.m, ppval.m, spline.m, union.m, bicg.m, bicgstab.m, cgs.m, eigs.m,
gmres.m, pcg.m, pcr.m, __alltohandles__.m, __sprand__.m, qmr.m, tfqmr.m,
betaincinv.m, cosint.m, gammainc.m, discrete_cdf.m, discrete_inv.m,
discrete_rnd.m, base2dec.m, strtok.m, compare_plot_demos.m,
html_compare_plot_demos.m, speed.m, test.m, weboptions.m, webread.m, webwrite.m:
Use Octave convention of space after function name and before '(' in scripts
directory.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 13 Oct 2020 18:17:29 -0700 |
parents | d8318c12d903 |
children | 7854d5752dd2 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2000-2020 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{yi} =} ppval (@var{pp}, @var{xi}) ## Evaluate the piecewise polynomial structure @var{pp} at the points @var{xi}. ## ## If @var{pp} describes a scalar polynomial function, the result is an array ## of the same shape as @var{xi}. Otherwise, the size of the result is ## @code{[pp.dim, length(@var{xi})]} if @var{xi} is a vector, or ## @code{[pp.dim, size(@var{xi})]} if it is a multi-dimensional array. ## @seealso{mkpp, unmkpp, spline, pchip} ## @end deftypefn function yi = ppval (pp, xi) if (nargin != 2) print_usage (); endif if (! (isstruct (pp) && isfield (pp, "form") && strcmp (pp.form, "pp"))) error ("ppval: first argument must be a pp-form structure"); endif ## Extract info. [x, P, n, k, d] = unmkpp (pp); ## dimension checks sxi = size (xi); if (isvector (xi)) xi = xi(:).'; endif nd = length (d); ## Determine intervals. xn = numel (xi); idx = lookup (x, xi, "lr"); P = reshape (P, [d, n * k]); P = shiftdim (P, nd); P = reshape (P, [n, k, d]); Pidx = P(idx(:), :); # 2D matrix size: x = coefs*prod (d), y = prod (sxi) if (isvector (xi)) Pidx = reshape (Pidx, [xn, k, d]); Pidx = shiftdim (Pidx, 1); dimvec = [d, xn]; else Pidx = reshape (Pidx, [sxi, k, d]); Pidx = shiftdim (Pidx, length (sxi)); dimvec = [d, sxi]; endif ndv = length (dimvec); ## Offsets. dx = (xi - x(idx))(:)'; dx = repmat (dx, [prod(d), 1]); dx = reshape (dx, dimvec); dx = shiftdim (dx, ndv - 1); ## Use Horner scheme. if (k > 1) yi = shiftdim (reshape (Pidx(1,:), dimvec), ndv - 1); else yi = shiftdim (reshape (Pidx, dimvec), ndv - 1); endif for i = 2 : k yi .*= dx; yi += shiftdim (reshape (Pidx(i,:), dimvec), ndv - 1); endfor ## Adjust shape. if ((numel (xi) > 1) || (length (d) == 1)) yi = reshape (shiftdim (yi, 1), dimvec); endif if (isvector (xi) && (d == 1)) yi = reshape (yi, sxi); elseif (isfield (pp, "orient") && strcmp (pp.orient, "first")) yi = shiftdim (yi, nd); endif if (d == 1) yi = reshape (yi, sxi); endif endfunction %!shared b, c, pp, pp2, xi, abserr %! b = 1:3; %! c = ones (2); %! pp = mkpp (b, c); %! abserr = 1e-14; %! pp2 = mkpp (b, [c;c], 2); %! xi = [1.1 1.3 1.9 2.1]; %! %!assert (ppval (pp, 1.1), 1.1, abserr) %!assert (ppval (pp, 2.1), 1.1, abserr) %!assert (ppval (pp, xi), [1.1 1.3 1.9 1.1], abserr) %!assert (ppval (pp, xi.'), [1.1 1.3 1.9 1.1].', abserr) %!assert (ppval (pp2, 1.1), [1.1;1.1], abserr) %!assert (ppval (pp2, 2.1), [1.1;1.1], abserr) %!assert (ppval (pp2, xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr) %!assert (ppval (pp2, xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr) %!assert (size (ppval (pp2, [xi;xi])), [2 2 4]) %!assert (ppval (mkpp([0 1],1), magic (3)), ones(3,3)) %! %!test %! breaks = [0, 1, 2, 3]; %! coefs = rand (6, 4); %! pp = mkpp (breaks, coefs, 2); %! ret = zeros (2, 4, 2); %! ret(:,:,1) = ppval (pp, breaks'); %! ret(:,:,2) = ppval (pp, breaks'); %! assert (ppval (pp, [breaks',breaks']), ret); ## Test input validation %!error <Invalid call> ppval () %!error <Invalid call> ppval (1) %!error <argument must be a pp-form structure> ppval (1,2) %!error <argument must be a pp-form structure> ppval (struct ("a", 1), 2) %!error <argument must be a pp-form structure> ppval (struct ("form", "ab"), 2)