1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 from copy import deepcopy
24 from numpy import array, arange, asarray, int32, float64, max, ones, pi
25 from unittest import TestCase
26
27
28 from lib.dispersion.cr72 import r2eff_CR72
29
30 """The 1H gyromagnetic ratio."""
31 g1H = 26.7522212 * 1e7
32
33 """The 15N gyromagnetic ratio."""
34 g15N = -2.7126 * 1e7
35
36
38 """Unit tests for the lib.dispersion.cr72 relax module."""
39
41 """Set up for all unit tests."""
42
43
44 self.r2 = None
45 self.r2a = 5.0
46 self.r2b = 10.0
47 self.dw = 3.0
48 self.pA = 0.9
49 self.kex = 1000.0
50 self.spins_params = ['r2a', 'r2b', 'dw', 'pA', 'kex']
51
52
53 self.model = "CR72 full"
54 self.num_spins = 2
55 self.fields = array([800. * 1E6])
56
57
58 self.exp_type = ['SQ CPMG']
59 self.offset = [0]
60
61
62 self.relax_times = self.fields / (100 * 100. *1E6 )
63 self.ncycs = []
64 self.points = []
65 self.value = []
66 self.error = []
67 self.num_disp_points = []
68
69 for i in range(len(self.fields)):
70 ncyc = arange(2, 1000. * self.relax_times[i], 4)
71
72 self.ncycs.append(ncyc)
73
74
75 cpmg_point = ncyc / self.relax_times[i]
76
77 nr_cpmg_points = len(cpmg_point)
78 self.num_disp_points.append(nr_cpmg_points)
79
80 self.points.append(list(cpmg_point))
81 self.value.append([2.0]*len(cpmg_point))
82 self.error.append([1.0]*len(cpmg_point))
83
84
86 """Calculate and check the R2eff values."""
87
88
89 self.params = self.assemble_param_vector(r2=self.r2, r2a=self.r2a, r2b=self.r2b, dw=self.dw, pA=self.pA, kex=self.kex, spins_params=self.spins_params)
90
91
92 values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times, offsets = self.return_r2eff_arrays()
93
94
95
96 end_index = []
97
98 end_index.append(len(self.exp_type) * self.num_spins * len(self.fields))
99 if self.model in ["CR72 full"]:
100 end_index.append(2 * len(self.exp_type) * self.num_spins * len(self.fields))
101
102 end_index.append(end_index[-1] + self.num_spins)
103
104
105 R20 = self.params[:end_index[1]].reshape(self.num_spins*2, len(self.fields))
106 R20A = R20[::2].flatten()
107 R20B = R20[1::2].flatten()
108 dw = self.params[end_index[1]:end_index[2]]
109 pA = self.params[end_index[2]]
110 kex = self.params[end_index[2]+1]
111
112
113 self.back_calc = deepcopy(values)
114
115
116
117 back_calc_shape = list( asarray(self.back_calc).shape )[:4]
118
119
120
121 self.max_num_disp_points = max(self.num_disp_points)
122
123
124
125
126
127 self.R20A_a = ones(back_calc_shape + [self.max_num_disp_points])
128 self.R20B_a = ones(back_calc_shape + [self.max_num_disp_points])
129 self.dw_frq_a = ones(back_calc_shape + [self.max_num_disp_points])
130 self.cpmg_frqs_a = ones(back_calc_shape + [self.max_num_disp_points])
131 self.num_disp_points_a = ones(back_calc_shape + [self.max_num_disp_points])
132 self.back_calc_a = ones(back_calc_shape + [self.max_num_disp_points])
133
134
135 for si in range(self.num_spins):
136
137 for mi in range(len(self.fields)):
138
139 num_disp_points = self.num_disp_points[mi]
140
141
142 self.cpmg_frqs_a[0][si][mi][0][:num_disp_points] = cpmg_frqs[0][mi][0]
143 self.num_disp_points_a[0][si][mi][0][:num_disp_points] = self.num_disp_points[mi]
144
145
146
147
148 for si in range(self.num_spins):
149
150 for mi in range(len(self.fields)):
151
152 num_disp_points = len(cpmg_frqs[0][mi][0])
153
154
155 r20_index = mi + si*len(self.fields)
156
157
158 self.R20A_a[0][si][mi][0] = array( [R20A[r20_index]] * self.max_num_disp_points, float64)
159 self.R20B_a[0][si][mi][0] = array( [R20B[r20_index]] * self.max_num_disp_points, float64)
160
161
162 dw_frq = dw[si] * frqs[0][si][mi]
163
164
165 self.dw_frq_a[0][si][mi][0] = array( [dw_frq] * self.max_num_disp_points, float64)
166
167
168 r2eff_CR72(r20a_orig=self.R20A_a, r20b_orig=self.R20B_a, dw_orig=self.dw_frq_a, r20a=self.R20A_a, r20b=self.R20B_a, pA=pA, dw=self.dw_frq_a, kex=kex, cpmg_frqs=self.cpmg_frqs_a, back_calc=self.back_calc_a)
169
170
171
172
173 for si in range(self.num_spins):
174
175 for mi in range(len(self.fields)):
176
177 num_disp_points = self.num_disp_points[mi]
178
179
180 self.back_calc[0][si][mi][0][:] = self.back_calc_a[0][si][mi][0][:num_disp_points]
181
182
183 for di in range(num_disp_points):
184 self.assertAlmostEqual(self.back_calc[0][si][mi][0][di], self.R20A_a[0][si][mi][0][di])
185
186
188 """Return numpy arrays of the R2eff/R1rho values and errors.
189
190 @return: The numpy array structures of the R2eff/R1rho values, errors, missing data, and corresponding Larmor frequencies. For each structure, the first dimension corresponds to the experiment types, the second the spins of a spin block, the third to the spectrometer field strength, and the fourth is the dispersion points. For the Larmor frequency structure, the fourth dimension is omitted. For R1rho-type data, an offset dimension is inserted between the spectrometer field strength and the dispersion points.
191 @rtype: lists of numpy float arrays, lists of numpy float arrays, lists of numpy float arrays, numpy rank-2 int array
192 """
193
194
195 exp_types = []
196 values = []
197 errors = []
198 missing = []
199 frqs = []
200 frqs_H = []
201 relax_times = []
202 offsets = []
203 for ei in range(len(self.exp_type)):
204 values.append([])
205 errors.append([])
206 missing.append([])
207 frqs.append([])
208 frqs_H.append([])
209 relax_times.append([])
210 offsets.append([])
211 for si in range(self.num_spins):
212 values[ei].append([])
213 errors[ei].append([])
214 missing[ei].append([])
215 frqs[ei].append([])
216 frqs_H[ei].append([])
217 offsets[ei].append([])
218 for mi in range(len(self.fields)):
219 values[ei][si].append([])
220 errors[ei][si].append([])
221 missing[ei][si].append([])
222 frqs[ei][si].append(0.0)
223 frqs_H[ei][si].append(0.0)
224 offsets[ei][si].append([])
225 for oi in range(len(self.offset)):
226 values[ei][si][mi].append([])
227 errors[ei][si][mi].append([])
228 missing[ei][si][mi].append([])
229 offsets[ei][si][mi].append([])
230 for mi in range(len(self.fields)):
231 relax_times[ei].append(None)
232
233 cpmg_frqs = []
234 for ei in range(len(self.exp_type)):
235 cpmg_frqs.append([])
236 for mi in range(len(self.fields)):
237 cpmg_frqs[ei].append([])
238 for oi in range(len(self.offset)):
239
240 cpmg_frqs[ei][mi].append([])
241
242
243
244 si = 0
245 for spin_index in range(self.num_spins):
246 data_flag = True
247
248 for ei in range(len(self.exp_type)):
249 exp_type = self.exp_type[ei]
250
251 if exp_type not in exp_types:
252 exp_types.append(exp_type)
253
254 for mi in range(len(self.fields)):
255
256 frq = self.fields[mi]
257
258
259 frqs[ei][si][mi] = 2.0 * pi * frq / g1H * g15N * 1e-6
260
261
262 cpmg_frqs[ei][mi][oi] = self.points[mi]
263
264 for oi in range(len(self.offset)):
265 for di in range(len(self.points[mi])):
266
267 missing[ei][si][mi][oi].append(0)
268
269
270 values[ei][si][mi][oi].append(self.value[mi][di])
271
272 errors[ei][si][mi][oi].append(self.error[mi][di])
273
274
275
276 relax_time = self.relax_times[mi]
277
278
279 relax_times[ei][mi] = relax_time
280
281
282 si += 1
283
284
285 relax_times = array(relax_times, float64)
286 for ei in range(len(self.exp_type)):
287 for si in range(self.num_spins):
288 for mi in range(len(self.fields)):
289 for oi in range(len(self.offset)):
290 values[ei][si][mi][oi] = array(values[ei][si][mi][oi], float64)
291 errors[ei][si][mi][oi] = array(errors[ei][si][mi][oi], float64)
292 missing[ei][si][mi][oi] = array(missing[ei][si][mi][oi], int32)
293
294
295 return values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times, offsets
296
297
298 - def assemble_param_vector(self, r2=None, r2a=None, r2b=None, dw=None, pA=None, kex=None, spins_params=None):
299 """Assemble the dispersion relaxation dispersion curve fitting parameter vector.
300
301 @keyword r2: The transversal relaxation rate.
302 @type r2: float
303 @keyword r2a: The transversal relaxation rate for state A in the absence of exchange.
304 @type r2a: float
305 @keyword r2b: The transversal relaxation rate for state B in the absence of exchange.
306 @type r2b: float
307 @keyword dw: The chemical exchange difference between states A and B in ppm.
308 @type dw: float
309 @keyword pA: The population of state A.
310 @type pA: float
311 @keyword kex: The rate of exchange.
312 @type kex: float
313 @keyword spins_params: List of parameter strings used in dispersion model.
314 @type spins_params: array of strings
315 @return: An array of the parameter values of the dispersion relaxation model.
316 @rtype: numpy float array
317 """
318
319
320 param_vector = []
321
322
323 for param_name, spin_index, mi in self.loop_parameters(spins_params=spins_params):
324 if param_name == 'r2':
325 value = r2
326 value = value + mi + spin_index*0.1
327 elif param_name == 'r2a':
328 value = r2a
329 value = value + mi+ spin_index*0.1
330 elif param_name == 'r2b':
331 value = r2b
332 value = value + mi + spin_index*0.1
333 elif param_name == 'dw':
334 value = dw
335 elif param_name == 'pA':
336 value = pA
337 elif param_name == 'kex':
338 value = kex
339
340
341 param_vector.append(value)
342
343
344 return array(param_vector, float64)
345
346
348 """Generator function for looping of the parameters of the cluster.
349
350 @keyword spins_params: List of parameter strings used in dispersion model.
351 @type spins_params: array of strings
352 @return: The parameter name.
353 @rtype: str
354 """
355
356
357
358 for spin_index in range(self.num_spins):
359
360
361 if 'r2' in spins_params:
362 for ei in range(len(self.exp_type)):
363 for mi in range(len(self.fields)):
364 yield 'r2', spin_index, mi
365
366
367 if 'r2a' in spins_params:
368 for ei in range(len(self.exp_type)):
369 for mi in range(len(self.fields)):
370 yield 'r2a', spin_index, mi
371
372
373
374 if 'r2b' in spins_params:
375 for ei in range(len(self.exp_type)):
376 for mi in range(len(self.fields)):
377 yield 'r2b', spin_index, mi
378
379
380
381 for spin_index in range(self.num_spins):
382
383 if 'dw' in spins_params:
384 yield 'dw', spin_index, 0
385
386
387 for param in spins_params:
388 if not param in ['r2', 'r2a', 'r2b', 'phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2', 'dw', 'dw_AB', 'dw_BC', 'dw_AB', 'dwH', 'dwH_AB', 'dwH_BC', 'dwH_AB']:
389 if param == 'pA':
390 yield 'pA', 0, 0
391 elif param == 'kex':
392 yield 'kex', 0, 0
393
394
396 """Test the r2eff_cr72() function for no exchange when dw = 0.0."""
397
398
399 self.dw = 0.0
400
401
402 self.calc_r2eff()
403
404
406 """Test the r2eff_cr72() function for no exchange when pA = 1.0."""
407
408
409 self.pA = 1.0
410
411
412 self.calc_r2eff()
413
414
416 """Test the r2eff_cr72() function for no exchange when kex = 0.0."""
417
418
419 self.kex = 0.0
420
421
422 self.calc_r2eff()
423
424
426 """Test the r2eff_cr72() function for no exchange when dw = 0.0 and pA = 1.0."""
427
428
429 self.pA = 1.0
430 self.dw = 0.0
431
432
433 self.calc_r2eff()
434
435
437 """Test the r2eff_cr72() function for no exchange when dw = 0.0 and kex = 0.0."""
438
439
440 self.dw = 0.0
441 self.kex = 0.0
442
443
444 self.calc_r2eff()
445
446
448 """Test the r2eff_cr72() function for no exchange when pA = 1.0 and kex = 0.0."""
449
450
451 self.pA = 1.0
452 self.kex = 0.0
453
454
455 self.calc_r2eff()
456
457
459 """Test the r2eff_cr72() function for no exchange when dw = 0.0, pA = 1.0, and kex = 0.0."""
460
461
462 self.dw = 0.0
463 self.pA = 1.0
464 self.kex = 0.0
465
466
467 self.calc_r2eff()
468
469
471 """Test the r2eff_cr72() function for no exchange when kex = 1e7."""
472
473
474 self.kex = 1e7
475
476
477 self.calc_r2eff()
478