Source Code for Module lib.dispersion.matrix_power

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2014 Troels E. Linnet                                         # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """Module for the calculation of the matrix power, for higher dimensional data.""" 
 24   
 25  # Python module imports. 
 26  from numpy.lib.stride_tricks import as_strided 
 27  from numpy import float64, int16, zeros 
 28  from numpy.linalg import matrix_power 
 29   
 30   
 31 -def create_index(data): 
 32      """ Method to create the helper index matrix, to help figuring out the index to store in the data matrix. """ 
 33   
 34      # Extract shapes from data. 
 35      NE, NS, NM, NO, ND, Row, Col = data.shape 
 36   
 37      # Make array to store index. 
 38      index = zeros([NE, NS, NM, NO, ND, 5], int16) 
 39   
 40      for ei in range(NE): 
 41          for si in range(NS): 
 42              for mi in range(NM): 
 43                  for oi in range(NO): 
 44                      for di in range(ND): 
 45                          index[ei, si, mi, oi, di] = ei, si, mi, oi, di 
 46   
 47      return index 
 48   
 49   
 50 -def matrix_power_strided_rank_NE_NS_NM_NO_ND_x_x(data, power): 
 51      """Calculate the exact matrix power by striding through higher dimensional data.  This of dimension [NE][NS][NM][NO][ND][X][X]. 
 52   
 53      Here X is the Row and Column length, of the outer square matrix. 
 54   
 55      @param data:        The square matrix to calculate the matrix exponential of. 
 56      @type data:         numpy float array of rank [NE][NS][NM][NO][ND][X][X] 
 57      @keyword power:     The matrix exponential power array, which hold the power integer to which to raise the outer matrix X,X to. 
 58      @type power:        numpy int array of rank [NE][NS][NM][NO][ND] 
 59      @return:            The matrix pwer.  This will have the same dimensionality as the data matrix. 
 60      @rtype:             numpy float array of rank [NE][NS][NM][NO][ND][X][X] 
 61      """ 
 62   
 63      # Extract shapes from data. 
 64      NE, NS, NM, NO, ND, Row, Col = data.shape 
 65   
 66      # Make array to store results 
 67      calc = zeros([NE, NS, NM, NO, ND, Row, Col], float64) 
 68   
 69      # Get the data view, from the helper function. 
 70      data_view = stride_help_square_matrix_rank_NE_NS_NM_NO_ND_x_x(data) 
 71   
 72      # Get the power view, from the helper function. 
 73      power_view = stride_help_element_rank_NE_NS_NM_NO_ND(power) 
 74   
 75      # The index view could be pre formed in init. 
 76      index = create_index(data) 
 77      index_view = stride_help_array_rank_NE_NS_NM_NO_ND_x(index) 
 78   
 79      # Zip them together and iterate over them. 
 80      for data_i, power_i, index_i in zip(data_view, power_view, index_view): 
 81          # Do power calculation with numpy method. 
 82          data_power_i = matrix_power(data_i, int(power_i)) 
 83   
 84          # Extract index from index_view. 
 85          ei, si, mi, oi, di = index_i 
 86   
 87          # Store the result. 
 88          calc[ei, si, mi, oi, di, :] = data_power_i 
 89   
 90      return calc 
 91   
 92   
 93 -def stride_help_array_rank_NE_NS_NM_NO_ND_x(data): 
 94      """ Method to stride through the data matrix, extracting the outer array with nr of elements as Column length. """ 
 95   
 96      # Extract shapes from data. 
 97      NE, NS, NM, NO, ND, Col = data.shape 
 98   
 99      # Calculate how many small matrices. 
100      Nr_mat = NE * NS * NM * NO * ND 
101   
102      # Define the shape for the stride view. 
103      shape = (Nr_mat, Col) 
104   
105      # Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16. 
106      itz = data.itemsize 
107   
108      # Bytes_between_elements 
109      bbe = 1 * itz 
110   
111      # Bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize. 
112      bbr = Col * itz 
113   
114      # Make a tuple of the strides. 
115      strides = (bbr, bbe) 
116   
117      # Make the stride view. 
118      data_view = as_strided(data, shape=shape, strides=strides) 
119   
120      return data_view 
121   
122   
123 -def stride_help_element_rank_NE_NS_NM_NO_ND(data): 
124      """ Method to stride through the data matrix, extracting the outer element. """ 
125   
126      # Extract shapes from data. 
127      NE, NS, NM, NO, Col = data.shape 
128   
129      # Calculate how many small matrices. 
130      Nr_mat = NE * NS * NM * NO * Col 
131   
132      # Define the shape for the stride view. 
133      shape = (Nr_mat, 1) 
134   
135      # Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16. 
136      itz = data.itemsize 
137   
138      # FIXME: Explain this. 
139      bbe = Col * itz 
140   
141      # FIXME: Explain this. 
142      bbr = 1 * itz 
143   
144      # Make a tuple of the strides. 
145      strides = (bbr, bbe) 
146   
147      # Make the stride view. 
148      data_view = as_strided(data, shape=shape, strides=strides) 
149   
150      return data_view 
151   
152   
153 -def stride_help_square_matrix_rank_NE_NS_NM_NO_ND_x_x(data): 
154      """ Helper function calculate the strides to go through the data matrix, extracting the outer Row X Col matrix. """ 
155   
156      # Extract shapes from data. 
157      NE, NS, NM, NO, ND, Row, Col = data.shape 
158   
159      # Calculate how many small matrices. 
160      Nr_mat = NE * NS * NM * NO * ND 
161   
162      # Define the shape for the stride view. 
163      shape = (Nr_mat, Row, Col) 
164   
165      # Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16. 
166      itz = data.itemsize 
167   
168      # Bytes_between_elements 
169      bbe = 1 * itz 
170   
171      # Bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize. 
172      bbr = Col * itz 
173   
174      # Bytes between matrices.  The byte distance is separated by the number of rows. 
175      bbm = Row * Col * itz 
176   
177      # Make a tuple of the strides. 
178      strides = (bbm, bbr, bbe) 
179   
180      # Make the stride view. 
181      data_view = as_strided(data, shape=shape, strides=strides) 
182   
183      return data_view 
184