View source for Numpy linalg

__TOC__

== 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

<source lang="python">
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)
</source>

=== 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.

<source lang="python">
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

</source>


== 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
<source lang="python">
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
</source>

=== Stride stride through outer matrix ===

<source lang="python">
#!/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()
</source>

== See also ==

[[Category:Development]]