Author: tlinnet
Date: Fri Jun 27 17:09:42 2014
New Revision: 24347
URL: http://svn.gna.org/viewcvs/relax?rev=24347&view=rev
Log:
Initial try to write up a 2x2 matrix by closed form.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/lib/dispersion/matrix_exponential.py
Modified: branches/disp_spin_speed/lib/dispersion/matrix_exponential.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/matrix_exponential.py?rev=24347&r1=24346&r2=24347&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/matrix_exponential.py
(original)
+++ branches/disp_spin_speed/lib/dispersion/matrix_exponential.py Fri
Jun 27 17:09:42 2014
@@ -23,11 +23,30 @@
"""Module for the calculation of the matrix exponential, for higher
dimensional data."""
# Python module imports.
-from numpy import array, any, einsum, eye, exp, iscomplex, newaxis,
multiply, tile
+from numpy import array, any, dot, einsum, eye, exp, iscomplex, int16,
float64, newaxis, multiply, tile, sqrt, zeros
+from numpy.lib.stride_tricks import as_strided
from numpy.linalg import eig, inv
# relax module imports.
from lib.check_types import is_complex
+
+
+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 matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(A, dtype=None):
@@ -169,4 +188,164 @@
# Return only the real part.
else:
- return array(eA.real)
+ 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].
+
+ 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][2][2]
+ @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][2][2]
+ """
+
+ # Extract shapes from data.
+ 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))
+
+ # Make array to store results
+ if dtype != None:
+ eA_mat = zeros([NS, NM, NO, ND, Row, Col], dtype)
+ else:
+ 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)
+
+ # 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)
+
+ # Zip them together and iterate over them.
+ for A_i, index_i in zip(A_view, index_view):
+ a11 = A_i[0, 0]
+ a12 = A_i[0, 1]
+ a21 = A_i[1, 0]
+ a22 = A_i[1, 1]
+
+ # Discriminant
+ a = 1
+ b = -a11 - a22
+ c = a11 * a22 - a12 * a21
+ dis = b**2 - 4*a*c
+
+ # If dis is positive: two distinct real roots
+ # If dis is zero: one distinct real roots
+ # If dis is negative: two complex roots
+
+ # Eigenvalues lambda_1, lambda_2.
+ l1 = (-b + dis) / (2*a)
+ l2 = (-b - dis) / (2*a)
+
+ # Define M: M = V^-1 * [ [l1 0], [0 l2] ] * V
+ W_m = array([ [l1, 0], [0, l2] ])
+
+ v1_1 = 1
+ v1_2 = (l1 - a11) / a12
+
+ v2_1 = 1
+ v2_2 = (l2 - a11) / a12
+
+ # normalized eigenvector
+ v1_1_norm = v1_1 / (sqrt(v1_1**2 + v1_2**2) )
+ v1_2_norm = v1_2 / (sqrt(v1_1**2 + v1_2**2) )
+ v2_1_norm = v2_1 / (sqrt(v2_1**2 + v2_2**2) )
+ v2_2_norm = v2_2 / (sqrt(v2_1**2 + v2_2**2) )
+
+ #V_m = array([ [v1_norm], [v2_norm] ])
+ V_m = array([ [v1_1_norm, v2_1_norm], [v1_2_norm, v2_2_norm] ])
+
+ V_inv_m = inv(V_m)
+
+ # Calculate the exact exponential.
+ dot_V_W = dot(V_m, W_m)
+ eA_m = dot(dot_V_W, V_inv_m)
+
+ # Save results.
+
+ # Extract index from index_view.
+ si, mi, oi, di = index_i
+
+ # Store the result.
+ eA_mat[si, mi, oi, di, :] = eA_m
+
+ return eA_mat
+
+
+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
r24347 - /branches/disp_spin_speed/lib/dispersion/matrix_exponential.py