Open main menu

Changes

Numpy linalg

3,407 bytes added, 13:53, 25 June 2014
</source>
And then=== 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/>. #
# #
###############################################################################
 
# 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+1)*arange(1, NDmi+1).astype(int) # Convert to intpower_arr = power_arr.astype(int)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# Make a stride function mat_view = as_strided(mat, shape=shape, that will make a view which is the chain of the outer Row x Col matrices.strides=cut_strides)
# Shape should end out in tuple.nr_of_mat = NE*NS*NM*NO*NDshape = tuple([nr_of_mat, Row, Col]) print "View of matrix has mat_view.shape:", mat_view.shape for mat_i in mat_view: print mat_i
# Get itemsizeNow for power print "now power" nr_of_mat = NE*NS*NM*NO*ND shape = tuple([nr_of_mat, Length 1]) print "View of one array element in bytes. Depends on dtype. float64=8matrix has shape:", complex128shape power_arr_shape =16asarray(power_arr.shape) print "power_arr_shape", power_arr_shapemat_itz power_itz = matpower_arr.itemsize print "itemsize is", mat_itz power_itz power_str = list(power_arr.strides) print "strides is", power_str, asarray(power_str)/power_itz
# Get strides of original matricbytes_between_elementsmat_str bbe = list(mat.strides)print "strides is", mat_str1 * power_itz
# Now we need to calculate how we move with the bytes between row. The distance in the original matrix, when we are creating a view with the# new matrix. It bytes to next row is easiest to start from the back number of the list of stridesColumns elements multiplied with itemsize.NE * NS * NM * NO bbr = ND * power_itz
# bytes_between_elements cut_strides = tuple([bbe, bbr])bbe = 1*mat_itz print cut_strides
# bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize.bbr = Col * mat_itz print power_arr
power_view = as_strided(power_arr, shape=shape, strides=cut_strides) # bytes betPower array for power_i in power_view: print power_i
collected_power_view = zeros([NE, NS, NM, NO, ND, Row, Col])
a # Zip them together and iterate i = zeros0 for mat_i, pow_i in zip([10mat_view, 2power_view): print mat_i, 2pow_i mat_power_view_i = matrix_power(mat_i, 2]int(pow_i))print a print a# Execute main function.itemsizeprint a.stridesif __name__ == "__main__":print asarray main(a.strides)/a.itemsize
</source>