Difference between revisions of "Numpy linalg"

From relax wiki
Jump to navigation Jump to search
Line 112: Line 112:
  
 
x=np.arange(row*col).reshape((row, col))
 
x=np.arange(row*col).reshape((row, col))
 +
print "shape is:", x.shape
 
print x
 
print x
  
Line 117: Line 118:
 
# Shape should end out in tuple.
 
# Shape should end out in tuple.
 
shape = tuple([row, col - window + 1,  window])
 
shape = tuple([row, col - window + 1,  window])
print "shape is:", shape
+
print "stride shape is:", shape
  
 
# Get strides
 
# Get strides
Line 127: Line 128:
 
print "itemsize is", x_str
 
print "itemsize is", x_str
  
# Again, strides should be in tuple
+
# 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])
 
cut_strides = tuple(list(x_str) + [x_itz])
 
print "cut strides are", cut_strides
 
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)
 
y = as_strided(x, shape=shape, strides=cut_strides)
 +
 +
print y
 
</source>
 
</source>
  
Line 139: Line 147:
 
from numpy.lib.stride_tricks import as_strided
 
from numpy.lib.stride_tricks import as_strided
 
import numpy as np
 
import numpy as np
 
+
 
NE, NS, NM, NO, ND, Row, Col = 1, 2, 2, 1, 2, 2, 2
 
NE, NS, NM, NO, ND, Row, Col = 1, 2, 2, 1, 2, 2, 2
 
+
mat = np.arange(1,NE*NS*NM*NO*ND*Row*Col+1).reshape(NE, NS, NM, NO, ND, Row, Col)
+
mat = np.arange(NE*NS*NM*NO*ND*Row*Col).reshape(NE, NS, NM, NO, ND, Row, Col)
 
print "mat is:"
 
print "mat is:"
print mat
 
  
sz = mat.itemsize
+
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):
 +
                    print("ei:%i, si:%i, mi:%i, oi:%i, di:%i"%(ei, si, mi, oi, di))
 +
                    print(mat[ei, si, mi, oi, di])
 +
 +
size = mat.itemsize
 
print "itemsize is:"
 
print "itemsize is:"
print sz
+
print size
  
 +
strides = mat.strides
 
print "strides is"
 
print "strides is"
print mat.strides
+
print strides
 
 
print "height and width are"
 
print h, w
 
 
 
bh,bw = Row,Col
 
 
 
shape = (h/bh, w/bw, bh, bw)
 
print shape
 
 
 
 
</source>
 
</source>

Revision as of 09:58, 21 June 2014

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/

http://stackoverflow.com/questions/14758283/is-there-a-numpy-scipy-dot-product-calculating-only-the-diagonal-entries-of-the

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/8070349/using-numpy-stride-tricks-to-get-non-overlapping-array-blocks

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://stackoverflow.com/questions/8070349/using-numpy-stride-tricks-to-get-non-overlapping-array-blocks

http://wiki.scipy.org/Cookbook/GameOfLifeStrides

http://wiki.scipy.org/Cookbook/SegmentAxis

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
import numpy as np
 
NE, NS, NM, NO, ND, Row, Col = 1, 2, 2, 1, 2, 2, 2
 
mat = np.arange(NE*NS*NM*NO*ND*Row*Col).reshape(NE, NS, NM, NO, ND, Row, Col)
print "mat is:"

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):
                    print("ei:%i, si:%i, mi:%i, oi:%i, di:%i"%(ei, si, mi, oi, di))
                    print(mat[ei, si, mi, oi, di])
 
size = mat.itemsize
print "itemsize is:"
print size

strides = mat.strides
print "strides is"
print strides