Numpy linalg

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

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

Stride stride through outer matrix

#!/usr/bin/env python

###############################################################################
#                                                                             #
# Copyright (C) 2014 Troels E. Linnet                                         #
#                                                                             #
# This file is part of the program relax (http://www.nmr-relax.com).          #
#                                                                             #
# This program is free software: you can redistribute it and/or modify        #
# it under the terms of the GNU General Public License as published by        #
# the Free Software Foundation, either version 3 of the License, or           #
# (at your option) any later version.                                         #
#                                                                             #
# This program is distributed in the hope that it will be useful,             #
# but WITHOUT ANY WARRANTY; without even the implied warranty of              #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               #
# GNU General Public License for more details.                                #
#                                                                             #
# You should have received a copy of the GNU General Public License           #
# along with this program.  If not, see <http://www.gnu.org/licenses/>.       #
#                                                                             #
###############################################################################

# Python module imports.
from numpy.lib.stride_tricks import as_strided
from numpy import arange, array, asarray, int16, zeros, uint8
from numpy.linalg import matrix_power

def main():
    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], int16)
    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):
                    for di in range(ND):
                        power_arr[ei, si, mi, oi, di] = 1+ei+si+mi+mi+di
    
    # 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 between matrices. The distance is now separeted by the number of rows.
    bbm = Row * Col * mat_itz

    # Make a tuble of the strides.
    cut_strides = tuple([bbm, bbr, bbe])
    print "cut_strides", cut_strides 

    # Get the view
    mat_view = as_strided(mat, shape=shape, strides=cut_strides)

    print "mat_view.shape", mat_view.shape
    for mat_i in mat_view:
        print mat_i

    # Now for power
    print "now power"
    nr_of_mat = NE*NS*NM*NO*ND
    shape = tuple([nr_of_mat, 1])
    print "View of matrix has shape:", shape
    power_arr_shape = asarray(power_arr.shape)
    print "power_arr_shape", power_arr_shape
    power_itz = power_arr.itemsize
    print "itemsize is", power_itz 
    power_str = list(power_arr.strides)
    print "strides is", power_str, asarray(power_str)/power_itz

    # bytes_between_elements
    bbe = ND * power_itz

    # bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize. NE * NS * NM * NO
    bbr = 1 * power_itz

    cut_strides = tuple([bbr, bbe])
    print "cut_strides", cut_strides

    print power_arr

    power_view = as_strided(power_arr, shape=shape, strides=cut_strides)
    # Power array
    for power_i in power_view:
        print power_i

    collected_power_view = zeros([NE, NS, NM, NO, ND, Row, Col])

    # Zip them together and iterate
    i = 0
    oi = 0
    mi = 0
    for mat_i, pow_i in zip(mat_view, power_view):
        mat_power_view_i = matrix_power(mat_i, int(pow_i))
        # Get the remainder after modulu division
        di = i % ND
        if di == 0 and i != 0:
            oi = (ND-i) % NO
        i += 1
        print mat_i, pow_i, oi, di
        

                


# Execute main function.
if __name__ == "__main__":
    main()