[Numpy-discussion] Problem with numpy's array reference assignment?

[Charles R Harris]

On Sun, Oct 6, 2013 at 11:30 AM, Bernhard Spinnler <
Bernhard.Spinnler at gmx.net> wrote:

> I have problems to get a piece of code to work with a new numpy/scipy
> version. The code essentially sets up a matrix Ryy and a vector Rya and
> solves the system of linear equations Ryy*c = Rya for c. Then it checks
> whether the resulting vector c satisfies the equation: Ryy*c must be equal
> to Rya.
>
> While this code runs fine on
> - python-2.7.5.p1, numpy-1.7.0, scipy-0.12.0.p0
> - python-2.7.3, numpy-1.7.1, scipy-0.12.0
>
> it fails on
> - python 2.7.2, numpy 1.9.0.dev-fde3dee, scipy 0.14.0.dev-4938da3
> with error:
>
>         AssertionError:
>         Arrays are not almost equal to 6 decimals
>
>         (mismatch 100.0%)
>          x: array([ 9.+0.j,  8.+0.j,  7.+0.j])
>          y: array([ 7.+0.j,  6.+0.j,  5.+0.j])
>
> The piece of code is:
> ------------------------------------
> import numpy as np
> from numpy.testing import assert_array_almost_equal
>
> lag0 = 5
> N = 3
> D = 2
> corr_ya = np.array([0,1,2,3,4,5,6,7,8,9+0j])
> corr_yy = np.array([0,1,2,3,4,5,4,3,2,1])
>
> Rya = corr_ya[lag0+D:lag0+D-N:-1]
> Ryy = np.zeros((N, N), complex)
> for i in np.arange(N):
>     Ryy[i,:] = corr_yy[lag0-i:lag0-i+N]
>
> c = np.linalg.solve(Ryy, Rya)
>
> # check result
> Ryy_x_c = np.dot(Ryy, c)
> print "Ryy*c =", repr(Ryy_x_c)
> print "Rya =", repr(Rya)
> assert_array_almost_equal(Ryy_x_c, Rya)
> ------------------------------------
>
> I guess that it has something to do with numpy assigning arrays by
> reference and not copying them, since the problem goes away when vector Rya
> is copied before passing it to solve, i.e. doing Rya_copy = copy(Rya) and
> passing Rya_copy to solve.
>
> The problem also does not occur when the initial array corr_ya is an
> integer array, e.g. by deleting the "+0j" from the last element of corr_ya
> above. (Normally, my initial arrays are complex. I just used artificially
> simplified arrays above to show the problem.) I could imagine that this is
> due to numpy copying the integers into a new complex array instead of
> creating a reference that is passed to solve.
>
> Could it be that a bug has been introduced in recent numpy/scipy version?
> Or am I misunderstanding something?
>
>
It is certainly possible, there have been changes in the linalg routines.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20131006/cf12ee34/attachment.html>