Source Code for Module lib.dispersion.matrix_power

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