lib.dispersion.matrix

1 ############################################################################### 2 # # 3 # Copyright (C) 2013 Edward d'Auvergne # 4 # Copyright (C) 2014 Troels E. Linnet # 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 exponential, for higher dimensional data.""" 25 26 # Python module imports. 27 from numpy import array, any, complex128, dot, einsum, eye, exp, iscomplex, int16, newaxis, multiply, tile, sqrt, version, zeros 28 from numpy.lib.stride_tricks import as_strided 29 from numpy.linalg import eig, inv 30 31

32 -def create_index(NE=None, NS=None, NM=None, NO=None, ND=None):

33 """Method to create the helper index numpy array, to help figuring out the indices to store in the exchange data matrix. 34 35 @keyword NE: The total number of experiment types. 36 @type NE: None or int 37 @keyword NS: The total number of spins of the spin cluster. 38 @type NS: int 39 @keyword NM: The total number of magnetic field strengths. 40 @type NM: int 41 @keyword NO: The total number of spin-lock offsets. 42 @type NO: int 43 @keyword ND: The total number of dispersion points (either the spin-lock field strength or the nu_CPMG frequency). 44 @type ND: int 45 @return: The numpy array for containing index indices for storing in the strided exchange data matrix. 46 @rtype: numpy int array of rank [NE][NS][NM][NO][ND][5] or [NS][NM][NO][ND][4]. 47 """ 48 49 # Make array to store index. 50 if NE != None: 51 index = zeros([NE, NS, NM, NO, ND, 5], int16) 52 53 else: 54 index = zeros([NS, NM, NO, ND, 4], int16) 55 56 # Make indices for storing in data matrix. 57 if NE != None: 58 for ei in range(NE): 59 for si in range(NS): 60 for mi in range(NM): 61 for oi in range(NO): 62 for di in range(ND): 63 index[ei, si, mi, oi, di] = ei, si, mi, oi, di 64 65 else: 66 for si in range(NS): 67 for mi in range(NM): 68 for oi in range(NO): 69 for di in range(ND): 70 index[si, mi, oi, di] = si, mi, oi, di 71 72 return index

73 74

75 -def data_view_via_striding_array_col(data_array=None):

76 """Method to stride through the data matrix, extracting the outer array with nr of elements as Column length. 77 78 @keyword data: The numpy data array to stride through. 79 @type data: numpy array of rank [NE][NS][NM][NO][ND][Col] or [NS][NM][NO][ND][Col]. 80 @return: The data view of the full numpy array, returned as a numpy array with number of small numpy arrays corresponding to Nr_mat=NE*NS*NM*NO*ND or Nr_mat=NS*NM*NO*ND, where each small array has size Col. 81 @rtype: numpy array of rank [NE*NS*NM*NO*ND][Col] or [NS*NM*NO*ND][Col]. 82 """ 83 84 # Get the expected shape of the higher dimensional column numpy array. 85 if len(data_array.shape) == 6: 86 # Extract shapes from data. 87 NE, NS, NM, NO, ND, Col = data_array.shape 88 89 else: 90 # Extract shapes from data. 91 NS, NM, NO, ND, Col = data_array.shape 92 93 # Set NE to 1. 94 NE = 1 95 96 # Calculate how many small matrices. 97 Nr_mat = NE * NS * NM * NO * ND 98 99 # Define the shape for the stride view. 100 shape = (Nr_mat, Col) 101 102 # Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16. 103 itz = data_array.itemsize 104 105 # Bytes_between_elements 106 bbe = 1 * itz 107 108 # Bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize. 109 bbr = Col * itz 110 111 # Make a tuple of the strides. 112 strides = (bbr, bbe) 113 114 # Make the stride view. 115 data_view = as_strided(data_array, shape=shape, strides=strides) 116 117 return data_view

118 119

120 -def data_view_via_striding_array_row_col(data_array=None):

121 """Method to stride through the data matrix, extracting the outer matrix with nr of elements as Row X Column length. Row and Col should have same size. 122 123 @keyword data: The numpy data array to stride through. 124 @type data: numpy array of rank [NE][NS][NM][NO][ND][Row][Col] or [NS][NM][NO][ND][Row][Col]. 125 @return: The data view of the full numpy array, returned as a numpy array with number of small numpy arrays corresponding to Nr_mat=NE*NS*NM*NO*ND or Nr_mat=NS*NM*NO*ND, where each small array has size Col. 126 @rtype: numpy array of rank [NE*NS*NM*NO*ND][Col] or [NS*NM*NO*ND][Col]. 127 """ 128 129 # Get the expected shape of the higher dimensional column numpy array. 130 if len(data_array.shape) == 7: 131 # Extract shapes from data. 132 NE, NS, NM, NO, ND, Row, Col = data_array.shape 133 134 else: 135 # Extract shapes from data. 136 NS, NM, NO, ND, Row, Col = data_array.shape 137 138 # Set NE to 1. 139 NE = 1 140 141 # Calculate how many small matrices. 142 Nr_mat = NE * NS * NM * NO * ND 143 144 # Define the shape for the stride view. 145 shape = (Nr_mat, Row, Col) 146 147 # Get itemsize, Length of one array element in bytes. Depends on dtype. float64=8, complex128=16. 148 itz = data_array.itemsize 149 150 # Bytes_between_elements 151 bbe = 1 * itz 152 153 # Bytes between row. The distance in bytes to next row is number of Columns elements multiplied with itemsize. 154 bbr = Col * itz 155 156 # Bytes between matrices. The byte distance is separated by the number of rows. 157 bbm = Row * bbr 158 159 # Make a tuple of the strides. 160 strides = (bbm, bbr, bbe) 161 162 # Make the stride view. 163 data_view = as_strided(data_array, shape=shape, strides=strides) 164 165 return data_view

166 167

168 -def matrix_exponential(A, dtype=None):

169 """Calculate the exact matrix exponential using the eigenvalue decomposition approach, for higher dimensional data. This of dimension [NE][NS][NM][NO][ND][X][X] or [NS][NM][NO][ND][X][X]. 170 171 Here X is the Row and Column length, of the outer square matrix. 172 173 @param A: The square matrix to calculate the matrix exponential of. 174 @type A: numpy float array of rank [NE][NS][NM][NO][ND][X][X] 175 @param dtype: If provided, forces the calculation to use the data type specified. 176 @type dtype: data-type, optional 177 @return: The matrix exponential. This will have the same dimensionality as the A matrix. 178 @rtype: numpy float array of rank [NE][NS][NM][NO][ND][X][X] 179 """ 180 181 # Get the expected shape of the higher dimensional column numpy array. 182 if len(A.shape) == 7: 183 # Extract shapes from data. 184 NE, NS, NM, NO, ND, Row, Col = A.shape 185 186 else: 187 # Extract shapes from data. 188 NS, NM, NO, ND, Row, Col = A.shape 189 190 # Set NE to None. 191 NE = None 192 193 # Convert dtype, if specified. 194 if dtype != None: 195 dtype_mat = A.dtype 196 197 # If the dtype is different from the input. 198 if dtype_mat != dtype: 199 # This needs to be made as a copy. 200 A = A.astype(dtype) 201 202 # Is the original matrix real? 203 complex_flag = any(iscomplex(A)) 204 205 # If numpy is under 1.8, there would be a need to do eig(A) per matrix. 206 if float(version.version[:3]) < 1.8: 207 # Make array to store results 208 if NE != None: 209 if dtype != None: 210 eA = zeros([NE, NS, NM, NO, ND, Row, Col], dtype=dtype) 211 else: 212 eA = zeros([NE, NS, NM, NO, ND, Row, Col], dtype=complex128) 213 else: 214 if dtype != None: 215 eA = zeros([NS, NM, NO, ND, Row, Col], dtype=dtype) 216 else: 217 eA = zeros([NS, NM, NO, ND, Row, Col], dtype=complex128) 218 219 # Get the data view, from the helper function. 220 A_view = data_view_via_striding_array_row_col(data_array=A) 221 222 # Create index view. 223 index = create_index(NE=NE, NS=NS, NM=NM, NO=NO, ND=ND) 224 index_view = data_view_via_striding_array_col(data_array=index) 225 226 # Zip them together and iterate over them. 227 for A_i, index_i in zip(A_view, index_view): 228 # The eigenvalue decomposition. 229 W_i, V_i = eig(A_i) 230 231 # Calculate the exponential. 232 W_exp_i = exp(W_i) 233 234 # Make a eye matrix. 235 eye_mat_i = eye(Row) 236 237 # Transform it to a diagonal matrix, with elements from vector down the diagonal. 238 # Use the dtype, if specified. 239 if dtype != None: 240 W_exp_diag_i = multiply(W_exp_i, eye_mat_i, dtype=dtype ) 241 else: 242 W_exp_diag_i = multiply(W_exp_i, eye_mat_i) 243 244 # Make dot product. 245 dot_V_W_i = dot( V_i, W_exp_diag_i) 246 247 # Compute the (multiplicative) inverse of a matrix. 248 inv_V_i = inv(V_i) 249 250 # Calculate the exact exponential. 251 eA_i = dot(dot_V_W_i, inv_V_i) 252 253 # Save results. 254 # Extract index from index_view. 255 if NE != None: 256 ei, si, mi, oi, di = index_i 257 eA[ei, si, mi, oi, di, :] = eA_i 258 else: 259 si, mi, oi, di = index_i 260 eA[si, mi, oi, di, :] = eA_i 261 262 else: 263 # The eigenvalue decomposition. 264 W, V = eig(A) 265 266 # W: The eigenvalues, each repeated according to its multiplicity. 267 # The eigenvalues are not necessarily ordered. 268 # The resulting array will be always be of complex type. Shape [NE][NS][NM][NO][ND][X] 269 # V: The normalized (unit 'length') eigenvectors, such that the column v[:,i] 270 # is the eigenvector corresponding to the eigenvalue w[i]. Shape [NE][NS][NM][NO][ND][X][X] 271 272 # Calculate the exponential of all elements in the input array. Shape [NE][NS][NM][NO][ND][X] 273 # Add one axis, to allow for broadcasting multiplication. 274 if NE != None: 275 W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1) 276 else: 277 W_exp = exp(W).reshape(NS, NM, NO, ND, Row, 1) 278 279 # Make a eye matrix, with Shape [NE][NS][NM][NO][ND][X][X] 280 if NE != None: 281 eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, newaxis, ...], (NE, NS, NM, NO, ND, 1, 1) ) 282 else: 283 eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, ...], (NS, NM, NO, ND, 1, 1) ) 284 285 # Transform it to a diagonal matrix, with elements from vector down the diagonal. 286 # Use the dtype, if specified. 287 if dtype != None: 288 W_exp_diag = multiply(W_exp, eye_mat, dtype=dtype ) 289 else: 290 W_exp_diag = multiply(W_exp, eye_mat) 291 292 # Make dot products for higher dimension. 293 # "...", the Ellipsis notation, is designed to mean to insert as many full slices (:) 294 # to extend the multi-dimensional slice to all dimensions. 295 dot_V_W = einsum('...ij, ...jk', V, W_exp_diag) 296 297 # Compute the (multiplicative) inverse of a matrix. 298 inv_V = inv(V) 299 300 # Calculate the exact exponential. 301 eA = einsum('...ij, ...jk', dot_V_W, inv_V) 302 303 # Return the complex matrix. 304 if complex_flag: 305 return array(eA) 306 307 # Return only the real part. 308 else: 309 return array(eA.real)

310 311

312 -def matrix_exponential_rank_NS_NM_NO_ND_2_2(A, dtype=None):

313 """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]. 314 315 Here X is the Row and Column length, of the outer square matrix. 316 317 @param A: The square matrix to calculate the matrix exponential of. 318 @type A: numpy float array of rank [NS][NM][NO][ND][2][2] 319 @param dtype: If provided, forces the calculation to use the data type specified. 320 @type dtype: data-type, optional 321 @return: The matrix exponential. This will have the same dimensionality as the A matrix. 322 @rtype: numpy float array of rank [NS][NM][NO][ND][2][2] 323 """ 324 325 # Extract shapes from data. 326 NS, NM, NO, ND, Row, Col = A.shape 327 328 # Convert dtype, if specified. 329 if dtype != None: 330 dtype_mat = A.dtype 331 332 # If the dtype is different from the input. 333 if dtype_mat != dtype: 334 # This needs to be made as a copy. 335 A = A.astype(dtype) 336 337 # Is the original matrix real? 338 complex_flag = any(iscomplex(A)) 339 340 # Make array to store results 341 if dtype != None: 342 eA_mat = zeros([NS, NM, NO, ND, Row, Col], dtype) 343 else: 344 eA_mat = zeros([NS, NM, NO, ND, Row, Col], dtype) 345 346 # Get the data view, from the helper function. 347 A_view = data_view_via_striding_array_row_col(data_array=A) 348 349 # The index view could be pre formed in init. 350 index = create_index(NS=NS, NM=NM, NO=NO, ND=ND) 351 index_view = data_view_via_striding_array_col(data_array=index) 352 353 # Zip them together and iterate over them. 354 for A_i, index_i in zip(A_view, index_view): 355 a11 = A_i[0, 0] 356 a12 = A_i[0, 1] 357 a21 = A_i[1, 0] 358 a22 = A_i[1, 1] 359 360 # Discriminant 361 a = 1 362 b = -a11 - a22 363 c = a11 * a22 - a12 * a21 364 dis = b**2 - 4*a*c 365 366 # If dis is positive: two distinct real roots 367 # If dis is zero: one distinct real roots 368 # If dis is negative: two complex roots 369 370 # Eigenvalues lambda_1, lambda_2. 371 l1 = (-b + dis) / (2*a) 372 l2 = (-b - dis) / (2*a) 373 374 # Define M: M = V^-1 * [ [l1 0], [0 l2] ] * V 375 W_m = array([ [l1, 0], [0, l2] ]) 376 377 v1_1 = 1 378 v1_2 = (l1 - a11) / a12 379 380 v2_1 = 1 381 v2_2 = (l2 - a11) / a12 382 383 # normalized eigenvector 384 v1_1_norm = v1_1 / (sqrt(v1_1**2 + v1_2**2) ) 385 v1_2_norm = v1_2 / (sqrt(v1_1**2 + v1_2**2) ) 386 v2_1_norm = v2_1 / (sqrt(v2_1**2 + v2_2**2) ) 387 v2_2_norm = v2_2 / (sqrt(v2_1**2 + v2_2**2) ) 388 389 #V_m = array([ [v1_norm], [v2_norm] ]) 390 V_m = array([ [v1_1_norm, v2_1_norm], [v1_2_norm, v2_2_norm] ]) 391 392 V_inv_m = inv(V_m) 393 394 # Calculate the exact exponential. 395 dot_V_W = dot(V_m, W_m) 396 eA_m = dot(dot_V_W, V_inv_m) 397 398 # Save results. 399 400 # Extract index from index_view. 401 si, mi, oi, di = index_i 402 403 # Store the result. 404 eA_mat[si, mi, oi, di, :] = eA_m 405 406 return eA_mat

407

Source Code for Module lib.dispersion.matrix_exponential