~ubuntu-branches/ubuntu/wily/octave/wily-proposed

« back to all changes in this revision

Viewing changes to scripts/general/interp1.m

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2014-03-10 17:32:29 UTC
  • mfrom: (1.2.6)
  • Revision ID: package-import@ubuntu.com-20140310173229-000oj77i9m6tjpvd
Tags: 3.8.1-1
* Imported Upstream version 3.8.1
  + libgui/src/main-window.cc: Fix problems with gui startup and focus issues
    (Closes: #738672)
  + scripts/pkg/private/fix_depends.m: Fix installing packages where
    dependency name contains '-' (Closes: #740984)
* debian/copyright: reflect upstream changes.
* octave.postinst: remove workaround for #681064 (package dir re-creation).
  (Closes: #681461)
* Remove patches applied upstream:
  + doc-compressed-info.diff
  + kfreebsd_tcgetattr.diff
  + octave-cli-manpage.diff
  + save-uint8-in-text-format.diff
* mkoctfile is now an executable binary instead of a shell script.
  + mkoctfile-mpi.diff: adapt for new format of mkoctfile.
  + d/control: add shlibs:Depends to liboctave-dev.
* Fix statement about the GUI in NEWS.Debian.
* debian/control: don't mention the FAQ in descriptions, it has been removed.
* Adapt to the fact that JIT is now disabled by default in ./configure script.
* Disable JIT on hurd-i386, LLVM not available there.

Show diffs side-by-side

added added

removed removed

Lines of Context:
22
22
## @deftypefnx {Function File} {@var{yi} =} interp1 (@var{y}, @var{xi})
23
23
## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{method})
24
24
## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, @var{extrap})
 
25
## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, "left")
 
26
## @deftypefnx {Function File} {@var{yi} =} interp1 (@dots{}, "right")
25
27
## @deftypefnx {Function File} {@var{pp} =} interp1 (@dots{}, "pp")
26
28
##
27
 
## One-dimensional interpolation.  Interpolates to determine the value of
28
 
## @var{yi} at the points, @var{xi}.  If not specified, @var{x} is taken
29
 
## to be the indices of @var{y}.  If @var{y} is a matrix or an N-dimensional
30
 
## array, the interpolation is performed on each column of @var{y}.
 
29
## One-dimensional interpolation.
 
30
##
 
31
## Interpolate input data to determine the value of @var{yi} at the points
 
32
## @var{xi}.  If not specified, @var{x} is taken to be the indices of @var{y}.
 
33
## If @var{y} is a matrix or an N-dimensional array, the interpolation is
 
34
## performed on each column of @var{y}.
31
35
##
32
36
## Method is one of:
33
37
##
34
38
## @table @asis
35
39
## @item @qcode{"nearest"}
36
 
## Return the nearest neighbor.
 
40
## Return the nearest neighbor
37
41
##
38
42
## @item @qcode{"linear"}
39
43
## Linear interpolation from nearest neighbors
49
53
## throughout the curve
50
54
## @end table
51
55
##
52
 
## Appending '*' to the start of the above method forces @code{interp1}
 
56
## Adding '*' to the start of any method above forces @code{interp1}
53
57
## to assume that @var{x} is uniformly spaced, and only @code{@var{x}(1)}
54
58
## and @code{@var{x}(2)} are referenced.  This is usually faster,
55
59
## and is never slower.  The default method is @qcode{"linear"}.
56
60
##
57
61
## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values
58
 
## beyond the endpoints.  If @var{extrap} is a number, replace values beyond
59
 
## the endpoints with that number.  If @var{extrap} is missing, assume NA.
 
62
## beyond the endpoints using the current @var{method}.  If @var{extrap} is a
 
63
## number, then replace values beyond the endpoints with that number.  When
 
64
## unspecified, @var{extrap} defaults to NA.
60
65
##
61
 
## If the string argument @qcode{"pp"} is specified, then @var{xi} should not be
62
 
## supplied and @code{interp1} returns the piecewise polynomial that
63
 
## can later be used with @code{ppval} to evaluate the interpolation.
 
66
## If the string argument @qcode{"pp"} is specified, then @var{xi} should not
 
67
## be supplied and @code{interp1} returns a piecewise polynomial object.  This 
 
68
## object can later be used with @code{ppval} to evaluate the interpolation.
64
69
## There is an equivalence, such that @code{ppval (interp1 (@var{x},
65
70
## @var{y}, @var{method}, @qcode{"pp"}), @var{xi}) == interp1 (@var{x}, @var{y},
66
71
## @var{xi}, @var{method}, @qcode{"extrap"})}.
71
76
## right-continuous.  If @var{x} is decreasing, the default discontinuous
72
77
## interpolant is left-continuous.
73
78
## The continuity condition of the interpolant may be specified by using
74
 
## the options, @qcode{"-left"} or @qcode{"-right"}, to select a left-continuous
 
79
## the options, @qcode{"left"} or @qcode{"right"}, to select a left-continuous
75
80
## or right-continuous interpolant, respectively.
76
81
## Discontinuous interpolation is only allowed for @qcode{"nearest"} and
77
82
## @qcode{"linear"} methods; in all other cases, the @var{x}-values must be
95
100
## @end group
96
101
## @end example
97
102
##
98
 
## @seealso{interpft}
 
103
## @seealso{interpft, interp2, interp3, interpn}
99
104
## @end deftypefn
100
105
 
101
106
## Author: Paul Kienzle
118
123
  xi = [];
119
124
  ispp = false;
120
125
  firstnumeric = true;
121
 
  rightcontinuous = [];
 
126
  rightcontinuous = NaN;
122
127
 
123
128
  if (nargin > 2)
124
129
    for i = 1:length (varargin)
129
134
          extrap = "extrap";
130
135
        elseif (strcmp ("pp", arg))
131
136
          ispp = true;
132
 
        elseif (any (strcmp ({"right", "-right"}, arg)))
 
137
        elseif (strcmp (arg, "right") || strcmp (arg, "-right"))
133
138
          rightcontinuous = true;
134
 
        elseif (any (strcmp ({"left", "-left"}, arg)))
 
139
        elseif (strcmp (arg, "left") || strcmp (arg, "-left"))
135
140
          rightcontinuous = false;
136
141
        else
137
142
          method = arg;
181
186
    y = y(p,:);
182
187
  endif
183
188
 
184
 
  if (isempty (rightcontinuous))
 
189
  if (isnan (rightcontinuous))
185
190
    ## If not specified, set the continuity condition
186
191
    if (x(end) < x(1))
187
192
      rightcontinuous = false;
205
210
    jumps = x(1:end-1) == x(2:end);
206
211
    have_jumps = any (jumps);
207
212
    if (have_jumps)
208
 
      if (any (strcmp (method, {"nearest", "linear"})))
 
213
      if (strcmp (method, "linear") || strcmp (method, ("nearest")))
209
214
        if (any (jumps(1:nx-2) & jumps(2:nx-1)))
210
 
          error ("interp1: extra points in discontinuities");
 
215
          warning ("interp1: multiple discontinuities at the same X value");
211
216
        endif
212
217
      else
213
218
        error ("interp1: discontinuities not supported for method '%s'", method);
217
222
 
218
223
  ## Proceed with interpolating by all methods.
219
224
  switch (method)
 
225
 
220
226
    case "nearest"
221
227
      pp = mkpp ([x(1); (x(1:nx-1)+x(2:nx))/2; x(nx)],
222
228
                 shiftdim (y, 1), szy(2:end));
227
233
      else
228
234
        yi = ppval (pp, reshape (xi, szx));
229
235
      endif
 
236
 
230
237
    case "*nearest"
231
238
      pp = mkpp ([x(1), x(1)+[0.5:(nx-1)]*dx, x(nx)],
232
239
                 shiftdim (y, 1), szy(2:end));
236
243
      else
237
244
        yi = ppval (pp, reshape (xi, szx));
238
245
      endif
 
246
 
239
247
    case "linear"
240
 
      dy = diff (y);
241
 
      dx = diff (x);
242
 
      dx = repmat (dx, [1 size(dy)(2:end)]);
243
 
      coefs = [(dy./dx).'(:), y(1:nx-1, :).'(:)];
 
248
 
244
249
      xx = x;
245
 
 
 
250
      yy = y;
 
251
      nxx = nx;
246
252
      if (have_jumps)
247
253
        ## Omit zero-size intervals.
248
 
        coefs(jumps, :) = [];
 
254
        yy(jumps, :) = [];
249
255
        xx(jumps) = [];
 
256
        nxx = rows (xx);
250
257
      endif
251
258
 
 
259
      dy = diff (yy);
 
260
      dx = diff (xx);
 
261
      dx = repmat (dx, [1 size(dy)(2:end)]);
 
262
 
 
263
      coefs = [(dy./dx).'(:), yy(1:nxx-1, :).'(:)];
 
264
 
252
265
      pp = mkpp (xx, coefs, szy(2:end));
253
266
      pp.orient = "first";
254
267
 
286
299
          yi = shiftdim (yi, 1);
287
300
        endif
288
301
      endif
 
302
 
289
303
    case {"spline", "*spline"}
290
304
      if (nx == 2 || starmethod)
291
305
        x = linspace (x(1), x(nx), ny);
302
316
          yi = shiftdim (yi, 1);
303
317
        endif
304
318
      endif
 
319
 
305
320
    otherwise
306
321
      error ("interp1: invalid method '%s'", method);
 
322
 
307
323
  endswitch
308
324
 
309
 
  if (! ispp)
310
 
    if (! ischar (extrap))
311
 
      ## determine which values are out of range and set them to extrap,
312
 
      ## unless extrap == "extrap".
313
 
      minx = min (x(1), x(nx));
314
 
      maxx = max (x(1), x(nx));
 
325
  if (! ispp && ! ischar (extrap))
 
326
    ## determine which values are out of range and set them to extrap,
 
327
    ## unless extrap == "extrap".
 
328
    minx = min (x(1), x(nx));
 
329
    maxx = max (x(1), x(nx));
315
330
 
316
 
      outliers = xi < minx | ! (xi <= maxx); # this catches even NaNs
317
 
      if (size_equal (outliers, yi))
318
 
        yi(outliers) = extrap;
319
 
        yi = reshape (yi, szx);
320
 
      elseif (!isvector (yi))
321
 
        yi(outliers, :) = extrap;
322
 
      else
323
 
        yi(outliers.') = extrap;
324
 
      endif
 
331
    outliers = xi < minx | ! (xi <= maxx); # this even catches NaNs
 
332
    if (size_equal (outliers, yi))
 
333
      yi(outliers) = extrap;
 
334
      yi = reshape (yi, szx);
 
335
    elseif (! isvector (yi))
 
336
      yi(outliers, :) = extrap;
 
337
    else
 
338
      yi(outliers.') = extrap;
325
339
    endif
326
340
  endif
327
341
 
388
402
%! h = plot (x, interp1 (x1, y1, x), 'b', x1, y1, 'sb');
389
403
%! hold on
390
404
%! g = plot (x, interp1 (x2, y2, x), 'r', x2, y2, '*r');
391
 
%! xlim ([1 3])
392
 
%! legend ([h(1), g(1)], {'left-continous', 'right-continuous'}, ...
393
 
%!         'location', 'east')
 
405
%! axis ([0.5 3.5 -0.5 1.5])
 
406
%! legend ([h(1), g(1)], {'left-continuous', 'right-continuous'}, ...
 
407
%!         'location', 'northwest')
394
408
%! legend boxoff
395
409
%! %--------------------------------------------------------
396
 
%! % red curve is left-continuos and blue is right-continuous at x = 2
 
410
%! % red curve is left-continuous and blue is right-continuous at x = 2
397
411
 
398
412
##FIXME: add test for n-d arguments here
399
413
 
594
608
## ENDBLOCK
595
609
## ENDBLOCKTEST
596
610
 
597
 
%!# test linear extrapolation
 
611
## test extrapolation (linear)
598
612
%!assert (interp1 ([1:5],[3:2:11],[0,6],"linear","extrap"), [1, 13], eps)
599
613
%!assert (interp1 (xp, yp, [-1, max(xp)+1],"linear",5), [5, 5])
600
614
 
 
615
## Basic sanity checks
601
616
%!assert (interp1 (1:2,1:2,1.4,"nearest"), 1)
602
617
%!assert (interp1 (1:2,1:2,1.4,"linear"), 1.4)
603
618
%!assert (interp1 (1:4,1:4,1.4,"cubic"), 1.4)
604
 
%!assert (interp1 (1:2,1:2,1.1, "spline"), 1.1)
 
619
%!assert (interp1 (1:2,1:2,1.1,"spline"), 1.1)
605
620
%!assert (interp1 (1:3,1:3,1.4,"spline"), 1.4)
606
621
 
607
622
%!assert (interp1 (1:2:4,1:2:4,1.4,"*nearest"), 1)
612
627
 
613
628
%!assert (interp1 ([3,2,1],[3,2,2],2.5), 2.5)
614
629
 
615
 
%!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap"), [-2,0.5,4,3,1.5])
616
630
%!assert (interp1 ([4,4,3,2,0],[0,1,4,2,1],[1.5,4,4.5], "linear"), [1.75,1,NA])
617
631
%!assert (interp1 (0:4, 2.5), 1.5)
618
632
 
 
633
## Left and Right discontinuities
 
634
%!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap", "right"), [-8,2,4,3,1.5])
 
635
%!assert (interp1 ([1,2,2,3,4],[0,1,4,2,1],[-1,1.5,2,2.5,3.5], "linear", "extrap", "left"), [-2,0.5,1,1.5,1.5])
 
636
 
 
637
%% Test input validation
619
638
%!error interp1 ()
620
 
%!error interp1 (1,1,1, "linear")
621
 
%!error interp1 (1,1,1, "*nearest")
622
 
%!error interp1 (1,1,1, "*linear")
623
 
%!error interp1 (1:2,1:2,1, "bogus")
 
639
%!error interp1 (1,2,3,4,5,6,7)
 
640
%!error <table too short> interp1 (1,1,1, "linear")
 
641
%!error <table too short> interp1 (1,1,1, "*nearest")
 
642
%!error <table too short> interp1 (1,1,1, "*linear")
 
643
%!warning <multiple discontinuities> interp1 ([1 1 1 2], [1 2 3 4], 1);
 
644
%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "pchip")
 
645
%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "cubic")
 
646
%!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "spline")
 
647
%!error <invalid method> interp1 (1:2,1:2,1, "bogus")
624
648