On Sun, Oct 6, 2013 at 11:30 AM, Bernhard Spinnler < bernhard.spinn...@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
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
