<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Sun, Oct 6, 2013 at 11:30 AM, Bernhard Spinnler <span dir="ltr"><<a href="mailto:Bernhard.Spinnler@gmx.net" target="_blank">Bernhard.Spinnler@gmx.net</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">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.<br>

<br>
While this code runs fine on<br>
- python-2.7.5.p1, numpy-1.7.0, scipy-0.12.0.p0<br>
- python-2.7.3, numpy-1.7.1, scipy-0.12.0<br>
<br>
it fails on<br>
- python 2.7.2, numpy 1.9.0.dev-fde3dee, scipy 0.14.0.dev-4938da3<br>
with error:<br>
<br>
        AssertionError:<br>
        Arrays are not almost equal to 6 decimals<br>
<br>
        (mismatch 100.0%)<br>
         x: array([ 9.+0.j,  8.+0.j,  7.+0.j])<br>
         y: array([ 7.+0.j,  6.+0.j,  5.+0.j])<br>
<br>
The piece of code is:<br>
------------------------------------<br>
import numpy as np<br>
from numpy.testing import assert_array_almost_equal<br>
<br>
lag0 = 5<br>
N = 3<br>
D = 2<br>
corr_ya = np.array([0,1,2,3,4,5,6,7,8,9+0j])<br>
corr_yy = np.array([0,1,2,3,4,5,4,3,2,1])<br>
<br>
Rya = corr_ya[lag0+D:lag0+D-N:-1]<br>
Ryy = np.zeros((N, N), complex)<br>
for i in np.arange(N):<br>
    Ryy[i,:] = corr_yy[lag0-i:lag0-i+N]<br>
<br>
c = np.linalg.solve(Ryy, Rya)<br>
<br>
# check result<br>
Ryy_x_c = np.dot(Ryy, c)<br>
print "Ryy*c =", repr(Ryy_x_c)<br>
print "Rya =", repr(Rya)<br>
assert_array_almost_equal(Ryy_x_c, Rya)<br>
------------------------------------<br>
<br>
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.<br>

<br>
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.<br>

<br>
Could it be that a bug has been introduced in recent numpy/scipy version? Or am I misunderstanding something?<br>
<br></blockquote><div><br></div><div>It is certainly possible, there have been changes in the linalg routines.<br><br></div><div>Chuck </div></div></div></div>