bugGNU Octave - Bugs: bug #54241, Lookup documentation

 
 

bug #54241: Lookup documentation

Submitter:  Michael Leitner <mleitner>
Submitted:  Tue 03 Jul 2018 09:25:10 AM UTC
   
 
Category:  Documentation Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Documentation
Status:  None Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

  Spam posted by anonymous
  Spam posted by anonymous
Wed 04 Jul 2018 07:14:27 PM UTC, comment #12: 


2.2 GiB
octave:1> a=rand(1e8,1);
3.0 GiB
b=a;
3.0 GiB
b=a(2:end);
3.0 GiB
b=a(1:end-1);
3.0 GiB
octave:6> b=a(1:end-1).^2;
3.8 GiB
octave:9> b=a(1:end-1);
3.0 GiB


Yeah, it's all parsing within that loop causing the delay; so you are right about the memory.

In your second analysis, though, you've overlooked the extra math needed to adjust the index from 1-start to 0-start:


octave:16> clear all
octave:17> table=sort(rand(5,1));
octave:18> y=rand(100000000,1);
octave:19> tic;for i=1:4 lookup(table(1:5),y);end;toc
Elapsed time is 11.1749 seconds.
octave:20> tic;for i=1:4 lookup(table(1:4),y);end;toc
Elapsed time is 9.89338 seconds.
octave:21> tic;for i=1:4 lookup(table(1:5),y,'r');end;toc
Elapsed time is 9.84659 seconds.
[Here's a more apples-to-apples test...]
octave:24> tic;for i=1:4 lookup(table(2:5),y)+1;end;toc
Elapsed time is 12.8138 seconds.


This is why I've been examining the actual usages of "lookup()" in the script code to see how much post-function-call math is done.  The math manipulation in C is but an extra one or two CPU cycles because the value is already in a register--very low overhead (see the 'r' option above).  But if there is post processing, that means bringing the values back from memory which adds much more overhead.

I'm attaching a small diff file that uses "lr" in the file _interp_cube_.m instead of using the table truncation technique.  (That's faster, correct?)  The only place I see this routine used is isocolors.m, and the modification doesn't seem to affect the tests or the appearance "demo isocolors".

One last comment is that I don't see why the rule regarding NaN is imposed as such, i.e.:


If 'y(i) >= table(end)' or 'isnan (y(i))' then 'idx(i)' is N.


That is, to conflate those two conditions by default is a loss of information.  Why not, by default, just map y(i) values of NaN to index values of NaN?  Then, perhaps add an option lookup(...,"nan",VAL) where VAL could be anything, e.g., -1, 0, N, NaN?

(file #44506)

Dan Sebald <sebald>
Wed 04 Jul 2018 10:46:50 AM UTC, comment #11: 

Just one point: Octave indeed does handle the case of y=x(1:end-1); efficiently. Just try the following


a=rand(1e7,1);
b=a;
b=a(2:end);
b=a(1:end-1);


and observe the memory usage while you do that (in the operating system, not with "whos"). It does not increase, that is, also when assigning subranges only a reference is stored.

Your timing is affected by the parsing of the expression. But I think the point is that we want lookup to be efficient for large y and/or table, while we do not care for small overheads that apply once per function call. But I think that my previous conclusions were premature, because:


table=sort(rand(5,1));
y=rand(100000000,1);
tic;for i=1:4 lookup(table(1:5),y);end;toc
Elapsed time is 6.47775 seconds.
tic;for i=1:4 lookup(table(1:4),y);end;toc
Elapsed time is 5.31523 seconds.
tic;for i=1:4 lookup(table(1:5),y,'r');end;toc
Elapsed time is 6.2201 seconds.


So actually what I thought should be faster (looking up into a smaller table) indeed is faster than looking up into the full table and doing post-processing. This is something that could be also improved in the code, making it smaller and executing faster.

Michael Leitner <mleitner>
Wed 04 Jul 2018 09:02:02 AM UTC, comment #10: 

There aren't that many uses of lookup.  Just do a


grep "lookup" */*/*.m
grep "lookup" */*/*/*/*.m


lookup() is being used in several places in a way similar to what you are doing by truncating the input table vector.

Your assumption about the memory efficiency of x(1:end-1) is likely incorrect.  Consider these tests:


octave:3> tic; for i=1:100000; lookup(table,y,'r'); end; toc
Elapsed time is 3.49451 seconds.
octave:4> tic; for i=1:100000; lookup(table(1:end-1),y); end; toc
Elapsed time is 5.35454 seconds.
octave:5> tic; for i=1:100000; table; end; toc
Elapsed time is 0.0756271 seconds.
octave:6> tic; for i=1:100000; table(1:end-1); end; toc
Elapsed time is 1.99329 seconds.


My guess is that to implement what you are thinking, i.e., that x(1:end-1) is simply an input and hence no new vector or parsing consideration is needed, would probably require some fair amount of work on JWE's part.  There'd have to be some kind of reference to a vector for which the referencing "sub-vector" has a start and stop index specifying a subrange.  That info would be needed for checking conformance, internal loop limits, etc.

The option "m" (and hence "b") seems like it could be accomplished with the routine ismember(), which already has a "rows" option.  I can see the "rows" option as useful for lookup().  There might be big-data users who have a need for organizing in a lexicographic manner.  Lexicographic is already used for strings in cell-lookup anyway, so maybe it's not all that difficult to implement.

Yes, your definition isn't quite correct, hence neither mine.  They are close, but also the endpoint handling is a monkey wrench in a simple definition.  There is even a conflicting aspect of the current definition with this example:


octave:25> lookup([3],5*rand(5,1),'lr')
ans =

   1
   1
   1
   1
   1


I.e., this is a case where 'r' doesn't mean all induces are at most N-1, i.e., 0.  But that's really not a case where someone would use lookup().  Plus, the case of table monotonically decreasing.

I would say:

1) Fix the "table(idx(i+1)) should be table(idx(i)+1)".

2) Consider explicitly writing out the formula and endpoints for the monotonically decreasing table scenario.  It would be a matter of convenience for the user.  E.g.,

"
If the table is decreasing, then the tests are reversed.  I.e.,
<correct formula and endpoints here>

For non-strictly monotonic tables, empty intervals are
always skipped. The result is undefined if TABLE is not monotonic, or if TABLE contains a NaN.
"

3) Implement "rows" as an option.

Dan Sebald <sebald>
Wed 04 Jul 2018 07:32:23 AM UTC, comment #9: 

Dan: interesting that you had just the same feeling with regard to 'l' and 'r' as I had. Indeed, these are just syntactic sugar (with not much energy content at that), because you do not even need to explicitly extend your table to infinity: as I said,  lookup(table,y,'l') is always equal to lookup(table(2:end),y)+1, and lookup(table,y,'r') is equal to lookup(table(1:end-1),y). And of course, lookup(table,y,'lr') is equal to lookup(table(2:end-1),y)+1. So not many bytes of code can be saved by using these options, at the expense of having to maintain the code in lookup and making the documentation more intimidating. And I too find it counterintuitive that there are options that make a function explicitly ignore an element of the input, when I could equally have just removed the element from the input (note that with Octave's memory management writing table(2:end) does not take any significant memory overhead). But if there are a number of functions that use these options, it is probably really something we will have to live with. Further: as I read the code of lookup, actually the function do_numeric_lookup in libinterp/corefcn/lookup.cc, it seems to me that the behaviour of 'l' and friends is provided by post-processing. So I would have assumed that lookup(table,y,'r') takes longer than lookup(table,y), while my proposal of lookup(table(1:end-1),y) should have been faster, because it does not do post-processing and the table is smaller. But the timing results do not confirm this, actually using 'l' and friends is faster. Does anybody understand why this is?

As regards the case of unsorted tables: the behaviour is undefined, that is, anything can happen. What probably will happen is that lookup does interval bisection and stops once it has found an idx so that table(idx)<=y<table(idx+1), which is not necessarily the first. So in this case also my definition looks incorrect, but is not, since everything is under the condition of sorted tables.

As regards my "rows" proposal: yes, in my example I define table and y, and I would like to be able to do


idx=lookup(table,y,"rows");


I can get the desired behaviour by the three lines I indicated, but that is not something the average user would come up with, and it also took me ten minutes or so. I do not see inhowfar I would need a y(i).

Michael Leitner <mleitner>
Tue 03 Jul 2018 09:08:01 PM UTC, comment #8: 

Here's one case where lookup is used with non-increasing table:


@linux ~/octave/octave/octave $ grep "lookup" */*/*.m
[snip]
scripts/statistics/discrete_inv.m:  inv(k) = v(length (p) - lookup (sort (p, "descend"), x(k)) + 1);
[snip]


I would say, looking at all the uses of lookup, there is a lot of table truncating and post-call "+1" going on.  So, one wonders if either the calls to lookup() could make use of the "lr" option, or maybe the location that is calling lookup() could use find() or ismember() creatively to achieve the same effect.

Dan Sebald <sebald>
Tue 03 Jul 2018 07:57:04 PM UTC, comment #7: 

You have broken the contract of the function. That Octave doesn't check for whether the table is sorted or not is irrelevant. It doesn't check for an error, fine. But you are calling the function incorrectly. Any call to lookup with a table argument that is not sorted in ascending order is invalid.

Mike Miller <mtmiller>
Group Member
Tue 03 Jul 2018 07:49:25 PM UTC, comment #6: 

It's not invalid.  Type


octave:93> lookup([1 5 2 4 3], 2)
ans =  3


Octave doesn't complain.

Dan Sebald <sebald>
Tue 03 Jul 2018 07:46:06 PM UTC, comment #5: 

With regard to lexicographic comparisons, I wonder if it would be useful to introduce a new function, say, lexcmp() that has a behavior like

lexcmp([3 2 1], [1 2 3]) -> 1
lexcmp([1 2 3], [3 2 1]) -> -1
lexcmp([1 2 3], [1 2 3]) -> 0
lexcmp('abc', 'xyz') -> -1
lexcmp('abc', 'abc') -> 0
etc.

It would be convenient to have such a thing.  Then one could do something like

ls =sortrow(lexdata);
lexcmp([3 2 4], lexdata);

Dan Sebald <sebald>
Tue 03 Jul 2018 07:39:41 PM UTC, comment #4: 

The very first line of the doc string is


Lookup values in a *sorted* table.


so your example is invalid.

Mike Miller <mtmiller>
Group Member
Tue 03 Jul 2018 07:20:38 PM UTC, comment #3: 

Mike wrote: "But I am also unclear on lookup, I have to relearn how it works every time I look at it"

Yes, my impression from months back was that "table lookup" is a very simple concept but after reading the documentation one is left confused.  Maybe if the definition were simpler and reflected precisely what the routine is doing, which Michael's definition is more accurate.  What I see wrong with the current documentation definition is that its definition is implicitly dependent on table being ordered.  That is "y(i) < table(idx(i)+1)" isn't necessarily true with regard to this routine.  Consider:


octave:87> y = 2; table = 5:-1:1;
octave:88> idx = lookup(table,y)
idx =  4
octave:89> y < table(idx+1)
ans = 0


With this example:


octave:79> lookup([1 5 2 4 3], 2)
ans =  3


it becomes clear that Michael's definition is more correct.

So how about we make the "increasing" aspect be an ancillary comment.  That is:

"For each element of @var{y} the function returns the index of the last element in @var{table} that is smaller than or equal to the entry in @var{y}...  If table is sorted, the result may be used as a prelude to interpolation."

Dan Sebald <sebald>
Tue 03 Jul 2018 06:27:38 PM UTC, comment #2: 

I cleaned up this documentation a bit not too long ago, but I didn't really look at the details as I didn't want to make too many changes.  I did have comments about the very things you are mentioning.  See

https://savannah.gnu.org/bugs/?func=detailitem&item_id=52785#comment4

https://savannah.gnu.org/bugs/?func=detailitem&item_id=52785#comment7

Rik mentioned keeping 'lr' because of consistency with other routines.  I sort of like leaving the infinite extensions up to the user by use of [-inf <rest of table> inf].

You are correct that it should be idx(i)+1, with consideration for the end points.  So that should be fixed.  The formula in itself is kind of confusing to me, simply because it doesn't seem like a conventional math construct.  The <= means "satisfies", I guess.  Then the "for all" is switching context to y(i) meaning essentially element-by-element basis.  Typically, when one uses the expression "for all" in math, it is part of the definition referring so some independent variable, "for all x in set U".

In a search for a more verbal and more immediately comprehensible explanation, I would suggest the singular "index of" rather than "indices" of because this is a unique value.  That is, in set operations (which is sort of what this is) we can have conceptual operations that return a unique value or multiple values, e.g., find(rand(1,5) < 0.2) returns multiple values for each unique input value, if one thinks of 0.2.  I.e., "When table is increasing, for each element of @var{y} the function returns the index of the last element in @var{table} that is smaller than or equal to the entry in @var{y}, that is, ...".  There's always ambiguity in the descriptive sentence.

Regarding the lexicographic comparison, in your example shouldn't there be something equivalent to "y(i)".  That is, since your example has 3 dimensions, shouldn't it be something like, say

lookup(table, [1 1 1], 'rows')

?  Where is that y in your example?  yd?

Dan Sebald <sebald>
Tue 03 Jul 2018 06:03:58 PM UTC, comment #1: 

Can you provide a patch or hg changeset that fixes the documentation?

Note that the lookup function is used by several Octave Forge packages (communcations, database, general, geometry, image, interval, splines, statistics), including the "b", "m", and "l" options. If the interface is changed it should be coordinated with the appropriate packages.

But I am also unclear on lookup, I have to relearn how it works every time I look at it, so I can't say if the options are useful or unnecessary at the moment.

Mike Miller <mtmiller>
Group Member
Tue 03 Jul 2018 09:25:10 AM UTC, original submission:  

The documentation of lookup is incorrect, instead of


@code{table(idx(i)) <= y(i) < table(idx(i+1))}


it should read


@code{table(idx(i)) <= y(i) < table(idx(i)+1)}


While we are at it, I would also propose to describe its behaviour in natural speech a bit better, as I find that I always have to try it with a simple example in order to understand its syntax, and whenever I indicate to someone that the function would be useful for their particular problem, they cannot understand from the documentation what it's about. So I propose "For increasing table, the function returns the indices of the last elements in @var{table} that are smaller than or equal to the entries in @var{y}, that is, ...".

By the way, I do not get at all the purpose of the options 'l' and 'r', because, unless I am very much mistaken, we have that lookup(table,y,'l') is always equal to lookup(table(2:end),y)+1, and lookup(table,y,'r') is even equal to just lookup(table(1:end-1),y).  And the standard answer that this is for compatibility with Matlab does not apply here, as lookup is not provided in Matlab.

On the other hand, I would find an option "rows", corresponding to lexicographic comparisons in analogy to unique, useful. As I see it, computing this efficiently at present would require the following evil trickery:


table=unique(round(3*rand(10,3)),"rows");
y=round(3*rand(10,3));

[yd,j]=sortrows(y);
[~,i]=sortrows([table; yd]);
idx(j,1)=cumsum(i<=size(table,1))(i>size(table,1));


Michael Leitner <mleitner>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by sebald (Posted a comment)
  • -email is unavailable- added by mleitner (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-04-26 mtmiller Carbon-CopyRemoved 80942 -
    2018-07-04 sebald Attached File- Added interp_cube_lookup_lr-djs2018jul04.diff, #44506

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code