Author: tlinnet
Date: Wed Jun 25 20:41:47 2014
New Revision: 24325
URL: http://svn.gna.org/viewcvs/relax?rev=24325&view=rev
Log:
First try to implement function that will calculate the matrix exponential
by striding through data.
Interestingly, it does not work. Theses systemtests will fail.
test_hansen_cpmg_data_to_ns_cpmg_2site_3D
test_hansen_cpmg_data_to_ns_cpmg_2site_3D_full
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Added:
branches/disp_spin_speed/lib/dispersion/matrix_power.py
Added: branches/disp_spin_speed/lib/dispersion/matrix_power.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/matrix_power.py?rev=24325&view=auto
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/matrix_power.py (added)
+++ branches/disp_spin_speed/lib/dispersion/matrix_power.py Wed Jun 25
20:41:47 2014
@@ -0,0 +1,184 @@
+###############################################################################
+#
#
+# 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/>.
#
+#
#
+###############################################################################
+
+# Module docstring.
+"""Module for the calculation of the matrix power, for higher dimensional
data."""
+
+# Python module imports.
+from numpy.lib.stride_tricks import as_strided
+from numpy import arange, array, asarray, float64, int16, sum, zeros
+from numpy.linalg import matrix_power
+
+
+def create_index(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 matrix_power_strided_rank_NE_NS_NM_NO_ND_x_x(data, power):
+ """Calculate the exact matrix power by striding through higher
dimensional data. This of dimension [NE][NS][NM][NO][ND][X][X].
+
+ Here X is the Row and Column length, of the outer square matrix.
+
+ @param data: The square matrix to calculate the matrix
exponential of.
+ @type data: numpy float array of rank
[NE][NS][NM][NO][ND][X][X]
+ @keyword power: The matrix exponential power array, which hold the
power integer to which to raise the outer matrix X,X to.
+ @type power: numpy int array of rank [NE][NS][NM][NO][ND]
+ @return: The matrix pwer. This will have the same
dimensionality as the data matrix.
+ @rtype: numpy float array of rank
[NE][NS][NM][NO][ND][X][X]
+ """
+
+ # 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], float64)
+
+ # Get the data view, from the helper function.
+ data_view = stride_help_square_matrix_rank_NE_NS_NM_NO_ND_x_x(data)
+
+ # Get the power view, from the helper function.
+ power_view = stride_help_element_rank_NE_NS_NM_NO_ND(power)
+
+ # The index view could be pre formed in init.
+ index = create_index(data)
+ index_view = stride_help_array_rank_NE_NS_NM_NO_ND_x(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 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_element_rank_NE_NS_NM_NO_ND(data):
+ """ Method to stride through the data matrix, extracting the outer
element. """
+
+ # Extract shapes from data.
+ NE, NS, NM, NO, Col = data.shape
+
+ # Calculate how many small matrices.
+ Nr_mat = NE * NS * NM * NO * Col
+
+ # Define the shape for the stride view.
+ shape = (Nr_mat, 1)
+
+ # Get itemsize, Length of one array element in bytes. Depends on
dtype. float64=8, complex128=16.
+ itz = data.itemsize
+
+ # FIXME: Explain this.
+ bbe = Col * itz
+
+ # FIXME: Explain this.
+ bbr = 1 * 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_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
+
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
relax-commits@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits
Re: r24325 - /branches/disp_spin_speed/lib/dispersion/matrix_power.py