|
7017
|
1 ## Copyright (C) 2000, 2006, 2007 Paul Kienzle |
|
5827
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
7016
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
5827
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
5827
|
18 |
|
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {@var{idx} =} lookup (@var{table}, @var{y}) |
|
|
21 ## Lookup values in a sorted table. Usually used as a prelude to |
|
|
22 ## interpolation. |
|
|
23 ## |
|
|
24 ## If table is strictly increasing and @code{idx = lookup (table, y)}, then |
|
|
25 ## @code{table(idx(i)) <= y(i) < table(idx(i+1))} for all @code{y(i)} |
|
|
26 ## within the table. If @code{y(i)} is before the table, then |
|
|
27 ## @code{idx(i)} is 0. If @code{y(i)} is after the table then |
|
|
28 ## @code{idx(i)} is @code{table(n)}. |
|
|
29 ## |
|
|
30 ## If the table is strictly decreasing, then the tests are reversed. |
|
|
31 ## There are no guarantees for tables which are non-monotonic or are not |
|
|
32 ## strictly monotonic. |
|
|
33 ## |
|
|
34 ## To get an index value which lies within an interval of the table, |
|
|
35 ## use: |
|
|
36 ## |
|
|
37 ## @example |
|
|
38 ## idx = lookup (table(2:length(table)-1), y) + 1 |
|
|
39 ## @end example |
|
|
40 ## |
|
|
41 ## @noindent |
|
|
42 ## This expression puts values before the table into the first |
|
|
43 ## interval, and values after the table into the last interval. |
|
|
44 ## @end deftypefn |
|
|
45 |
|
|
46 ## Changed from binary search to sort. |
|
|
47 ## Thanks to Kai Habel <kai.habel@gmx.de> for the suggestion. |
|
|
48 |
|
|
49 ## TODO: sort-based lookup is significantly slower given a large table |
|
|
50 ## TODO: and small lookup vector. This shouldn't be a problem since |
|
|
51 ## TODO: interpolation (the reason for the table lookup in the first |
|
|
52 ## TODO: place) usually involves subsampling of an existing table. The |
|
|
53 ## TODO: other use of interpolation, looking up values one at a time, is |
|
|
54 ## TODO: unfortunately significantly slower for large tables. |
|
|
55 ## TODO: sort is order O((lt+lx)*log(lt+lx)) |
|
|
56 ## TODO: search is order O(lx*log(lt)) |
|
|
57 ## TODO: Clearly, search is asymptotically better than sort, but sort |
|
|
58 ## TODO: is compiled and search is not. Could support both, or recode |
|
|
59 ## TODO: search in C++, or assume things are good enough as they stand. |
|
|
60 |
|
|
61 function idx = lookup (table, xi) |
|
|
62 if (nargin == 2) |
|
|
63 if (isempty (table)) |
|
|
64 idx = zeros (size (xi)); |
|
|
65 elseif (isvector (table)) |
|
|
66 [nr, nc] = size (xi); |
|
|
67 lt = length (table); |
|
|
68 if (table(1) > table(lt)) |
|
|
69 ## decreasing table |
|
|
70 [v, p] = sort ([xi(:); table(:)]); |
|
|
71 idx(p) = cumsum (p > nr*nc); |
|
|
72 idx = lt - idx(1:nr*nc); |
|
|
73 else |
|
|
74 ## increasing table |
|
|
75 [v, p] = sort ([table(:); xi(:) ]); |
|
|
76 idx(p) = cumsum (p <= lt); |
|
|
77 idx = idx(lt+1:lt+nr*nc); |
|
|
78 endif |
|
|
79 idx = reshape (idx, nr, nc); |
|
|
80 else |
|
|
81 error ("lookup: table must be a vector"); |
|
|
82 endif |
|
|
83 else |
|
|
84 print_usage (); |
|
|
85 endif |
|
|
86 endfunction |
|
|
87 |
|
|
88 %!assert (lookup(1:3, 0.5), 0) # value before table |
|
|
89 %!assert (lookup(1:3, 3.5), 3) # value after table error |
|
|
90 %!assert (lookup(1:3, 1.5), 1) # value within table error |
|
|
91 %!assert (lookup(1:3, [3,2,1]), [3,2,1]) |
|
|
92 %!assert (lookup([1:4]', [1.2, 3.5]'), [1, 3]'); |
|
|
93 %!assert (lookup([1:4], [1.2, 3.5]'), [1, 3]'); |
|
|
94 %!assert (lookup([1:4]', [1.2, 3.5]), [1, 3]); |
|
|
95 %!assert (lookup([1:4], [1.2, 3.5]), [1, 3]); |
|
|
96 %!assert (lookup(1:3, [3, 2, 1]), [3, 2, 1]); |
|
|
97 %!assert (lookup([3:-1:1], [3.5, 3, 1.2, 2.5, 2.5]), [0, 1, 2, 1, 1]) |
|
|
98 %!assert (isempty(lookup([1:3], []))) |
|
|
99 %!assert (isempty(lookup([1:3]', []))) |
|
|
100 %!assert (lookup(1:3, [1, 2; 3, 0.5]), [1, 2; 3, 0]); |
