Mercurial > hg > octave-nkf

comparison scripts/general/interp1.m @ 9754:4219e5cf773d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve interp1 and pchip
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 22 Oct 2009 10:12:54 +0200
parents e9dc2ed2ec0f
children 9a1c4fe44af8
comparison
equal deleted inserted replaced
9753:892e2aa7bc75 9754:4219e5cf773d
124 endfor 124 endfor
125 endif 125 endif
126 126
127 ## reshape matrices for convenience 127 ## reshape matrices for convenience
128 x = x(:); 128 x = x(:);
129 nx = size (x, 1); 129 nx = rows (x);
130 if (isvector(y) && size (y, 1) == 1) 130 if (isvector (y))
131 y = y(:); 131 y = y(:);
132 endif 132 endif
133 ndy = ndims (y);
134 szy = size (y); 133 szy = size (y);
135 ny = szy(1); 134 y = y(:,:);
136 nc = prod (szy(2:end)); 135 [ny, nc] = size (y);
137 y = reshape (y, ny, nc);
138 szx = size (xi); 136 szx = size (xi);
139 xi = xi(:); 137 xi = xi(:);
140 138
141 ## determine sizes 139 ## determine sizes
142 if (nx < 2 || ny < 2) 140 if (nx < 2 || ny < 2)
143 error ("interp1: table too short"); 141 error ("interp1: table too short");
144 endif 142 endif
145 143
146 ## determine which values are out of range and set them to extrap, 144 ## check whether x is sorted; sort if not.
147 ## unless extrap == "extrap" in which case, extrapolate them like we 145 if (! issorted (x))
148 ## should be doing in the first place. 146 [x, p] = sort (x);
149 minx = x(1); 147 y = y(p,:);
150 maxx = x(nx); 148 endif
151 if (minx > maxx) 149
152 tmp = minx; 150 starmethod = method(1) == "*";
153 minx = maxx; 151
154 maxx = tmp; 152 if (starmethod)
155 endif
156 if (method(1) == "*")
157 dx = x(2) - x(1); 153 dx = x(2) - x(1);
158 endif 154 else
159 155 if (any (x(1:nx-1) == x(2:nx)))
160 if (! pp) 156 error ("interp1: table must be strictly monotonic");
161 if (ischar (extrap) && strcmp (extrap, "extrap")) 157 endif
162 range = 1:size (xi, 1); 158 endif
163 yi = zeros (size (xi, 1), size (y, 2)); 159
164 else 160 ## Proceed with interpolating by all methods.
165 range = find (xi >= minx & xi <= maxx); 161
166 yi = extrap * ones (size (xi, 1), size (y, 2)); 162 switch (method)
167 if (isempty (range)) 163 case "nearest"
168 if (! isvector (y) && length (szx) == 2
169 && (szx(1) == 1 || szx(2) == 1))
170 if (szx(1) == 1)
171 yi = reshape (yi, [szx(2), szy(2:end)]);
172 else
173 yi = reshape (yi, [szx(1), szy(2:end)]);
174 endif
175 else
176 yi = reshape (yi, [szx, szy(2:end)]);
177 endif
178 return;
179 endif
180 xi = xi(range);
181 endif
182 endif
183
184 if (strcmp (method, "nearest"))
185 if (pp) 164 if (pp)
186 yi = mkpp ([x(1); (x(1:end-1)+x(2:end))/2; x(end)], y, szy(2:end)); 165 yi = mkpp ([x(1); (x(1:end-1)+x(2:end))/2; x(end)], y, szy(2:end));
187 else 166 else
188 idx = lookup (0.5*(x(1:nx-1)+x(2:nx)), xi) + 1; 167 idx = lookup (0.5*(x(1:nx-1)+x(2:nx)), xi) + 1;
189 yi(range,:) = y(idx,:); 168 yi = y(idx,:);
190 endif 169 endif
191 elseif (strcmp (method, "*nearest")) 170 case "*nearest"
192 if (pp) 171 if (pp)
193 yi = mkpp ([x(1); x(1)+[0.5:(ny-1)]'*dx; x(nx)], y, szy(2:end)); 172 yi = mkpp ([x(1); x(1)+[0.5:(ny-1)]'*dx; x(nx)], y, szy(2:end));
194 else 173 else
195 idx = max (1, min (ny, floor((xi-x(1))/dx+1.5))); 174 idx = max (1, min (ny, floor((xi-x(1))/dx+1.5)));
196 yi(range,:) = y(idx,:); 175 yi = y(idx,:);
197 endif 176 endif
198 elseif (strcmp (method, "linear")) 177 case "linear"
199 dy = y(2:ny,:) - y(1:ny-1,:); 178 dy = diff (y);
200 dx = x(2:nx) - x(1:nx-1); 179 dx = diff (x);
201 if (pp) 180 if (pp)
202 yi = mkpp (x, [dy./dx, y(1:end-1)], szy(2:end)); 181 yi = mkpp (x, [dy./dx, y(1:end-1)], szy(2:end));
203 else 182 else
204 ## find the interval containing the test point 183 ## find the interval containing the test point
205 idx = lookup (x, xi, "lr"); 184 idx = lookup (x, xi, "lr");
206 ## use the endpoints of the interval to define a line 185 ## use the endpoints of the interval to define a line
207 s = (xi - x(idx))./dx(idx); 186 s = (xi - x(idx))./dx(idx);
208 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); 187 yi = bsxfun (@times, s, dy(idx,:)) + y(idx,:);
209 endif 188 endif
210 elseif (strcmp (method, "*linear")) 189 case "*linear"
190 dy = diff (y);
211 if (pp) 191 if (pp)
212 dy = [y(2:ny,:) - y(1:ny-1,:)];
213 yi = mkpp (x(1) + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end)); 192 yi = mkpp (x(1) + [0:ny-1]*dx, [dy./dx, y(1:end-1)], szy(2:end));
214 else 193 else
215 ## find the interval containing the test point 194 ## find the interval containing the test point
216 t = (xi - x(1))/dx + 1; 195 t = (xi - x(1))/dx + 1;
217 idx = max(1,min(ny,floor(t))); 196 idx = max (1, min (ny - 1, floor (t)));
218 197
219 ## use the endpoints of the interval to define a line 198 ## use the endpoints of the interval to define a line
220 dy = [y(2:ny,:) - y(1:ny-1,:); y(ny,:) - y(ny-1,:)];
221 s = t - idx; 199 s = t - idx;
222 yi(range,:) = s(:,ones(1,nc)).*dy(idx,:) + y(idx,:); 200 yi = bsxfun (@times, s, dy(idx,:)) + y(idx,:);
223 endif 201 endif
224 elseif (strcmp (method, "pchip") || strcmp (method, "*pchip")) 202 case {"pchip", "*pchip"}
225 if (nx == 2 || method(1) == "*") 203 if (nx == 2 || starmethod)
226 x = linspace (x(1), x(nx), ny); 204 x = linspace (x(1), x(nx), ny);
227 endif 205 endif
228 ## Note that pchip's arguments are transposed relative to interp1 206 ## Note that pchip's arguments are transposed relative to interp1
229 if (pp) 207 if (pp)
230 yi = pchip (x.', y.'); 208 yi = pchip (x.', y.');
231 yi.d = szy(2:end); 209 yi.d = szy(2:end);
232 else 210 else
233 yi(range,:) = pchip (x.', y.', xi.').'; 211 yi = pchip (x.', y.', xi.').';
234 endif 212 endif
235 213
236 elseif (strcmp (method, "cubic") || (strcmp (method, "*cubic") && pp)) 214 case {"cubic", "*cubic"}
237 ## FIXME Is there a better way to treat pp return return and *cubic
238 if (method(1) == "*")
239 x = linspace (x(1), x(nx), ny).';
240 nx = ny;
241 endif
242
243 if (nx < 4 || ny < 4) 215 if (nx < 4 || ny < 4)
244 error ("interp1: table too short"); 216 error ("interp1: table too short");
245 endif 217 endif
246 idx = lookup (x(2:nx-1), xi, "lr"); 218
247 219 ## FIXME Is there a better way to treat pp return and *cubic
248 ## Construct cubic equations for each interval using divided 220 if (starmethod && ! pp)
249 ## differences (computation of c and d don't use divided differences 221 ## From: Miloje Makivic
250 ## but instead solve 2 equations for 2 unknowns). Perhaps 222 ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html
251 ## reformulating this as a lagrange polynomial would be more efficient. 223 t = (xi - x(1))/dx + 1;
252 i = 1:nx-3; 224 idx = max (min (floor (t), ny-2), 2);
253 J = ones (1, nc); 225 t = t - idx;
254 dx = diff (x); 226 t2 = t.*t;
255 dx2 = x(i+1).^2 - x(i).^2; 227 tp = 1 - 0.5*t;
256 dx3 = x(i+1).^3 - x(i).^3; 228 a = (1 - t2).*tp;
257 a = diff (y, 3)./dx(i,J).^3/6; 229 b = (t2 + t).*tp;
258 b = (diff (y(1:nx-1,:), 2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2; 230 c = (t2 - t).*tp/3;
259 c = (diff (y(1:nx-2,:), 1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J); 231 d = (t2 - 1).*t/6;
260 d = y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J); 232 J = ones (1, nc);
261 233
262 if (pp) 234 yi = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ...
263 xs = [x(1);x(3:nx-2)]; 235 + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:);
264 yi = mkpp ([x(1);x(3:nx-2);x(nx)], 236 else
265 [a(:), (b(:) + 3.*xs(:,J).*a(:)), ... 237 if (starmethod)
266 (c(:) + 2.*xs(:,J).*b(:) + 3.*xs(:,J)(:).^2.*a(:)), ... 238 x = linspace (x(1), x(nx), ny).';
267 (d(:) + xs(:,J).*c(:) + xs(:,J).^2.*b(:) + ... 239 nx = ny;
268 xs(:,J).^3.*a(:))], szy(2:end)); 240 endif
269 else 241
270 yi(range,:) = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ... 242 idx = lookup (x(2:nx-1), xi, "lr");
271 + c(idx,:)).*xi(:,J) + d(idx,:); 243
272 endif 244 ## Construct cubic equations for each interval using divided
273 elseif (strcmp (method, "*cubic")) 245 ## differences (computation of c and d don't use divided differences
274 if (nx < 4 || ny < 4) 246 ## but instead solve 2 equations for 2 unknowns). Perhaps
275 error ("interp1: table too short"); 247 ## reformulating this as a lagrange polynomial would be more efficient.
276 endif 248 i = 1:nx-3;
277 249 J = ones (1, nc);
278 ## From: Miloje Makivic 250 dx = diff (x);
279 ## http://www.npac.syr.edu/projects/nasa/MILOJE/final/node36.html 251 dx2 = x(i+1).^2 - x(i).^2;
280 t = (xi - x(1))/dx + 1; 252 dx3 = x(i+1).^3 - x(i).^3;
281 idx = max (min (floor (t), ny-2), 2); 253 a = diff (y, 3)./dx(i,J).^3/6;
282 t = t - idx; 254 b = (diff (y(1:nx-1,:), 2)./dx(i,J).^2 - 6*a.*x(i+1,J))/2;
283 t2 = t.*t; 255 c = (diff (y(1:nx-2,:), 1) - a.*dx3(:,J) - b.*dx2(:,J))./dx(i,J);
284 tp = 1 - 0.5*t; 256 d = y(i,:) - ((a.*x(i,J) + b).*x(i,J) + c).*x(i,J);
285 a = (1 - t2).*tp; 257
286 b = (t2 + t).*tp; 258 if (pp)
287 c = (t2 - t).*tp/3; 259 xs = [x(1);x(3:nx-2)];
288 d = (t2 - 1).*t/6; 260 yi = mkpp ([x(1);x(3:nx-2);x(nx)],
289 J = ones (1, nc); 261 [a(:), (b(:) + 3.*xs(:,J).*a(:)), ...
290 262 (c(:) + 2.*xs(:,J).*b(:) + 3.*xs(:,J)(:).^2.*a(:)), ...
291 yi(range,:) = a(:,J) .* y(idx,:) + b(:,J) .* y(idx+1,:) ... 263 (d(:) + xs(:,J).*c(:) + xs(:,J).^2.*b(:) + ...
292 + c(:,J) .* y(idx-1,:) + d(:,J) .* y(idx+2,:); 264 xs(:,J).^3.*a(:))], szy(2:end));
293 265 else
294 elseif (strcmp (method, "spline") || strcmp (method, "*spline")) 266 yi = ((a(idx,:).*xi(:,J) + b(idx,:)).*xi(:,J) ...
295 if (nx == 2 || method(1) == "*") 267 + c(idx,:)).*xi(:,J) + d(idx,:);
268 endif
269 endif
270 case {"spline", "*spline"}
271 if (nx == 2 || starmethod)
296 x = linspace(x(1), x(nx), ny); 272 x = linspace(x(1), x(nx), ny);
297 endif 273 endif
298 ## Note that spline's arguments are transposed relative to interp1 274 ## Note that spline's arguments are transposed relative to interp1
299 if (pp) 275 if (pp)
300 yi = spline (x.', y.'); 276 yi = spline (x.', y.');
301 yi.d = szy(2:end); 277 yi.d = szy(2:end);
302 else 278 else
303 yi(range,:) = spline (x.', y.', xi.').'; 279 yi = spline (x.', y.', xi.').';
304 endif 280 endif
305 else 281 otherwise
306 error ("interp1: invalid method '%s'", method); 282 error ("interp1: invalid method '%s'", method);
307 endif 283 endswitch
308 284
309 if (! pp) 285 if (! pp)
286 if (! ischar (extrap))
287 ## determine which values are out of range and set them to extrap,
288 ## unless extrap == "extrap".
289 minx = min (x(1), x(nx));
290 maxx = max (x(1), x(nx));
291
292 outliers = xi < minx | ! (xi <= maxx); # this catches even NaNs
293 yi(outliers, :) = extrap;
294 endif
295
310 if (! isvector (y) && length (szx) == 2 && (szx(1) == 1 || szx(2) == 1)) 296 if (! isvector (y) && length (szx) == 2 && (szx(1) == 1 || szx(2) == 1))
311 if (szx(1) == 1) 297 if (szx(1) == 1)
312 yi = reshape (yi, [szx(2), szy(2:end)]); 298 yi = reshape (yi, [szx(2), szy(2:end)]);
313 else 299 else
314 yi = reshape (yi, [szx(1), szy(2:end)]); 300 yi = reshape (yi, [szx(1), szy(2:end)]);