Mercurial > hg > octave-lyh

comparison scripts/polynomial/ppval.m @ 12608:59e2460acae1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make piecewise polynomial (pp) functions more compatible
author Kai Habel <kai.habel@gmx.de>
date Wed, 23 Feb 2011 08:11:40 +0100
parents c792872f8942
children e81ddf9cacd5
comparison
equal deleted inserted replaced
12606:3047363c376d 12608:59e2460acae1
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi}) 20 ## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi})
21 ## Evaluate piecewise polynomial @var{pp} at the points @var{xi}. 21 ## Evaluate piece-wise polynomial structure @var{pp} at the points @var{xi}.
22 ## If @var{pp} is scalar-valued, the result is an array of the same shape as 22 ## If @var{pp} describes a scalar polynomial function, the result is an
23 ## @var{xi}. 23 ## array of the same shape as @var{xi}.
24 ## Otherwise, the size of the result is @code{[pp.d, length(@var{xi})]} if 24 ## Otherwise, the size of the result is @code{[pp.dim, length(@var{xi})]} if
25 ## @var{xi} is a vector, or @code{[pp.d, size(@var{xi})]} if it is a 25 ## @var{xi} is a vector, or @code{[pp.dim, size(@var{xi})]} if it is a
26 ## multi-dimensional array. If pp.orient is 1, the dimensions are permuted as 26 ## multi-dimensional array.
27 ##
28 ##, the dimensions are permuted as
27 ## in interp1, to 29 ## in interp1, to
28 ## @code{[pp.d, length(@var{xi})]} and @code{[pp.d, size(@var{xi})]} 30 ## @code{[pp.d, length(@var{xi})]} and @code{[pp.d, size(@var{xi})]}
29 ## respectively. 31 ## respectively.
30 ## @seealso{mkpp, unmkpp, spline} 32 ## @seealso{mkpp, unmkpp, spline, pchip, interp1}
31 ## @end deftypefn 33 ## @end deftypefn
32 34
33 function yi = ppval (pp, xi) 35 function yi = ppval (pp, xi)
34 36
35 if (nargin != 2) 37 if (nargin != 2)
36 print_usage (); 38 print_usage ();
37 endif 39 endif
38 if (! isstruct (pp)) 40 if (! isstruct (pp) && strcmp (pp.form, "pp"))
39 error ("ppval: PP must be a structure"); 41 error ("ppval: expects a pp-form structure");
40 endif 42 endif
41 43
42 ## Extract info. 44 ## Extract info.
43 x = pp.x; 45 [x, P, n, k, d] = unmkpp (pp);
44 P = pp.P; 46
45 d = pp.d; 47 ## dimension checks
46 k = size (P, 3); 48 sxi = size (xi);
47 nd = size (P, 1); 49 if (isvector (xi))
50 xi = xi(:).';
51 endif
52
53 nd = length (d);
48 54
49 ## Determine resulting shape. 55 ## Determine intervals.
50 if (d == 1) # scalar case 56 xn = numel (xi);
51 yisz = size (xi); 57 idx = lookup (x, xi, "lr");
52 elseif (isvector (xi)) # this is special 58
53 yisz = [d, length(xi)]; 59 P = reshape (P, [d, n * k]);
54 else # general 60 P = shiftdim (P, nd);
55 yisz = [d, size(xi)]; 61 P = reshape (P, [n, k, d]);
62 Pidx = P(idx(:), :);#2d matrix size x: coefs*prod(d) y: prod(sxi)
63
64 if (isvector(xi))
65 Pidx = reshape (Pidx, [xn, k, d]);
66 Pidx = shiftdim (Pidx, 1);
67 dimvec = [d, xn];
68 else
69 Pidx = reshape (Pidx, [sxi, k, d]);
70 Pidx = shiftdim (Pidx, length (sxi));
71 dimvec = [d, sxi];
72 end
73 ndv = length (dimvec);
74
75 ## Offsets.
76 dx = (xi - x(idx));
77 dx = repmat (dx, [prod(d), 1]);
78 dx = reshape (dx, dimvec);
79 dx = shiftdim (dx, ndv - 1);
80
81 ## Use Horner scheme.
82 yi = Pidx;
83 if (k > 1)
84 yi = shiftdim (reshape (Pidx(1,:), dimvec), ndv - 1);
85 endif
86
87 for i = 2 : k;
88 yi .*= dx;
89 yi += shiftdim (reshape (Pidx(i,:), dimvec), ndv - 1);
90 endfor
91
92 ## Adjust shape.
93 if ((numel (xi) > 1) || (length (d) == 1))
94 yi = reshape (shiftdim (yi, 1), dimvec);
56 endif 95 endif
57 96
58 ## Determine intervals. 97 if (isvector (xi) && (d == 1))
59 xi = xi(:); 98 yi = reshape (yi, sxi);
60 xn = numel (xi); 99 elseif (isfield (pp, "orient") && strcmp (pp.orient, "first"))
61 100 yi = shiftdim(yi, nd);
62 idx = lookup (x, xi, "lr");
63
64 ## Offsets.
65 dx = (xi - x(idx)).';
66 dx = dx(ones (1, nd), :); # spread (do nothing in 1D)
67
68 ## Use Horner scheme.
69 yi = P(:,idx,1);
70 for i = 2:k;
71 yi .*= dx;
72 yi += P(:,idx,i);
73 endfor
74
75 ## Adjust shape.
76 yi = reshape (yi, yisz);
77 if (d != 1 && pp.orient == 1)
78 ## Switch dimensions to match interp1 order.
79 yi = shiftdim (yi, length (d));
80 endif 101 endif
81 102
103 ##
104 #if (d == 1)
105 # yi = reshape (yi, sxi);
106 #endif
107
82 endfunction 108 endfunction
109
110 %!shared b,c,pp,pp2,xi,abserr
111 %! b = 1:3; c = ones(2); pp=mkpp(b,c);abserr = 1e-14;pp2=mkpp(b,[c;c],2);
112 %! xi = [1.1 1.3 1.9 2.1];
113 %!assert (ppval(pp,1.1), 1.1, abserr);
114 %!assert (ppval(pp,2.1), 1.1, abserr);
115 %!assert (ppval(pp,xi), [1.1 1.3 1.9 1.1], abserr);
116 %!assert (ppval(pp,xi.'), [1.1 1.3 1.9 1.1].', abserr);
117 %!assert (ppval(pp2,1.1), [1.1;1.1], abserr);
118 %!assert (ppval(pp2,2.1), [1.1;1.1], abserr);
119 %!assert (ppval(pp2,xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr);
120 %!assert (ppval(pp2,xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr);
121 %!assert (size(ppval(pp2,[xi;xi])), [2 2 4]);