On Oct 29, 12:05 am, Manny <[EMAIL PROTECTED]> wrote:
> In sage, the arguments of numpy.dot can be changed when they
> shouldn't. The following illustrates the problem:
>
> import numpy as npy
> # With sage.rings.integer.Integer:
> x = npy.array([1])
> print x
>   # [1]
> npy.dot(x, npy.array([0]))
> print x
>   # [0]
>
> # This also happens with QQ arrays, but not
> # with any of the finite fields I tried.
> y = MatrixSpace(QQ,2).identity_matrix().numpy()
> print y # [[1 0] [0 1]]
> npy.dot(y, npy.array(VectorSpace(QQ,2).basis()[0]))
> print y # [[1 0] [0 0]]
>
> I am using SAGE version 2.8.7.
>
> Explanations and workarounds are appreciated.

Thank you for reporting this problem!  I've opened a ticket here:
http://sagetrac.org/sage_trac/ticket/1038

Here's the (slightly technical) explanation.  Recently, Robert
Bradshaw added a significant optimization to arithmetic in Sage.  The
optimization is that if you're about to compute a*b, and you know that
the old value of a will no longer be used after this multiplication,
you can reuse the old storage to hold the result of the
multiplication.  This significantly cuts down on the number of memory
allocations and deallocations to perform basic arithmetic.  However,
for this to be safe, we need added restrictions on Python's C API,
above and beyond what Python normally requires.  All of the core Sage
code does follow these restrictions (because it's all generated using
Cython, and Robert added a command-line argument to Cython to generate
compliant code).  However, numpy's C code, unsurprisingly, does not.

The simplest workaround is not to try to store Sage values in numpy
arrays.  Is there a reason you can't use Sage matrices?  If there is
some matrix functionality that numpy has and Sage doesn't, maybe we
could add it to Sage.

At this point, I can only think of two other workarounds; both are
fairly ugly.  The simplest, very ugly and fairly inefficient
workaround is to replace:
x = npy.dot(a, b)
with:
copies = copy(a), copy(b)
x = npy.dot(a, b)

The extra references to the values of a and b held by "copies" should
keep the optimization from triggering.

The other workaround is to remove the optimization from Sage.  Find
the file sage/structure/element.pyx, and search for blocks of code
that look like this:

            if  (<RefPyObject *>left).ob_refcnt < inplace_threshold:
                return _iadd_c(<ModuleElement>left,
<ModuleElement>right)
            else:
                return _add_c(<ModuleElement>left,
<ModuleElement>right)

Remove (or comment out) the first three lines of each block.  There
are eight such blocks, and you can find them by searching for
"inplace_threshold".

Then run "sage -b".  This will recompile that file.

Carl Witty


--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at 
http://groups.google.com/group/sage-support
URLs: http://sage.math.washington.edu/sage/ and http://sage.scipy.org/sage/
-~----------~----~----~----~------~----~------~--~---

Reply via email to