Author: tlinnet
Date: Fri Aug 1 18:09:35 2014
New Revision: 24909
URL: http://svn.gna.org/viewcvs/relax?rev=24909&view=rev
Log:
Removed all un-used helper functions, and matrix exponential functions.
They are now condensed to the fewest possible functions.
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=24909&r1=24908&r2=24909&view=diff
==============================================================================
--- trunk/lib/dispersion/matrix_exponential.py (original)
+++ trunk/lib/dispersion/matrix_exponential.py Fri Aug 1 18:09:35 2014
@@ -71,43 +71,6 @@
return index
-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):
- """ Method to create the helper index matrix, to help figuring out the
index to store in the data matrix. """
-
- # Extract shapes from data.
- NS, NM, NO, ND, Row, Col = data.shape
-
- # Make array to store index.
- index = zeros([NS, NM, NO, ND, 4], int16)
-
- for si in range(NS):
- for mi in range(NM):
- for oi in range(NO):
- for di in range(ND):
- index[si, mi, oi, di] = si, mi, oi, di
-
- return index
-
-
def data_view_via_striding_array_col(data_array=None):
"""Method to stride through the data matrix, extracting the outer array
with nr of elements as Column length.
@@ -239,7 +202,7 @@
complex_flag = any(iscomplex(A))
# If numpy is under 1.8, there would be a need to do eig(A) per matrix.
- if float(version.version[:3]) < 1.9:
+ if float(version.version[:3]) < 1.8:
# Make array to store results
if NE != None:
if dtype != None:
@@ -343,248 +306,6 @@
return array(eA.real)
-def matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(A, dtype=None):
- """Calculate the exact matrix exponential using the eigenvalue
decomposition approach, for higher dimensional data. This of dimension
[NS][NS][NM][NO][ND][X][X].
-
- Here X is the Row and Column length, of the outer square matrix.
-
- @param A: The square matrix to calculate the matrix exponential of.
- @type A: numpy float array of rank [NE][NS][NM][NO][ND][X][X]
- @param dtype: If provided, forces the calculation to use the data type
specified.
- @type type: data-type, optional
- @return: The matrix exponential. This will have the same
dimensionality as the A matrix.
- @rtype: numpy float array of rank [NE][NS][NM][NO][ND][X][X]
- """
-
- # Extract shape.
- NE, NS, NM, NO, ND, Row, Col = A.shape
-
- # Convert dtype, if specified.
- if dtype != None:
- dtype_mat = A.dtype
-
- # If the dtype is different from the input.
- if dtype_mat != dtype:
- # This needs to be made as a copy.
- A = A.astype(dtype)
-
- # Is the original matrix real?
- complex_flag = any(iscomplex(A))
-
- # 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:
- # 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:
- return array(eA)
-
- # Return only the real part.
- else:
- return array(eA.real)
-
-
-def matrix_exponential_rank_NS_NM_NO_ND_x_x(A, dtype=None):
- """Calculate the exact matrix exponential using the eigenvalue
decomposition approach, for higher dimensional data. This of dimension
[NS][NM][NO][ND][X][X].
-
- Here X is the Row and Column length, of the outer square matrix.
-
- @param A: The square matrix to calculate the matrix exponential of.
- @type A: numpy float array of rank [NS][NM][NO][ND][X][X]
- @param dtype: If provided, forces the calculation to use the data type
specified.
- @type type: data-type, optional
- @return: The matrix exponential. This will have the same
dimensionality as the A matrix.
- @rtype: numpy float array of rank [NS][NM][NO][ND][X][X]
- """
-
- # Extract shape.
- NS, NM, NO, ND, Row, Col = A.shape
-
- # Convert dtype, if specified.
- if dtype != None:
- dtype_mat = A.dtype
-
- # If the dtype is different from the input.
- if dtype_mat != dtype:
- # This needs to be made as a copy.
- A = A.astype(dtype)
-
- # Is the original matrix real?
- complex_flag = any(iscomplex(A))
-
- # 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([NS, NM, NO, ND, Row, Col], dtype)
- else:
- eA = zeros([NS, NM, NO, ND, Row, Col], dtype)
-
- # Get the data view, from the helper function.
- A_view = stride_help_square_matrix_rank_NS_NM_NO_ND_x_x(A)
-
- # Create index view.
- index = create_index_rank_NS_NM_NO_ND_x_x(A)
- index_view = stride_help_array_rank_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.
- si, mi, oi, di = index_i
-
- # Store the result.
- eA[si, mi, oi, di, :] = eA_i
-
- else:
- # 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
[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
[NS][NM][NO][ND][X][X]
-
- # Calculate the exponential of all elements in the input array.
Shape [NS][NM][NO][ND][X]
- # Add one axis, to allow for broadcasting multiplication.
- W_exp = exp(W).reshape(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, ...],
(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:
- return array(eA)
-
- # Return only the real part.
- else:
- return array(eA.real)
-
-
def matrix_exponential_rank_NS_NM_NO_ND_2_2(A, dtype=None):
"""Calculate the exact matrix exponential using the closed form in terms
of the matrix elements., for higher dimensional data. This is of dimension
[NS][NM][NO][ND][2][2].
@@ -620,11 +341,11 @@
eA_mat = zeros([NS, NM, NO, ND, Row, Col], dtype)
# Get the data view, from the helper function.
- A_view = stride_help_square_matrix_rank_NS_NM_NO_ND_x_x(A)
+ A_view = data_view_via_striding_array_row_col(data_array=A)
# The index view could be pre formed in init.
- index = create_index_rank_NS_NM_NO_ND_x_x(A)
- index_view = stride_help_array_rank_NS_NM_NO_ND_x(index)
+ index = create_index(NS=NS, NM=NM, NO=NO, ND=ND)
+ index_view = data_view_via_striding_array_col(data_array=index)
# Zip them together and iterate over them.
for A_i, index_i in zip(A_view, index_view):
@@ -680,129 +401,3 @@
eA_mat[si, mi, oi, di, :] = eA_m
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. """
-
- # Extract shapes from data.
- NS, NM, NO, ND, Col = data.shape
-
- # Calculate how many small matrices.
- Nr_mat = 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_square_matrix_rank_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.
- NS, NM, NO, ND, Row, Col = data.shape
-
- # Calculate how many small matrices.
- Nr_mat = 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_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
r24909 - /trunk/lib/dispersion/matrix_exponential.py