Explanation of issue with multiplying matrices element-wise and numpy ndarray views. This is what caused to the difference between the numpy and fortran pulse propagator functions in my blog post (http://themodernscientist.com/posts/2013/2013-06-09-simulation_of_nmr_shaped_pulses/)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
a = np.eye(3) | |
b = np.eye(3) | |
c = np.eye(3) | |
for i in range(10): | |
d = np.random.random((3,3)) | |
# copy each matrix element individually | |
d00 = d[0,0].copy() | |
d01 = d[0,1].copy() | |
d02 = d[0,2].copy() | |
d10 = d[1,0].copy() | |
d11 = d[1,1].copy() | |
d12 = d[1,2].copy() | |
d20 = d[2,0].copy() | |
d21 = d[2,1].copy() | |
d22 = d[2,2].copy() | |
c00 = c[0,0].copy() | |
c01 = c[0,1].copy() | |
c02 = c[0,2].copy() | |
c10 = c[1,0].copy() | |
c11 = c[1,1].copy() | |
c12 = c[1,2].copy() | |
c20 = c[2,0].copy() | |
c21 = c[2,1].copy() | |
c22 = c[2,2].copy() | |
# a = d * a | |
a = np.dot(d,a) | |
# b = d * b --> this operates on variables that are views | |
# which causes their values to change during the multiplication | |
b[0,0] = d[0,0]*b[0,0] + d[0,1]*b[1,0] + d[0,2]*b[2,0] | |
b[0,1] = d[0,0]*b[0,1] + d[0,1]*b[1,1] + d[0,2]*b[2,1] | |
b[0,2] = d[0,0]*b[0,2] + d[0,1]*b[1,2] + d[0,2]*b[2,2] | |
b[1,0] = d[1,0]*b[0,0] + d[1,1]*b[1,0] + d[1,2]*b[2,0] | |
b[1,1] = d[1,0]*b[0,1] + d[1,1]*b[1,1] + d[1,2]*b[2,1] | |
b[1,2] = d[1,0]*b[0,2] + d[1,1]*b[1,2] + d[1,2]*b[2,2] | |
b[2,0] = d[2,0]*b[0,0] + d[2,1]*b[1,0] + d[2,2]*b[2,0] | |
b[2,1] = d[2,0]*b[0,1] + d[2,1]*b[1,1] + d[2,2]*b[2,1] | |
b[2,2] = d[2,0]*b[0,2] + d[2,1]*b[1,2] + d[2,2]*b[2,2] | |
# c = d * c | |
c[0,0] = d00*c00 + d01*c10 + d02*c20 | |
c[0,1] = d00*c01 + d01*c11 + d02*c21 | |
c[0,2] = d00*c02 + d01*c12 + d02*c22 | |
c[1,0] = d10*c00 + d11*c10 + d12*c20 | |
c[1,1] = d10*c01 + d11*c11 + d12*c21 | |
c[1,2] = d10*c02 + d11*c12 + d12*c22 | |
c[2,0] = d20*c00 + d21*c10 + d22*c20 | |
c[2,1] = d20*c01 + d21*c11 + d22*c21 | |
c[2,2] = d20*c02 + d21*c12 + d22*c22 | |
print 'max(abs(a-b)) ', np.max(np.abs(a-b)) | |
print 'max(abs(a-c)) ', np.max(np.abs(a-c)) | |
print 'a', a | |
print 'b', b #oops | |
print 'c', c |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Here are the results: