1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 """Module for the calculation of the matrix exponential, for higher dimensional data."""
25
26
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
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
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
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
85 if len(data_array.shape) == 6:
86
87 NE, NS, NM, NO, ND, Col = data_array.shape
88
89 else:
90
91 NS, NM, NO, ND, Col = data_array.shape
92
93
94 NE = 1
95
96
97 Nr_mat = NE * NS * NM * NO * ND
98
99
100 shape = (Nr_mat, Col)
101
102
103 itz = data_array.itemsize
104
105
106 bbe = 1 * itz
107
108
109 bbr = Col * itz
110
111
112 strides = (bbr, bbe)
113
114
115 data_view = as_strided(data_array, shape=shape, strides=strides)
116
117 return data_view
118
119
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
130 if len(data_array.shape) == 7:
131
132 NE, NS, NM, NO, ND, Row, Col = data_array.shape
133
134 else:
135
136 NS, NM, NO, ND, Row, Col = data_array.shape
137
138
139 NE = 1
140
141
142 Nr_mat = NE * NS * NM * NO * ND
143
144
145 shape = (Nr_mat, Row, Col)
146
147
148 itz = data_array.itemsize
149
150
151 bbe = 1 * itz
152
153
154 bbr = Col * itz
155
156
157 bbm = Row * bbr
158
159
160 strides = (bbm, bbr, bbe)
161
162
163 data_view = as_strided(data_array, shape=shape, strides=strides)
164
165 return data_view
166
167
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
182 if len(A.shape) == 7:
183
184 NE, NS, NM, NO, ND, Row, Col = A.shape
185
186 else:
187
188 NS, NM, NO, ND, Row, Col = A.shape
189
190
191 NE = None
192
193
194 if dtype != None:
195 dtype_mat = A.dtype
196
197
198 if dtype_mat != dtype:
199
200 A = A.astype(dtype)
201
202
203 complex_flag = any(iscomplex(A))
204
205
206 if float(version.version[:3]) < 1.8:
207
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
220 A_view = data_view_via_striding_array_row_col(data_array=A)
221
222
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
227 for A_i, index_i in zip(A_view, index_view):
228
229 W_i, V_i = eig(A_i)
230
231
232 W_exp_i = exp(W_i)
233
234
235 eye_mat_i = eye(Row)
236
237
238
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
245 dot_V_W_i = dot( V_i, W_exp_diag_i)
246
247
248 inv_V_i = inv(V_i)
249
250
251 eA_i = dot(dot_V_W_i, inv_V_i)
252
253
254
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
264 W, V = eig(A)
265
266
267
268
269
270
271
272
273
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
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
286
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
293
294
295 dot_V_W = einsum('...ij, ...jk', V, W_exp_diag)
296
297
298 inv_V = inv(V)
299
300
301 eA = einsum('...ij, ...jk', dot_V_W, inv_V)
302
303
304 if complex_flag:
305 return array(eA)
306
307
308 else:
309 return array(eA.real)
310
311
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
326 NS, NM, NO, ND, Row, Col = A.shape
327
328
329 if dtype != None:
330 dtype_mat = A.dtype
331
332
333 if dtype_mat != dtype:
334
335 A = A.astype(dtype)
336
337
338 complex_flag = any(iscomplex(A))
339
340
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
347 A_view = data_view_via_striding_array_row_col(data_array=A)
348
349
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
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
361 a = 1
362 b = -a11 - a22
363 c = a11 * a22 - a12 * a21
364 dis = b**2 - 4*a*c
365
366
367
368
369
370
371 l1 = (-b + dis) / (2*a)
372 l2 = (-b - dis) / (2*a)
373
374
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
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
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
395 dot_V_W = dot(V_m, W_m)
396 eA_m = dot(dot_V_W, V_inv_m)
397
398
399
400
401 si, mi, oi, di = index_i
402
403
404 eA_mat[si, mi, oi, di, :] = eA_m
405
406 return eA_mat
407