Author: tlinnet
Date: Fri Aug 1 18:09:19 2014
New Revision: 24900
URL: http://svn.gna.org/viewcvs/relax?rev=24900&view=rev
Log:
Implemented first try to stride through data, when computing the eig() of
higher dimensional data.
Systemtest test_cpmg_synthetic_b14_to_ns3d_cluster survived this
transformation.
The systemtest will go from about 11 seconds to 22 seconds.
Modified:
trunk/lib/dispersion/matrix_exponential.py
Modified: trunk/lib/dispersion/matrix_exponential.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/lib/dispersion/matrix_exponential.py?rev=24900&r1=24899&r2=24900&view=diff
==============================================================================
--- trunk/lib/dispersion/matrix_exponential.py (original)
+++ trunk/lib/dispersion/matrix_exponential.py Fri Aug 1 18:09:19 2014
@@ -23,9 +23,28 @@
"""Module for the calculation of the matrix exponential, for higher
dimensional data."""
# Python module imports.
-from numpy import array, any, dot, einsum, eye, exp, iscomplex, int16,
newaxis, multiply, tile, sqrt, zeros
+from numpy import array, any, dot, einsum, eye, exp, iscomplex, int16,
newaxis, multiply, tile, sqrt, version, zeros
from numpy.lib.stride_tricks import as_strided
from numpy.linalg import eig, inv
+
+
+def create_index_rank_NE_NS_NM_NO_ND_x_x(data):
+ """ Method to create the helper index matrix, to help figuring out the
index to store in the data matrix. """
+
+ # Extract shapes from data.
+ NE, NS, NM, NO, ND, Row, Col = data.shape
+
+ # 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_index_rank_NS_NM_NO_ND_x_x(data):
@@ -74,39 +93,89 @@
# Is the original matrix real?
complex_flag = any(iscomplex(A))
- # The eigenvalue decomposition.
- W, V = eig(A)
-
- # W: The eigenvalues, each repeated according to its multiplicity.
- # The eigenvalues are not necessarily ordered.
- # The resulting array will be always be of complex type. Shape
[NE][NS][NM][NO][ND][X]
- # V: The normalized (unit 'length') eigenvectors, such that the column
v[:,i]
- # is the eigenvector corresponding to the eigenvalue w[i]. Shape
[NE][NS][NM][NO][ND][X][X]
-
- # Calculate the exponential of all elements in the input array. Shape
[NE][NS][NM][NO][ND][X]
- # Add one axis, to allow for broadcasting multiplication.
- W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1)
-
- # Make a eye matrix, with Shape [NE][NS][NM][NO][ND][X][X]
- eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, newaxis,
...], (NE, NS, NM, NO, ND, 1, 1) )
-
- # Transform it to a diagonal matrix, with elements from vector down the
diagonal.
- # Use the dtype, if specified.
- if dtype != None:
- W_exp_diag = multiply(W_exp, eye_mat, dtype=dtype )
+ # If numpy is under 1.8, there would be a need to do eig(A) per matrix.
+ if float(version.version[:3]) < 1.8:
+ # Make array to store results
+ if dtype != None:
+ eA = zeros([NE, NS, NM, NO, ND, Row, Col], dtype)
+ else:
+ eA = zeros([NE, NS, NM, NO, ND, Row, Col], dtype)
+
+ # Get the data view, from the helper function.
+ A_view = stride_help_square_matrix_rank_NE_NS_NM_NO_ND_x_x(A)
+
+ # Create index view.
+ index = create_index_rank_NE_NS_NM_NO_ND_x_x(A)
+ index_view = stride_help_array_rank_NE_NS_NM_NO_ND_x(index)
+
+ # Zip them together and iterate over them.
+ for A_i, index_i in zip(A_view, index_view):
+ # The eigenvalue decomposition.
+ W_i, V_i = eig(A_i)
+
+ # Calculate the exponential.
+ W_exp_i = exp(W_i)
+
+ # Make a eye matrix.
+ eye_mat_i = eye(Row)
+
+ # Transform it to a diagonal matrix, with elements from vector
down the diagonal.
+ # Use the dtype, if specified.
+ if dtype != None:
+ W_exp_diag_i = multiply(W_exp_i, eye_mat_i, dtype=dtype )
+ else:
+ W_exp_diag_i = multiply(W_exp_i, eye_mat_i)
+
+ # Make dot product.
+ dot_V_W_i = dot( V_i, W_exp_diag_i)
+
+ # Compute the (multiplicative) inverse of a matrix.
+ inv_V_i = inv(V_i)
+
+ # Calculate the exact exponential.
+ eA_i = dot(dot_V_W_i, inv_V_i)
+
+ # Save results.
+ # Extract index from index_view.
+ ei, si, mi, oi, di = index_i
+
+ # Store the result.
+ eA[ei, si, mi, oi, di, :] = eA_i
+
else:
- W_exp_diag = multiply(W_exp, eye_mat)
-
- # Make dot products for higher dimension.
- # "...", the Ellipsis notation, is designed to mean to insert as many
full slices (:)
- # to extend the multi-dimensional slice to all dimensions.
- dot_V_W = einsum('...ij, ...jk', V, W_exp_diag)
-
- # Compute the (multiplicative) inverse of a matrix.
- inv_V = inv(V)
-
- # Calculate the exact exponential.
- eA = einsum('...ij, ...jk', dot_V_W, inv_V)
+ # The eigenvalue decomposition.
+ W, V = eig(A)
+
+ # W: The eigenvalues, each repeated according to its multiplicity.
+ # The eigenvalues are not necessarily ordered.
+ # The resulting array will be always be of complex type. Shape
[NE][NS][NM][NO][ND][X]
+ # V: The normalized (unit 'length') eigenvectors, such that the
column v[:,i]
+ # is the eigenvector corresponding to the eigenvalue w[i]. Shape
[NE][NS][NM][NO][ND][X][X]
+
+ # Calculate the exponential of all elements in the input array.
Shape [NE][NS][NM][NO][ND][X]
+ # Add one axis, to allow for broadcasting multiplication.
+ W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1)
+
+ # Make a eye matrix, with Shape [NE][NS][NM][NO][ND][X][X]
+ eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, newaxis,
...], (NE, NS, NM, NO, ND, 1, 1) )
+
+ # Transform it to a diagonal matrix, with elements from vector down
the diagonal.
+ # Use the dtype, if specified.
+ if dtype != None:
+ W_exp_diag = multiply(W_exp, eye_mat, dtype=dtype )
+ else:
+ W_exp_diag = multiply(W_exp, eye_mat)
+
+ # Make dot products for higher dimension.
+ # "...", the Ellipsis notation, is designed to mean to insert as
many full slices (:)
+ # to extend the multi-dimensional slice to all dimensions.
+ dot_V_W = einsum('...ij, ...jk', V, W_exp_diag)
+
+ # Compute the (multiplicative) inverse of a matrix.
+ inv_V = inv(V)
+
+ # Calculate the exact exponential.
+ eA = einsum('...ij, ...jk', dot_V_W, inv_V)
# Return the complex matrix.
if complex_flag:
@@ -285,6 +354,36 @@
return eA_mat
+def stride_help_array_rank_NE_NS_NM_NO_ND_x(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_array_rank_NS_NM_NO_ND_x(data):
""" Method to stride through the data matrix, extracting the outer array
with nr of elements as Column length. """
@@ -346,3 +445,36 @@
data_view = as_strided(data, shape=shape, strides=strides)
return data_view
+
+
+def stride_help_square_matrix_rank_NE_NS_NM_NO_ND_x_x(data):
+ """ Helper function calculate the strides to go 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
r24900 - /trunk/lib/dispersion/matrix_exponential.py