Numpy linalg
Jump to navigation
Jump to search
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
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/>. #
# #
###############################################################################
"""
This script is for testing how to stride over matrices, in a large data array, when they are positioned in the end of the dimensions.
It also serves as a validation tool, and an efficient way to profile the calculation.
It is optimal to eventually try to implement a faster matrix power calculation.
"""
# Python module imports.
import cProfile
import pstats
import tempfile
from numpy.lib.stride_tricks import as_strided
from numpy import arange, array, asarray, int16, sum, zeros
from numpy.linalg import matrix_power
def main():
# Nr of iterations.
nr_iter = 50
# Print statistics.
verbose = True
# Calc for single_normal.
if True:
#if False:
# Define filename
sn_filename = tempfile.NamedTemporaryFile(delete=False).name
# Profile for a single spin.
cProfile.run('single_normal(iter=%s)'%nr_iter, sn_filename)
# Read all stats files into a single object
sn_stats = pstats.Stats(sn_filename)
# Clean up filenames for the report
sn_stats.strip_dirs()
# Sort the statistics by the cumulative time spent in the function. cumulative, time, calls
sn_stats.sort_stats('cumulative')
# Print report for single.
if verbose:
print("This is the report for single normal.")
sn_stats.print_stats()
# Calc for single_strided.
if True:
#if False:
# Define filename
ss_filename = tempfile.NamedTemporaryFile(delete=False).name
# Profile for a single spin.
cProfile.run('single_strided(iter=%s)'%nr_iter, ss_filename)
# Read all stats files into a single object
ss_stats = pstats.Stats(ss_filename)
# Clean up filenames for the report
ss_stats.strip_dirs()
# Sort the statistics by the cumulative time spent in the function. cumulative, time, calls
ss_stats.sort_stats('cumulative')
# Print report for single.
if verbose:
print("This is the report for single strided.")
ss_stats.print_stats()
# Calc for cluster_normal.
if True:
#if False:
# Define filename
cn_filename = tempfile.NamedTemporaryFile(delete=False).name
# Profile for a cluster spin.
cProfile.run('cluster_normal(iter=%s)'%nr_iter, cn_filename)
# Read all stats files into a single object
cn_stats = pstats.Stats(cn_filename)
# Clean up filenames for the report
cn_stats.strip_dirs()
# Sort the statistics by the cumulative time spent in the function. cumulative, time, calls
cn_stats.sort_stats('cumulative')
# Print report for cluster.
if verbose:
print("This is the report for cluster normal.")
cn_stats.print_stats()
# Calc for cluster_strided.
if True:
#if False:
# Define filename
cs_filename = tempfile.NamedTemporaryFile(delete=False).name
# Profile for a cluster spin.
cProfile.run('cluster_strided(iter=%s)'%nr_iter, cs_filename)
# Read all stats files into a single object
cs_stats = pstats.Stats(cs_filename)
# Clean up filenames for the report
cs_stats.strip_dirs()
# Sort the statistics by the cumulative time spent in the function. cumulative, time, calls
cs_stats.sort_stats('cumulative')
# Print report for cluster.
if verbose:
print("This is the report for cluster strided.")
cs_stats.print_stats()
class Profile():
"""
Class Profile inherits the Dispersion container class object.
"""
def __init__(self, NE=1, NS=3, NM=2, NO=1, ND=20, Row=7, Col=7):
# Setup the size of data array:
self.NE = NE
self.NS = NS
self.NM = NM
self.NO = NO
self.ND = ND
self.Row = Row
self.Col = Col
# Create the data matrix
self.data = self.create_data()
# Create the index matrix
self.index = self.create_index()
# Create the power matrix
self.power = self.create_power()
# Create the validation matrix
self.vali = self.create_vali()
def create_data(self):
""" Method to create the imagninary data structure"""
NE = self.NE
NS = self.NS
NM = self.NM
NO = self.NO
ND = self.ND
Row = self.Row
Col = self.Col
# Create the data matrix for testing.
data = arange(NE*NS*NM*NO*ND*Row*Col).reshape(NE, NS, NM, NO, ND, Row, Col)
return data
def create_index(self):
""" Method to create the helper index matrix, to help figuring out the index to store in the data matrix. """
NE = self.NE
NS = self.NS
NM = self.NM
NO = self.NO
ND = self.ND
# Make array to store index.
index = zeros([NE, NS, NM, NO, ND, 5], int16)
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):
index[ei, si, mi, oi, di] = ei, si, mi, oi, di
return index
def create_power(self):
""" Method to create the test power data array. """
NE = self.NE
NS = self.NS
NM = self.NM
NO = self.NO
ND = self.ND
power = zeros([NE, NS, NM, NO, ND], int16)
# 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[ei, si, mi, oi, di] = 1+ei+si+mi+mi+di
return power
def create_vali(self):
""" Create validation matrix for power array calculation. """
NE = self.NE
NS = self.NS
NM = self.NM
NO = self.NO
ND = self.ND
Row = self.Row
Col = self.Col
power = self.power
data = self.data
# Make array to store results
vali = 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):
# Get the outer square data matrix,
data_i = data[ei, si, mi, oi, di]
# Get the power.
power_i = power[ei, si, mi, oi, di]
# Do power calculation with numpy method.
data_power_i = matrix_power(data_i, power_i)
# Store result.
vali[ei, si, mi, oi, di] = data_power_i
return vali
def stride_help_square_matrix(self, data):
""" Method to stride through the data matrix, extracting the outer Row X Col matrix. """
# Extract shapes from data.
NE, NS, NM, NO, ND, Row, Col = data.shape
# Calculate how many small matrices.
Nr_mat = NE * NS * NM * NO * ND
# Define the shape for the stride view.
shape = (Nr_mat, Row, Col)
# Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16.
itz = data.itemsize
# Bytes_between_elements
bbe = 1 * itz
# Bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize.
bbr = Col * itz
# Bytes between matrices. The byte distance is separated by the number of rows.
bbm = Row * Col * itz
# Make a tuple of the strides.
strides = (bbm, bbr, bbe)
# Make the stride view.
data_view = as_strided(data, shape=shape, strides=strides)
return data_view
def stride_help_array(self, data):
""" Method to stride through the data matrix, extracting the outer array with nr of elements as Column length. """
# Extract shapes from data.
NE, NS, NM, NO, ND, Col = data.shape
# Calculate how many small matrices.
Nr_mat = NE * NS * NM * NO * ND
# Define the shape for the stride view.
shape = (Nr_mat, Col)
# Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16.
itz = data.itemsize
# Bytes_between_elements
bbe = 1 * itz
# Bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize.
bbr = Col * itz
# Make a tuple of the strides.
strides = (bbr, bbe)
# Make the stride view.
data_view = as_strided(data, shape=shape, strides=strides)
return data_view
def stride_help_element(self, data):
""" Method to stride through the data matrix, extracting the outer element. """
# Extract shapes from data.
NE, NS, NM, NO, Col = data.shape
# Calculate how many small matrices.
Nr_mat = NE * NS * NM * NO * Col
# Define the shape for the stride view.
shape = (Nr_mat, 1)
# Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16.
itz = data.itemsize
# FIXME: Explain this.
bbe = Col * itz
# FIXME: Explain this.
bbr = 1 * itz
# Make a tuple of the strides.
strides = (bbr, bbe)
# Make the stride view.
data_view = as_strided(data, shape=shape, strides=strides)
return data_view
def calc_normal(self, data, power):
""" The normal way to do the calculation """
# Extract shapes from data.
NE, NS, NM, NO, ND, Row, Col = data.shape
# Make array to store results
calc = 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):
# Get the outer square data matrix,
data_i = data[ei, si, mi, oi, di]
# Get the power.
power_i = power[ei, si, mi, oi, di]
# Do power calculation with numpy method.
data_power_i = matrix_power(data_i, power_i)
# Store result.
calc[ei, si, mi, oi, di] = data_power_i
return calc
def calc_strided(self, data, power):
""" The strided way to do the calculation """
# Extract shapes from data.
NE, NS, NM, NO, ND, Row, Col = data.shape
# Make array to store results
calc = zeros([NE, NS, NM, NO, ND, Row, Col])
# Get the data view, from the helper function.
data_view = self.stride_help_square_matrix(data)
# Get the power view, from the helpwer function.
power_view = self.stride_help_element(power)
# The index view could be pre formed in init.
index = self.index
index_view = self.stride_help_array(index)
# Zip them together and iterate over them.
for data_i, power_i, index_i in zip(data_view, power_view, index_view):
# Do power calculation with numpy method.
data_power_i = matrix_power(data_i, int(power_i))
# Extract index from index_view.
ei, si, mi, oi, di = index_i
# Store the result.
calc[ei, si, mi, oi, di] = data_power_i
return calc
def validate(self):
""" The validation of calculations """
data = self.data
power = self.power
# Calculate by normal way.
calc_normal = self.calc_normal(data, power)
# Find the difference to the validated method.
diff_normal = calc_normal - self.vali
if sum(diff_normal) != 0.0:
print("The normal method is different from the validated data")
# Calculate by strided way.
calc_strided = self.calc_strided(data, power)
# Find the difference to the validated method.
diff_strided = calc_strided - self.vali
if sum(diff_strided) != 0.0:
print("The strided method is different from the validated data")
def single_normal(NS=1, iter=None):
# Instantiate class
MC = Profile(NS=NS)
# Get the init data.
data = MC.data
power = MC.power
# Loop 100 times for each spin in the clustered analysis (to make the timing numbers equivalent).
for spin_index in xrange(100):
# Repeat the function call, to simulate minimisation.
for i in xrange(iter):
calc = MC.calc_normal(data, power)
def single_strided(NS=1, iter=None):
# Instantiate class
MC = Profile(NS=NS)
# Get the init data.
data = MC.data
power = MC.power
# Loop 100 times for each spin in the clustered analysis (to make the timing numbers equivalent).
for spin_index in xrange(100):
# Repeat the function call, to simulate minimisation.
for i in xrange(iter):
calc = MC.calc_strided(data, power)
def cluster_normal(NS=100, iter=None):
# Instantiate class
MC = Profile(NS=NS)
# Get the init data.
data = MC.data
power = MC.power
# Repeat the function call, to simulate minimisation.
for i in xrange(iter):
calc = MC.calc_normal(data, power)
def cluster_strided(NS=100, iter=None):
# Instantiate class
MC = Profile(NS=NS)
# Get the init data.
data = MC.data
power = MC.power
# Repeat the function call, to simulate minimisation.
for i in xrange(iter):
calc = MC.calc_strided(data, power)
# First validate
#Initiate My Class.
MC = Profile()
# Validate all calculations.
MC.validate()
# Execute main function.
if __name__ == "__main__":
# Initiate cProfiling.
main()
