Open main menu

Changes

Numpy linalg

9,921 bytes added, 17:26, 25 June 2014
# #
###############################################################################
 
"""
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, uint8
from numpy.linalg import matrix_power
def main():
NE# 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, NStime, NMcalls 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, NOcs_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, NDtime, Rowcalls 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, Col NE= 1, NS=3, NM=2, NO=1, 4ND=20, 2Row=7, 2Col=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
mat # 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   power_arr def create_index(self): """ Method to create the helper index matrix, to help figuring out the index to store in the data matrix. """  NE = zeros([self.NE, NS = self.NS, NM = self.NM, NO = self.NO, ND = self.ND], int16) print "mat is:"
# 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_arr 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)
# Make Get itemsize, Length of one array to store resultselement in bytes. Depends on dtype. float64=8, complex128=16. collected_power itz = zeros([NE, NS, NM, NO, ND, Row, Col])data.itemsize
# Normal way, by loop of loops. for ei in range(NE):Bytes_between_elements for si in range(NS): for mi in range(NM): for oi in range(NO): for di in range(ND): power_i bbe = 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_i1 * itz
mat_power_i = matrix_power(mat_i, power_i) print mat_power_i # Store powerBytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize. collected_power[ei, si, mi, oi, di] bbr = mat_power_iCol * 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   printdef stride_help_array(self, data): """ Method to stride through the data matrix, extracting the outer array with nr of elements as Column length. ""Collected power shape is: %s") print collected_power # Extract shapes from data.shape print collected_power[ei NE, NS, NM, siNO, miND, oiCol = data.shape  # Calculate how many small matrices. Nr_mat = NE * NS * NM * NO * ND  # Define the shape for the stride view. shape = (Nr_mat, di]Col)
# Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16.
itz = data.itemsize
## #Bytes_between_elements # Make a stride function, that will make a view which is the chain of the outer Row x Col matrices. bbe = 1 * itz
# Shape should end out Bytes between row. The distance in tuplebytes to next row is number of Columns elements multiplied with itemsize. nr_of_mat bbr = NE*NSCol *NM*NO*ND shape = tuple([nr_of_mat, Row, Col]) print "View of matrix has shape:", shapeitz
# Get itemsize, Length Make a tuple of one array element in bytesthe strides. Depends on dtype strides = (bbr, bbe)  # Make the stride view. float64 data_view = as_strided(data, shape=8shape, complex128strides=16.strides)  return data_view   mat_itz = matdef stride_help_element(self, data): """ Method to stride through the data matrix, extracting the outer element.itemsize print "itemsize is""  # Extract shapes from data. NE, mat_itz NS, NM, NO, Col = data.shape # Get strides of original matricCalculate how many small matrices. mat_str Nr_mat = list(matNE * NS * NM * NO * Col  # Define the shape for the stride view.strides) print "strides is" shape = (Nr_mat, mat_str1)
# Now we need to calculate how we move with the Get itemsize, Length of one array element in bytes in the original matrix. Depends on dtype. float64=8, when we are creating a view with thecomplex128=16. # new matrix. It is easiest to start from the back of the list of strides itz = data.itemsize
# bytes_between_elementsFIXME: Explain this. bbe = 1Col *mat_itzitz
# bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsizeFIXME: Explain this. bbr = Col 1 * mat_itzitz
# bytes between matrices. The distance is now separeted by Make a tuple of the number of rowsstrides. bbm strides = Row * Col * mat_itz(bbr, bbe)
# Make a tuble of the stridesstride view. cut_strides data_view = tupleas_strided([bbmdata, bbrshape=shape, bbe]strides=strides) print "cut_strides", cut_strides
# Get the view mat_view = as_strided(mat, shape=shape, strides=cut_strides) return data_view
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 = tupledef calc_normal([nr_of_matself, data, 1]power) 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 isThe normal way to do the calculation ", power_itz power_str = list(power_arr.strides) print "strides is", power_str, asarray(power_str)/power_itz
# bytes_between_elementsExtract shapes from data. bbe NE, NS, NM, NO, ND, Row, Col = ND * power_itzdata.shape
# bytes between row. The distance in bytes Make array to next row is number of Columns elements multiplied with itemsize. store results calc = zeros([NE * , NS * , NM * , NO bbr = 1 * power_itz, ND, Row, Col])
cut_strides # 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 = tuple(data[bbrei, bbesi, mi, oi, di]) print "cut_strides", cut_strides
print power_arr # Get the power. power_i = power[ei, si, mi, oi, di]
power_view # Do power calculation with numpy method. data_power_i = as_stridedmatrix_power(power_arrdata_i, shape=shape, strides=cut_stridespower_i) # Power array for power_i in power_view: print power_i
collected_power_view = zeros( # Store result. calc[NE, NS, NMei, NOsi, NDmi, Rowoi, Coldi])= data_power_i
# 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, direturn 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>