[SciPy-dev] Problem in linalg functions with overwrite check

[Pearu Peterson]

On 28 Feb 2003, Travis Oliphant wrote:

> 
> I've found a potential problem in the linear algebra functions related
> to the overwrite_a parameter.

I agree that what follows is a problem.

> Many routines alter the overwrite_a parameter after conversion of the
> input object to an array.  The check is not very robust however when the
> input is an object that uses an array as it's main core (like a Matrix
> object).
> 
> For example. Code looks something like this
> 
> a1 = asarray(a)   # Convert a to an array object
> 
> overwrite_a = overwrite_a or (a1 is not a)
> 
> So, if a was not already an array then overwrite_a is modified to
> overwrite the a1 array (which is presumed to be a copy anyway).   
> 
> The problem with this strategy is that if a is an object that has an
> __array__ method.  The asarray function will call it to obtain the
> object a as an array.  Often, this call does not return a new copy of
> the data but just a reference to the array forming the underlying
> storage of a.  So, while a and a1 are different objects and so a1 is not
> a, they use the same storage space, so overwriting a1 changes the object
> a leading to strange, unexpected behavior.
> 
> I got caught while passing Matrix objects to the linalg.det routine and
> getting different answers everytime the function ran (because the
> underlying matrix was being changed everytime).  
> 
> I'm asking for suggestions on how to fix this.  
> 
> We could for example, try to call the __array__ method on our own and
> then if it succeeds never alter the overwrite_a parameter.  Or more
> elegantly perhaps we could try to write some function
> share_storage(a1,a) --- maybe partly in C  that tries to determine
> whether or not a1 and a actually share the same storage.

Correct me if I am wrong but share_storage(a1,a) will always try to call 
a.__array__ and that will result exactly a1 if succesful. Consequently,
share_storage(a1,a) returns True whenever a has __array__ method
(assuming that subsequent calls to a.__array__ will return identical
objects). So, share_storage(a1,a) boils down to checking
hasattr(a,'__array__').
This analysis suggest the following fix:

if not hasattr(a,'__array__'):
  overwrite_a = overwrite_a or (a1 is not a)

or equivalent, but more efficient one (saves few expensive attribute
lookups):

  overwrite_a = overwrite_a \
                or (a1 is not a and not hasattr(a,'__array__'))

Pearu