Numpy linalg
Contents
How to transpose higher dimension arrays
http://jameshensman.wordpress.com/2010/06/14/multiple-matrix-multiplication-in-numpy/
Faster dot product using BLAS
http://www.huyng.com/posts/faster-numpy-dot-product/
http://stackoverflow.com/questions/5990577/speeding-up-numpy-dot
http://wiki.scipy.org/PerformanceTips
http://thread.gmane.org/gmane.comp.python.numeric.general/28135/
Multi dot
http://wiki.scipy.org/Cookbook/MultiDot
Einsum
http://chintaksheth.wordpress.com/2013/07/31/numpy-the-tricks-of-the-trade-part-ii/
a = np.arange(4).reshape(2,2)
print a
print "np.einsum('ii', a), row i multiplied downwards"
print np.einsum('ii', a)
print "np.einsum('ij', a), same matrix ?"
print np.einsum('ij', a)
print "np.einsum('ji', a), transpose"
print np.einsum('ji', a)
print "np.einsum('ij,jk', a, a), dot product"
print np.einsum('ij,jk', a, a)
print np.dot(a, a)
Ellipsis broadcasting in numpy.einsum
http://stackoverflow.com/questions/16591696/ellipsis-broadcasting-in-numpy-einsum
http://comments.gmane.org/gmane.comp.python.numeric.general/53705
http://stackoverflow.com/questions/118370/how-do-you-use-the-ellipsis-slicing-syntax-in-python
http://stackoverflow.com/questions/772124/what-does-the-python-ellipsis-object-do
"..." Is designed to mean at this point, insert as many full slices (:) to extend the multi-dimensional slice to all dimensions.
print ""
a = np.arange(4).reshape(2,2)
print "a is"
print a
print "dot a"
print np.dot(a, a)
# Expand one axis in start, and tile up 2 times.
a2 = np.tile(a[None,:], (2, 1, 1))
print "a2 shape", a2.shape
print "einsum dot product over higher dimensions"
a2_e = np.einsum('...ij,...jk', a2, a2)
print a2_e
# Expand one axis in start, and tile up 2 times.
a3 = np.tile(a2[None,:], (2, 1, 1, 1))
print "a3 shape", a3.shape
print "einsum dot product over higher dimensions"
a3_e = np.einsum('...ij,...jk', a3, a3)
print a3_e
# With Ellipsis and axis notation
a = np.arange(4).reshape(2,2)
a = np.arange(4).reshape(2,2)
print "a is"
print a
print "dot a"
print np.dot(a, a)
# Expand one axis in start, and tile up 2 times.
a2 = np.tile(a[None,:], (2, 1, 1))
print "a2 shape", a2.shape
print "einsum dot product over higher dimensions"
a2_e = np.einsum(a2, [Ellipsis, 0, 1], a2, [Ellipsis, 1, 2])
print a2_e
# Expand one axis in start, and tile up 2 times.
a3 = np.tile(a2[None,:], (2, 1, 1, 1))
print "a3 shape", a3.shape
print "einsum dot product over higher dimensions"
a3_e = np.einsum(a3, [Ellipsis, 0, 1], a3, [Ellipsis, 1, 2])
print a3_e
Stride tricks
http://chintaksheth.wordpress.com/2013/07/31/numpy-the-tricks-of-the-trade-part-ii/
http://stackoverflow.com/questions/4936620/using-strides-for-an-efficient-moving-average-filter
http://www.rigtorp.se/2011/01/01/rolling-statistics-numpy.html
http://wiki.scipy.org/Cookbook/GameOfLifeStrides
http://wiki.scipy.org/Cookbook/SegmentAxis
http://scipy-lectures.github.io/advanced/advanced_numpy/
Simple, with window 3
row, col = 3, 5
x=np.arange(row*col).reshape((row, col))
print "shape is:", x.shape
print x
window = 3
# Shape should end out in tuple.
shape = tuple([row, col - window + 1, window])
print "stride shape is:", shape
# Get strides
x_str = x.strides
print "strides is", x_str
# Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16.
x_itz = x.itemsize
print "itemsize is", x_str
# Again, strides should be in tuple.
# Strides tuble is a way to tell, how far is to next element.
# Next element is first to next row in bytes = col * itemsize
# and then how many bytes to next element in the row.
cut_strides = tuple(list(x_str) + [x_itz])
print "cut strides are", cut_strides
cut_strides_cal = tuple([col*x_itz] + [x_itz] + [x_itz])
print "Can be calculated as [col*x_itz] + [x_itz]:", cut_strides_cal
y = as_strided(x, shape=shape, strides=cut_strides)
print y
And then
from numpy.lib.stride_tricks import as_strided
from numpy import arange, asarray, int16, zeros
from numpy.linalg import matrix_power
NE, NS, NM, NO, ND, Row, Col = 1, 3, 2, 1, 4, 2, 2
mat = arange(NE*NS*NM*NO*ND*Row*Col).reshape(NE, NS, NM, NO, ND, Row, Col)
power_arr = zeros([NE, NS, NM, NO, ND])
print "mat is:"
# Preform the power array, to be outside profiling test.
for ei in range(NE):
for si in range(NS):
for mi in range(NM):
for oi in range(NO):
power_arr[ei, si, mi, oi, :] = (mi+1)*arange(1, ND+1).astype(int)
# Convert to int
power_arr = power_arr.astype(int)
# Make array to store results
collected_power = zeros([NE, NS, NM, NO, ND, Row, Col])
# Normal way, by loop of loops.
for ei in range(NE):
for si in range(NS):
for mi in range(NM):
for oi in range(NO):
for di in range(ND):
power_i = power_arr[ei, si, mi, oi, di]
print("ei:%i, si:%i, mi:%i, oi:%i, di:%i, power:%i"%(ei, si, mi, oi, di, power_i))
mat_i = mat[ei, si, mi, oi, di]
#print mat_i
mat_power_i = matrix_power(mat_i, power_i)
print mat_power_i
# Store power
collected_power[ei, si, mi, oi, di] = mat_power_i
print("Collected power shape is: %s")
print collected_power.shape
print collected_power[ei, si, mi, oi, di]
###
# Make a stride function, that will make a view which is the chain of the outer Row x Col matrices.
# Shape should end out in tuple.
nr_of_mat = NE*NS*NM*NO*ND
shape = tuple([nr_of_mat, Row, Col])
print "View of matrix has shape:", shape
# Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16.
mat_itz = mat.itemsize
print "itemsize is", mat_itz
# Get strides of original matric
mat_str = list(mat.strides)
print "strides is", mat_str
# Now we need to calculate how we move with the bytes in the original matrix, when we are creating a view with the
# new matrix. It is easiest to start from the back of the list of strides.
# bytes_between_elements
bbe = 1*mat_itz
# bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize.
bbr = Col * mat_itz
# bytes bet
a = zeros([10, 2, 2, 2])
print a
print a.itemsize
print a.strides
print asarray(a.strides)/a.itemsize
