Re: [Octave-bug-tracker] [bug #40341] Logical indexing into sparse matrices causes OOM

[David Bateman]

On 10/22/2013 08:45 PM, Jordi Gutiérrez Hermoso wrote:
> 2) numel() should have a return type of octave_idx_type and defined
> as the total number of zero and non-zero elements if the matrix and
> so limited to [0 intmax()]
> This one I disagree with. It obviously has to be true for full
> matrices because we need to be able to linearly index as much as
> numel() elements, but it's an unnatural restriction for sparse
> matrices. Instead, we should patch numel() to return double 

For both sparse and full matrices I presume and not just sparse
matrices, otherwise I don't see that this makes sense

> and handle
> the locations where this fails by letting the compiler handle this (as
> an intermediate step, we can pick another type so that there are no
> implicit cast from double to octave_idx_type and let the compiler find
> these locations for us).

This is not clear what you mean, can you give an example.

>
> The problem is that you wrote sparse with some assumptions, and when
> your assumptions are no longer valid, you're blaming the calling code
> for not following the assumptions of your internal implementation
> details. 

All code is written with assumptions. Replacing the numel return type
with a double just moves the assumption from 2^31 to 2^53. If this is
the case why bother at all, when we start compiling with 64bit integers
then all the code as it stands will have intmax defined as 2^63

> In order to fix this, you're making your assumptions leak out
> of your implementation into the calling code, breaking encapsulation,
> playing whack-a-mole with each location that violates your
> assumptions,
The assumptions are always there, even with your solution unless you
want to use an arbitrary precision variable as the return type of numel.
Either you treat these assumptions and have code that works or don't and
you don't have code that works when it runs into your assumptions. Next
you're going to tell me that with a double return type for numel we'll
never see a matrix that large and so don't have to expose "your"
assumption, ? You are then refusing to respect the first principal.

1) The row and column indices of the sparse matrix should be of the type
octave_idx_type, limited to [0 intmax()] and not artificially limited to
something smaller

And look at

http://www.cise.ufl.edu/research/sparse/matrices/list_by_dimension.html

Sorry Jordi your assumption is already passed and your "encapsulation"
is already broken. There is a matrix of order 118142155 in Tim Davis'
library or greater than 2^26 but nearer to 2^27 and guess what n^2 is
greater than bitmax("double") and of course numel will fail on this
matrix for your return type of double as well.

> and increasing the maintenance burden on everyone,
> including yourself, and introducing some new bugs as you go along (and
> fixing it with j++, seriously?).
If you mean the changeset

http://hg.savannah.gnu.org/hgweb/octave/rev/e3ac73be7894

then "seriously" how is this related to numel. It was an off-by-one
error as the code I original copied from the old version of sprand to
reinclude here was written like that, but the following code had
slightly changed. Yes j++ is the solution, though you might prefer to
change floor to ceil and a couple of other things and end up with the
same solution. How exactly does this increase the "maintenance burden"
on everyone.

"Seriously" Jordi this is just another irrelevant attack from you

>
> What are you going to do when someone actually does call numel() for a
> sparse matrix from the Octave interpreter? 
What do you mean, you already can, and it fails as you'd expect

octave:21> numel(speye(10))
ans =  100
octave:22> numel(speye(10000))
ans =  100000000
octave:23> numel(speye(100000))
error: out of memory or dimension too large for Octave's index type

and you'd treat the case numel(speye(118142155)) in what manner exactly
if numel had a return type of double ?

> You can't leak your failed
> encapsulation all the way to user-submitted code. I'm waiting for
> David Spies to run into that bug.

And you'd do what if numel returns a double type?

>
> By the way, it's been a while since you've been actively involved in
> Octave development, and I haven't contributed much myself either
> lately in actual code, but nowadays when we fix a bug, we include a
> test for it. I'll go later and add tests to each mole you've whacked,
> if you don't beat me to it.

Jordi, you're trying to give the impression that you seem to think I'm a
bad coder and don't add tests for my patches. You do realise of course
know that I ported Paul Kienzle's test/assert functions from
Octave-Forge and deprecated the DejaGnu test code and have written a
large number of tests myself (see changelog entries from 2005), even in
my recent patches when I thought it was necessary, see

http://hg.savannah.gnu.org/hgweb/octave/rev/8a54a481ecb5
http://hg.savannah.gnu.org/hgweb/octave/rev/8031fc73f291

and the original change to sprand in june also had a test for the case
discussed when creating the patch

http://hg.savannah.gnu.org/hgweb/octave/rev/8fce0ed4894a

So frankly Jordi get your facts straight before you try to attack me.
This is the third time you and I have had this flame war, and it serves
no purpose except that it allows you Jordi to express your agression
towards me.

And what is even worse Bug #40341 has nothing to do with numel. Giving
numel a double return type won't fix it. Using 64bit indexing will
however, at least until David Spies tries

n = 2^31 + 1;
s = sprandn (n, n,  1 / n);
s (s < 0) = 0;

The suggestion I made in

https://savannah.gnu.org/bugs/?40341

Would fix even that case. John please decide what you think the return
type of numel should be and I'll align behind it.

D.