1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 """Target functions for relaxation dispersion."""
26
27
28 from copy import deepcopy
29 from numpy import all, arctan2, cos, dot, float64, int16, isfinite, max, multiply, ones, rollaxis, pi, sin, sum, zeros
30 from numpy.ma import masked_equal
31
32
33 from lib.dispersion.b14 import r2eff_B14
34 from lib.dispersion.cr72 import r2eff_CR72
35 from lib.dispersion.dpl94 import r1rho_DPL94
36 from lib.dispersion.it99 import r2eff_IT99
37 from lib.dispersion.lm63 import r2eff_LM63
38 from lib.dispersion.lm63_3site import r2eff_LM63_3site
39 from lib.dispersion.m61 import r1rho_M61
40 from lib.dispersion.m61b import r1rho_M61b
41 from lib.dispersion.mp05 import r1rho_MP05
42 from lib.dispersion.mmq_cr72 import r2eff_mmq_cr72
43 from lib.dispersion.ns_cpmg_2site_3d import r2eff_ns_cpmg_2site_3D
44 from lib.dispersion.ns_cpmg_2site_expanded import r2eff_ns_cpmg_2site_expanded
45 from lib.dispersion.ns_cpmg_2site_star import r2eff_ns_cpmg_2site_star
46 from lib.dispersion.ns_mmq_3site import r2eff_ns_mmq_3site_mq, r2eff_ns_mmq_3site_sq_dq_zq
47 from lib.dispersion.ns_mmq_2site import r2eff_ns_mmq_2site_mq, r2eff_ns_mmq_2site_sq_dq_zq
48 from lib.dispersion.ns_r1rho_2site import ns_r1rho_2site
49 from lib.dispersion.ns_r1rho_3site import ns_r1rho_3site
50 from lib.dispersion.ns_matrices import r180x_3d
51 from lib.dispersion.tp02 import r1rho_TP02
52 from lib.dispersion.tap03 import r1rho_TAP03
53 from lib.dispersion.tsmfk01 import r2eff_TSMFK01
54 from lib.dispersion.variables import EXP_TYPE_CPMG_DQ, EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_ZQ, EXP_TYPE_LIST_CPMG, EXP_TYPE_R1RHO, MODEL_B14, MODEL_B14_FULL, MODEL_CR72, MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, MODEL_LIST_CPMG, MODEL_LIST_FULL, MODEL_LIST_DW_MIX_DOUBLE, MODEL_LIST_DW_MIX_QUADRUPLE, MODEL_LIST_INV_RELAX_TIMES, MODEL_LIST_R20B, MODEL_LIST_MMQ, MODEL_LIST_MQ_CPMG, MODEL_LIST_R1RHO, MODEL_LIST_R1RHO_OFF_RES, MODEL_LM63, MODEL_LM63_3SITE, MODEL_M61, MODEL_M61B, MODEL_MP05, MODEL_MMQ_CR72, MODEL_NOREX, MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_EXPANDED, MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_MMQ_2SITE, MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR, MODEL_NS_R1RHO_2SITE, MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR, MODEL_TAP03, MODEL_TP02, MODEL_TSMFK01
55 from lib.errors import RelaxError
56 from lib.float import isNaN
57 from target_functions.chi2 import chi2_rankN
58
59
61 - def __init__(self, model=None, num_params=None, num_spins=None, num_frq=None, exp_types=None, values=None, errors=None, missing=None, frqs=None, frqs_H=None, cpmg_frqs=None, spin_lock_nu1=None, chemical_shifts=None, offset=None, tilt_angles=None, r1=None, relax_times=None, scaling_matrix=None, recalc_tau=True, r1_fit=False):
62 """Relaxation dispersion target functions for optimisation.
63
64 Models
65 ======
66
67 The following analytic models are currently supported:
68
69 - 'No Rex': The model for no chemical exchange relaxation.
70 - 'LM63': The Luz and Meiboom (1963) 2-site fast exchange model.
71 - 'LM63 3-site': The Luz and Meiboom (1963) 3-site fast exchange model.
72 - 'CR72': The reduced Carver and Richards (1972) 2-site model for all time scales with R20A = R20B.
73 - 'CR72 full': The full Carver and Richards (1972) 2-site model for all time scales.
74 - 'IT99': The Ishima and Torchia (1999) 2-site model for all time scales with skewed populations (pA >> pB).
75 - 'TSMFK01': The Tollinger et al. (2001) 2-site very-slow exchange model, range of microsecond to second time scale.
76 - 'B14': The Baldwin (2014) 2-site exact solution model for all time scales with R20A = R20B..
77 - 'B14 full': The Baldwin (2014) 2-site exact solution model for all time scales.
78 - 'M61': The Meiboom (1961) 2-site fast exchange model for R1rho-type experiments.
79 - 'DPL94': The Davis, Perlman and London (1994) 2-site fast exchange model for R1rho-type experiments.
80 - 'M61 skew': The Meiboom (1961) on-resonance 2-site model with skewed populations (pA >> pB) for R1rho-type experiments.
81 - 'TP02': The Trott and Palmer (2002) 2-site exchange model for R1rho-type experiments.
82 - 'TAP03': The Trott, Abergel and Palmer (2003) 2-site exchange model for R1rho-type experiments.
83 - 'MP05': The Miloushev and Palmer (2005) off-resonance 2-site exchange model for R1rho-type experiments.
84
85 The following numerical models are currently supported:
86
87 - 'NS CPMG 2-site 3D': The reduced numerical solution for the 2-site Bloch-McConnell equations for CPMG data using 3D magnetisation vectors with R20A = R20B.
88 - 'NS CPMG 2-site 3D full': The full numerical solution for the 2-site Bloch-McConnell equations for CPMG data using 3D magnetisation vectors.
89 - 'NS CPMG 2-site star': The reduced numerical solution for the 2-site Bloch-McConnell equations for CPMG data using complex conjugate matrices with R20A = R20B.
90 - 'NS CPMG 2-site star full': The full numerical solution for the 2-site Bloch-McConnell equations for CPMG data using complex conjugate matrices.
91 - 'NS CPMG 2-site expanded': The numerical solution for the 2-site Bloch-McConnell equations for CPMG data expanded using Maple by Nikolai Skrynnikov.
92 - 'NS MMQ 2-site': The numerical solution for the 2-site Bloch-McConnell equations for combined proton-heteronuclear SQ, ZD, DQ, and MQ CPMG data with R20A = R20B.
93 - 'NS MMQ 3-site linear': The numerical solution for the 3-site Bloch-McConnell equations linearised with kAC = kCA = 0 for combined proton-heteronuclear SQ, ZD, DQ, and MQ CPMG data with R20A = R20B = R20C.
94 - 'NS MMQ 3-site': The numerical solution for the 3-site Bloch-McConnell equations for combined proton-heteronuclear SQ, ZD, DQ, and MQ CPMG data with R20A = R20B = R20C.
95 - 'NS R1rho 2-site': The numerical solution for the 2-site Bloch-McConnell equations for R1rho data with R20A = R20B.
96 - 'NS R1rho 3-site linear': The numerical solution for the 3-site Bloch-McConnell equations linearised with kAC = kCA = 0 for R1rho data with R20A = R20B = R20C.
97 - 'NS R1rho 3-site': The numerical solution for the 3-site Bloch-McConnell equations for R1rho data with R20A = R20B = R20C.
98
99
100 Indices
101 =======
102
103 The data structures used in this target function class consist of many different index types. These are:
104
105 - Ei: The index for each experiment type.
106 - Si: The index for each spin of the spin cluster.
107 - Mi: The index for each magnetic field strength.
108 - Oi: The index for each spin-lock offset. In the case of CPMG-type data, this index is currently always zero.
109 - Di: The index for each dispersion point (either the spin-lock field strength or the nu_CPMG frequency).
110 - Ti: The index for each time point. This is currently unused but might change in the future.
111
112
113 Counts
114 ======
115
116 The indices in the previous section have a corresponding count:
117
118 - NE: The total number of experiment types.
119 - NS: The total number of spins of the spin cluster.
120 - NM: The total number of magnetic field strengths.
121 - NO: The total number of spin-lock offsets.
122 - ND: The total number of dispersion points (either the spin-lock field strength or the nu_CPMG frequency).
123 - NT: The total number of time points. This is currently unused but might change in the future.
124
125
126 @keyword model: The relaxation dispersion model to fit.
127 @type model: str
128 @keyword num_param: The number of parameters in the model.
129 @type num_param: int
130 @keyword num_spins: The number of spins in the cluster.
131 @type num_spins: int
132 @keyword num_frq: The number of spectrometer field strengths.
133 @type num_frq: int
134 @keyword exp_types: The list of experiment types. The dimensions are {Ei}.
135 @type exp_types: list of str
136 @keyword values: The R2eff/R1rho values. The dimensions are {Ei, Si, Mi, Oi, Di}.
137 @type values: rank-4 list of numpy rank-1 float arrays
138 @keyword errors: The R2eff/R1rho errors. The dimensions are {Ei, Si, Mi, Oi, Di}.
139 @type errors: rank-4 list of numpy rank-1 float arrays
140 @keyword missing: The data structure indicating missing R2eff/R1rho data. The dimensions are {Ei, Si, Mi, Oi, Di}.
141 @type missing: rank-4 list of numpy rank-1 int arrays
142 @keyword frqs: The spin Larmor frequencies (in MHz*2pi to speed up the ppm to rad/s conversion). The dimensions are {Ei, Si, Mi}.
143 @type frqs: rank-3 list of floats
144 @keyword frqs_H: The proton spin Larmor frequencies for the MMQ-type models (in MHz*2pi to speed up the ppm to rad/s conversion). The dimensions are {Ei, Si, Mi}.
145 @type frqs_H: rank-3 list of floats
146 @keyword cpmg_frqs: The CPMG frequencies in Hertz. This will be ignored for R1rho experiments. The dimensions are {Ei, Mi, Oi}.
147 @type cpmg_frqs: rank-3 list of floats
148 @keyword spin_lock_nu1: The spin-lock field strengths in Hertz. This will be ignored for CPMG experiments. The dimensions are {Ei, Mi, Oi}.
149 @type spin_lock_nu1: rank-3 list of floats
150 @keyword chemical_shifts: The chemical shifts in rad/s. This is only used for off-resonance R1rho models. The ppm values are not used to save computation time, therefore they must be converted to rad/s by the calling code. The dimensions are {Ei, Si, Mi}.
151 @type chemical_shifts: rank-3 list of floats
152 @keyword offset: The structure of spin-lock or hard pulse offsets in rad/s. This is only currently used for off-resonance R1rho models. The dimensions are {Ei, Si, Mi, Oi}.
153 @type offset: rank-4 list of floats
154 @keyword tilt_angles: The spin-lock rotating frame tilt angle. This is only used for off-resonance R1rho models. The dimensions are {Ei, Si, Mi, Oi, Di}.
155 @type tilt_angles: rank-5 list of floats
156 @keyword r1: The R1 relaxation rates. This is only used for off-resonance R1rho models. The dimensions are {Si, Mi}.
157 @type r1: rank-2 list of floats
158 @keyword relax_times: The experiment specific fixed time period for relaxation (in seconds). The dimensions are {Ei, Mi, Oi, Di, Ti}.
159 @type relax_times: rank-4 list of floats
160 @keyword scaling_matrix: The square and diagonal scaling matrix.
161 @type scaling_matrix: numpy rank-2 float array
162 @keyword recalc_tau: A flag which if True will cause tau_CPMG to be recalculated to remove user input truncation.
163 @type recalc_tau: bool
164 @keyword r1_fit: A flag which if True will allow R1 values to be optimised. If False, preloaded R1 values will be used instead.
165 @type r1_fit: bool
166 """
167
168
169 if model not in MODEL_LIST_FULL:
170 raise RelaxError("The model '%s' is unknown." % model)
171 if values == None:
172 raise RelaxError("No values have been supplied to the target function.")
173 if errors == None:
174 raise RelaxError("No errors have been supplied to the target function.")
175 if missing == None:
176 raise RelaxError("No missing data information has been supplied to the target function.")
177 if model in MODEL_LIST_R1RHO_OFF_RES:
178 if chemical_shifts == None:
179 raise RelaxError("Chemical shifts must be supplied for the '%s' R1rho off-resonance dispersion model." % model)
180 if not r1_fit and r1 is None:
181 raise RelaxError("R1 relaxation rates must be supplied for the '%s' R1rho off-resonance dispersion model when not fitting the values." % model)
182
183
184 self.model = model
185 self.num_params = num_params
186 self.exp_types = exp_types
187 self.scaling_matrix = scaling_matrix
188 self.values_orig = values
189 self.cpmg_frqs_orig = cpmg_frqs
190 self.spin_lock_nu1_orig = spin_lock_nu1
191
192
193
194
195 self.NE = len(self.exp_types)
196 self.NS = num_spins
197 self.NM = num_frq
198
199
200 max_NO = 1
201 for ei in range(self.NE):
202 for si in range(self.NS):
203 for mi in range(self.NM):
204 nr_NO = len(offset[ei][si][mi])
205 if nr_NO > max_NO:
206 max_NO = nr_NO
207
208
209 self.NO = max_NO
210
211
212 max_ND = 1
213 for ei in range(self.NE):
214 for si in range(self.NS):
215 for mi in range(self.NM):
216 for oi in range(self.NO):
217 if cpmg_frqs != None and len(cpmg_frqs[ei][mi][oi]):
218 nr_ND = len(cpmg_frqs[ei][mi][oi])
219 if nr_ND > max_ND:
220 max_ND = nr_ND
221 elif spin_lock_nu1 != None and len(spin_lock_nu1[ei][mi][oi]):
222 nr_ND = len(spin_lock_nu1[ei][mi][oi])
223 if nr_ND > max_ND:
224 max_ND = nr_ND
225
226
227 self.NO = max_NO
228 self.ND = max_ND
229
230
231 self.numpy_array_shape = [self.NE, self.NS, self.NM, self.NO, self.ND]
232
233
234 numpy_array_zeros = zeros(self.numpy_array_shape, float64)
235 numpy_array_ones = ones(self.numpy_array_shape, float64)
236
237
238 self.no_nd_ones = ones([self.NO, self.ND], float64)
239 self.nm_no_nd_ones = ones([self.NM, self.NO, self.ND], float64)
240
241
242 self.r1_struct = deepcopy(numpy_array_zeros)
243 self.r1rho_prime_struct = deepcopy(numpy_array_zeros)
244 self.r20_struct = deepcopy(numpy_array_zeros)
245 self.r20a_struct = deepcopy(numpy_array_zeros)
246 self.r20b_struct = deepcopy(numpy_array_zeros)
247 self.r20c_struct = deepcopy(numpy_array_zeros)
248
249
250 self.dw_struct = deepcopy(numpy_array_zeros)
251 self.dwH_struct = deepcopy(numpy_array_zeros)
252 self.dw_AB_struct = deepcopy(numpy_array_zeros)
253 self.dw_BC_struct = deepcopy(numpy_array_zeros)
254 self.dwH_AB_struct = deepcopy(numpy_array_zeros)
255 self.dwH_BC_struct = deepcopy(numpy_array_zeros)
256 self.phi_ex_struct = deepcopy(numpy_array_zeros)
257 self.phi_ex_B_struct = deepcopy(numpy_array_zeros)
258 self.phi_ex_C_struct = deepcopy(numpy_array_zeros)
259
260
261 self.values = deepcopy(numpy_array_zeros)
262 self.errors = deepcopy(numpy_array_ones)
263 self.missing = deepcopy(numpy_array_zeros)
264 self.disp_struct = deepcopy(numpy_array_zeros)
265
266
267 self.cpmg_frqs = deepcopy(numpy_array_ones)
268 self.frqs = deepcopy(numpy_array_zeros)
269 self.frqs_squared = deepcopy(numpy_array_zeros)
270 self.frqs_H = deepcopy(numpy_array_zeros)
271 self.relax_times = deepcopy(numpy_array_zeros)
272 if model in MODEL_LIST_INV_RELAX_TIMES:
273 self.inv_relax_times = deepcopy(numpy_array_zeros)
274 self.tau_cpmg = deepcopy(numpy_array_zeros)
275 self.power = deepcopy(numpy_array_zeros)
276 self.r1 = deepcopy(numpy_array_zeros)
277 self.spin_lock_omega1 = deepcopy(numpy_array_zeros)
278 self.spin_lock_omega1_squared = deepcopy(numpy_array_zeros)
279 self.offset = deepcopy(numpy_array_zeros)
280 self.chemical_shifts = deepcopy(numpy_array_zeros)
281 self.tilt_angles = deepcopy(numpy_array_zeros)
282 self.num_offsets = zeros([self.NE, self.NS, self.NM], int16)
283 self.num_disp_points = zeros([self.NE, self.NS, self.NM, self.NO], int16)
284
285
286 self.has_missing = False
287
288
289 for ei in range(self.NE):
290 for si in range(self.NS):
291 for mi in range(self.NM):
292
293 frq = frqs[ei][si][mi]
294 self.frqs[ei, si, mi, :] = frq
295 self.frqs_squared[ei, si, mi, :] = frq**2
296 if frqs_H != None:
297 frq_H = frqs_H[ei][si][mi]
298 self.frqs_H[ei, si, mi, :] = frq_H
299
300
301 if r1 is not None:
302 r1_l = r1[si][mi]
303 self.r1[ei, si, mi, :] = r1_l
304
305
306 if chemical_shifts != None:
307 chemical_shift = chemical_shifts[ei][si][mi]
308 self.chemical_shifts[ei, si, mi, :] = chemical_shift
309
310
311 if len(offset[ei][si][mi]):
312 self.num_offsets[ei, si, mi] = len(self.offset[ei, si, mi])
313 else:
314 self.num_offsets[ei, si, mi] = 0
315
316
317 for oi in range(self.NO):
318 if cpmg_frqs != None and len(cpmg_frqs[ei][mi][oi]):
319 cpmg_frqs_list = cpmg_frqs[ei][mi][oi]
320 num_disp_points = len(cpmg_frqs_list)
321 self.cpmg_frqs[ei, si, mi, oi, :num_disp_points] = cpmg_frqs_list
322
323 for di in range(num_disp_points):
324 cpmg_frq = cpmg_frqs[ei][mi][oi][di]
325
326
327 relax_time = max(relax_times[ei][mi][oi][di])
328
329 if isNaN(relax_time):
330 power = 0
331
332
333 else:
334 power = int(round(cpmg_frq * relax_time))
335 self.power[ei, si, mi, oi, di] = power
336
337
338 if recalc_tau:
339 tau_cpmg = 0.25 * relax_time / power
340 else:
341 tau_cpmg = 0.25 / cpmg_frq
342 self.tau_cpmg[ei, si, mi, oi, di] = tau_cpmg
343
344 elif spin_lock_nu1 != None and len(spin_lock_nu1[ei][mi][oi]):
345 num_disp_points = len( spin_lock_nu1[ei][mi][oi] )
346 else:
347 num_disp_points = 0
348
349 self.num_disp_points[ei, si, mi, oi] = num_disp_points
350
351
352 self.values[ei, si, mi, oi, :num_disp_points] = values[ei][si][mi][oi]
353 self.errors[ei, si, mi, oi, :num_disp_points] = errors[ei][si][mi][oi]
354 self.disp_struct[ei, si, mi, oi, :num_disp_points] = ones(num_disp_points)
355
356
357 if offset != None and len(offset[ei][si][mi]):
358 self.offset[ei, si, mi, oi] = offset[ei][si][mi][oi]
359
360
361 for di in range(num_disp_points):
362 if missing[ei][si][mi][oi][di]:
363 self.has_missing = True
364 self.missing[ei, si, mi, oi, di] = 1.0
365
366
367 if tilt_angles != None and di < len(tilt_angles[ei][si][mi][oi]):
368 self.tilt_angles[ei, si, mi, oi, di] = tilt_angles[ei][si][mi][oi][di]
369
370
371 if spin_lock_nu1 != None and len(spin_lock_nu1[ei][mi][oi]):
372 self.spin_lock_omega1[ei, si, mi, oi, di] = 2.0 * pi * spin_lock_nu1[ei][mi][oi][di]
373 self.spin_lock_omega1_squared[ei, si, mi, oi, di] = self.spin_lock_omega1[ei, si, mi, oi, di] ** 2
374
375
376
377 if relax_times != None and len(relax_times[ei][mi][oi]):
378 relax_time = max(relax_times[ei][mi][oi][di])
379 self.relax_times[ei, si, mi, oi, di] = relax_time
380
381
382 if model in MODEL_LIST_INV_RELAX_TIMES:
383 self.inv_relax_times[ei, si, mi, oi, di] = 1.0 / relax_time
384
385
386 if model in MODEL_LIST_INV_RELAX_TIMES:
387
388 if sum(self.inv_relax_times) == 0.0:
389 raise RelaxError("The inverted relaxation time data array all contains 0.0. Please check your setup.")
390
391
392 if not all(isfinite(self.inv_relax_times)):
393 raise RelaxError("The inverted relaxation time data array contains not finite values. Please check your setup.")
394
395
396 self.back_calc = deepcopy(self.values)
397
398
399
400 self.mask_replace_blank = masked_equal(self.missing, 1.0)
401
402
403 self.experiment_type_setup()
404
405
406 self.scaling_flag = False
407 if self.scaling_matrix is not None:
408 self.scaling_flag = True
409
410
411 self.end_index = []
412
413
414 if r1_fit:
415
416 self.end_index.append(self.NE * self.NS * self.NM)
417
418 self.end_index.append(self.end_index[-1] + self.NE * self.NS * self.NM)
419
420
421 else:
422
423 self.end_index.append(self.NE * self.NS * self.NM)
424
425 if model in MODEL_LIST_R20B:
426 self.end_index.append(2 * self.NE * self.NS * self.NM)
427
428
429 self.end_index.append(self.end_index[-1] + self.NS)
430
431
432 if model in MODEL_LIST_DW_MIX_DOUBLE:
433 self.end_index.append(self.end_index[-1] + self.NS)
434
435 elif model in MODEL_LIST_DW_MIX_QUADRUPLE:
436 self.end_index.append(self.end_index[-1] + self.NS)
437 self.end_index.append(self.end_index[-1] + self.NS)
438 self.end_index.append(self.end_index[-1] + self.NS)
439
440
441 if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL]:
442 self.r180x = r180x_3d()
443
444
445 if model in [MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_MMQ_2SITE]:
446 self.M0 = zeros(2, float64)
447
448 if model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]:
449 self.M0 = zeros(3, float64)
450
451 if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL]:
452 M0_0 = zeros( [self.NE, self.NS, self.NM, self.NO, self.ND, 7, 1], float64)
453 M0_0[:, :, :, :, :, 0, 0] = 0.5
454 self.M0 = M0_0
455
456
457 self.M0_T = rollaxis(self.M0, 6, 5)
458
459 if model in [MODEL_NS_R1RHO_2SITE]:
460
461 da_mat = self.chemical_shifts - self.offset
462
463
464 theta_mat = arctan2(self.spin_lock_omega1, da_mat)
465 M0_0 = zeros([6, 1], float64)
466 M0_0[0, 0] = 1
467 M0_sin = multiply.outer( sin(theta_mat), M0_0 )
468 M0_2 = zeros([6, 1], float64)
469 M0_2[2, 0] = 1
470 M0_cos = multiply.outer( cos(theta_mat), M0_2 )
471 self.M0 = M0_sin + M0_cos
472
473
474 self.M0_T = rollaxis(self.M0, 6, 5)
475
476 if model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]:
477 self.M0 = zeros(9, float64)
478
479
480 da_mat = self.chemical_shifts - self.offset
481
482
483 theta_mat = arctan2(self.spin_lock_omega1, da_mat)
484 M0_0 = zeros([9, 1], float64)
485 M0_0[0, 0] = 1
486
487
488 M0_sin = multiply.outer( sin(theta_mat), M0_0 )
489 M0_2 = zeros([9, 1], float64)
490 M0_2[2, 0] = 1
491
492
493 M0_cos = multiply.outer( cos(theta_mat), M0_2 )
494 self.M0 = M0_sin + M0_cos
495
496
497 self.M0_T = rollaxis(self.M0, 6, 5)
498
499
500 if model == MODEL_NOREX:
501
502 if self.exp_types[0] in EXP_TYPE_LIST_CPMG:
503 self.func = self.func_NOREX
504 else:
505 if r1_fit:
506 self.func = self.func_NOREX_R1RHO_FIT_R1
507 else:
508 self.func = self.func_NOREX_R1RHO
509 if model == MODEL_LM63:
510 self.func = self.func_LM63
511 if model == MODEL_LM63_3SITE:
512 self.func = self.func_LM63_3site
513 if model == MODEL_CR72_FULL:
514 self.func = self.func_CR72_full
515 if model == MODEL_CR72:
516 self.func = self.func_CR72
517 if model == MODEL_IT99:
518 self.func = self.func_IT99
519 if model == MODEL_TSMFK01:
520 self.func = self.func_TSMFK01
521 if model == MODEL_B14:
522 self.func = self.func_B14
523 if model == MODEL_B14_FULL:
524 self.func = self.func_B14_full
525 if model == MODEL_NS_CPMG_2SITE_3D_FULL:
526 self.func = self.func_ns_cpmg_2site_3D_full
527 if model == MODEL_NS_CPMG_2SITE_3D:
528 self.func = self.func_ns_cpmg_2site_3D
529 if model == MODEL_NS_CPMG_2SITE_EXPANDED:
530 self.func = self.func_ns_cpmg_2site_expanded
531 if model == MODEL_NS_CPMG_2SITE_STAR_FULL:
532 self.func = self.func_ns_cpmg_2site_star_full
533 if model == MODEL_NS_CPMG_2SITE_STAR:
534 self.func = self.func_ns_cpmg_2site_star
535 if model == MODEL_M61:
536 self.func = self.func_M61
537 if model == MODEL_M61B:
538 self.func = self.func_M61b
539 if model == MODEL_DPL94:
540 if r1_fit:
541 self.func = self.func_DPL94_fit_r1
542 else:
543 self.func = self.func_DPL94
544 if model == MODEL_TP02:
545 if r1_fit:
546 self.func = self.func_TP02_fit_r1
547 else:
548 self.func = self.func_TP02
549 if model == MODEL_TAP03:
550 if r1_fit:
551 self.func = self.func_TAP03_fit_r1
552 else:
553 self.func = self.func_TAP03
554 if model == MODEL_MP05:
555 if r1_fit:
556 self.func = self.func_MP05_fit_r1
557 else:
558 self.func = self.func_MP05
559 if model == MODEL_NS_R1RHO_2SITE:
560 if r1_fit:
561 self.func = self.func_ns_r1rho_2site_fit_r1
562 else:
563 self.func = self.func_ns_r1rho_2site
564 if model == MODEL_NS_R1RHO_3SITE:
565 self.func = self.func_ns_r1rho_3site
566 if model == MODEL_NS_R1RHO_3SITE_LINEAR:
567 self.func = self.func_ns_r1rho_3site_linear
568 if model == MODEL_MMQ_CR72:
569 self.func = self.func_mmq_CR72
570 if model == MODEL_NS_MMQ_2SITE:
571 self.func = self.func_ns_mmq_2site
572 if model == MODEL_NS_MMQ_3SITE:
573 self.func = self.func_ns_mmq_3site
574 if model == MODEL_NS_MMQ_3SITE_LINEAR:
575 self.func = self.func_ns_mmq_3site_linear
576
577
578 - def calc_B14_chi2(self, R20A=None, R20B=None, dw=None, pA=None, kex=None):
579 """Calculate the chi-squared value of the Baldwin (2014) 2-site exact solution model for all time scales.
580
581
582 @keyword R20A: The R2 value for state A in the absence of exchange.
583 @type R20A: list of float
584 @keyword R20B: The R2 value for state B in the absence of exchange.
585 @type R20B: list of float
586 @keyword dw: The chemical shift differences in ppm for each spin.
587 @type dw: list of float
588 @keyword pA: The population of state A.
589 @type pA: float
590 @keyword kex: The rate of exchange.
591 @type kex: float
592 @return: The chi-squared value.
593 @rtype: float
594 """
595
596
597 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
598
599
600 self.r20a_struct[:] = multiply.outer( R20A.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
601 self.r20b_struct[:] = multiply.outer( R20B.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
602
603
604 r2eff_B14(r20a=self.r20a_struct, r20b=self.r20b_struct, pA=pA, dw=self.dw_struct, dw_orig=dw, kex=kex, ncyc=self.power, inv_tcpmg=self.inv_relax_times, tcp=self.tau_cpmg, back_calc=self.back_calc)
605
606
607 self.back_calc = self.back_calc*self.disp_struct
608
609
610 if self.has_missing:
611
612 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
613
614
615 return chi2_rankN(self.values, self.back_calc, self.errors)
616
617
618 - def calc_CR72_chi2(self, R20A=None, R20B=None, dw=None, pA=None, kex=None):
619 """Calculate the chi-squared value of the Carver and Richards (1972) 2-site exchange model on all time scales.
620
621 @keyword R20A: The R2 value for state A in the absence of exchange.
622 @type R20A: list of float
623 @keyword R20B: The R2 value for state B in the absence of exchange.
624 @type R20B: list of float
625 @keyword dw: The chemical shift differences in ppm for each spin.
626 @type dw: list of float
627 @keyword pA: The population of state A.
628 @type pA: float
629 @keyword kex: The rate of exchange.
630 @type kex: float
631 @return: The chi-squared value.
632 @rtype: float
633 """
634
635
636 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
637
638
639 self.r20a_struct[:] = multiply.outer( R20A.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
640 self.r20b_struct[:] = multiply.outer( R20B.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
641
642
643 r2eff_CR72(r20a=self.r20a_struct, r20a_orig=R20A, r20b=self.r20b_struct, r20b_orig=R20B, pA=pA, dw=self.dw_struct, dw_orig=dw, kex=kex, cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc)
644
645
646 self.back_calc = self.back_calc*self.disp_struct
647
648
649 if self.has_missing:
650
651 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
652
653
654 return chi2_rankN(self.values, self.back_calc, self.errors)
655
656
657 - def calc_DPL94(self, R1=None, r1rho_prime=None, phi_ex=None, kex=None):
658 """Calculation function for the Davis, Perlman and London (1994) fast 2-site off-resonance exchange model for R1rho-type experiments.
659
660 @keyword R1: The R1 value.
661 @type R1: list of float
662 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
663 @type r1rho_prime: list of float
664 @keyword phi_ex: The fast exchange factor pA.pB.dw**2 (ppm).
665 @type phi_ex: list of float
666 @keyword kex: The rate of exchange.
667 @type kex: float
668 @return: The chi-squared value.
669 @rtype: float
670 """
671
672
673 multiply( multiply.outer( phi_ex.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_squared, out=self.phi_ex_struct )
674
675
676 self.r1rho_prime_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
677
678
679 r1rho_DPL94(r1rho_prime=self.r1rho_prime_struct, phi_ex=self.phi_ex_struct, kex=kex, theta=self.tilt_angles, R1=R1, spin_lock_fields2=self.spin_lock_omega1_squared, back_calc=self.back_calc)
680
681
682 self.back_calc = self.back_calc*self.disp_struct
683
684
685 if self.has_missing:
686
687 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
688
689
690 return chi2_rankN(self.values, self.back_calc, self.errors)
691
692
693 - def calc_MP05(self, R1=None, r1rho_prime=None, dw=None, pA=None, kex=None):
694 """Calculation function for the Miloushev and Palmer (2005) R1rho off-resonance 2-site model.
695
696 @keyword R1: The R1 value.
697 @type R1: list of float
698 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
699 @type r1rho_prime: list of float
700 @keyword dw: The chemical shift differences in ppm for each spin.
701 @type dw: list of float
702 @keyword pA: The population of state A.
703 @type pA: float
704 @keyword kex: The rate of exchange.
705 @type kex: float
706 @return: The chi-squared value.
707 @rtype: float
708 """
709
710
711 self.r1rho_prime_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
712
713
714 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
715
716
717 r1rho_MP05(r1rho_prime=self.r1rho_prime_struct, omega=self.chemical_shifts, offset=self.offset, pA=pA, dw=self.dw_struct, kex=kex, R1=R1, spin_lock_fields=self.spin_lock_omega1, spin_lock_fields2=self.spin_lock_omega1_squared, back_calc=self.back_calc)
718
719
720 self.back_calc = self.back_calc*self.disp_struct
721
722
723 if self.has_missing:
724
725 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
726
727
728 return chi2_rankN(self.values, self.back_calc, self.errors)
729
730
732 """Calculation function for no exchange, for R1rho off resonance models.
733
734 @keyword R1: The R1 value.
735 @type R1: list of float
736 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
737 @type r1rho_prime: list of float
738 @return: The chi-squared value.
739 @rtype: float
740 """
741
742
743 self.r1rho_prime_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
744
745
746 self.back_calc[:] = R1 * cos(self.tilt_angles)**2 + self.r1rho_prime_struct * sin(self.tilt_angles)**2
747
748
749 self.back_calc = self.back_calc*self.disp_struct
750
751
752 if self.has_missing:
753
754 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
755
756
757 return chi2_rankN(self.values, self.back_calc, self.errors)
758
759
761 """Calculate the chi-squared value of the 'NS CPMG 2-site' models.
762
763 @keyword R20A: The R2 value for state A in the absence of exchange.
764 @type R20A: list of float
765 @keyword R20B: The R2 value for state B in the absence of exchange.
766 @type R20B: list of float
767 @keyword dw: The chemical shift differences in ppm for each spin.
768 @type dw: list of float
769 @keyword pA: The population of state A.
770 @type pA: float
771 @keyword kex: The rate of exchange.
772 @type kex: float
773 @return: The chi-squared value.
774 @rtype: float
775 """
776
777
778 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
779
780
781 self.r20a_struct[:] = multiply.outer( R20A.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
782 self.r20b_struct[:] = multiply.outer( R20B.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
783
784
785 r2eff_ns_cpmg_2site_3D(r180x=self.r180x, M0=self.M0, M0_T=self.M0_T, r20a=self.r20a_struct, r20b=self.r20b_struct, pA=pA, dw=self.dw_struct, dw_orig=dw, kex=kex, inv_tcpmg=self.inv_relax_times, tcp=self.tau_cpmg, back_calc=self.back_calc, num_points=self.num_disp_points, power=self.power)
786
787
788 self.back_calc = self.back_calc*self.disp_struct
789
790
791 if self.has_missing:
792
793 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
794
795
796 return chi2_rankN(self.values, self.back_calc, self.errors)
797
798
800 """Calculate the chi-squared value of the 'NS CPMG 2-site star' models.
801
802 @keyword R20A: The R2 value for state A in the absence of exchange.
803 @type R20A: list of float
804 @keyword R20B: The R2 value for state B in the absence of exchange.
805 @type R20B: list of float
806 @keyword dw: The chemical shift differences in ppm for each spin.
807 @type dw: list of float
808 @keyword pA: The population of state A.
809 @type pA: float
810 @keyword kex: The rate of exchange.
811 @type kex: float
812 @return: The chi-squared value.
813 @rtype: float
814 """
815
816
817 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
818
819
820 self.r20a_struct[:] = multiply.outer( R20A.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
821 self.r20b_struct[:] = multiply.outer( R20B.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
822
823
824 r2eff_ns_cpmg_2site_star(M0=self.M0, r20a=self.r20a_struct, r20b=self.r20b_struct, pA=pA, dw=self.dw_struct, dw_orig=dw, kex=kex, inv_tcpmg=self.inv_relax_times, tcp=self.tau_cpmg, back_calc=self.back_calc, num_points=self.num_disp_points, power=self.power)
825
826
827 self.back_calc = self.back_calc*self.disp_struct
828
829
830 if self.has_missing:
831
832 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
833
834
835 return chi2_rankN(self.values, self.back_calc, self.errors)
836
837
838 - def calc_ns_mmq_3site_chi2(self, R20A=None, R20B=None, R20C=None, dw_AB=None, dw_BC=None, dwH_AB=None, dwH_BC=None, pA=None, pB=None, kex_AB=None, kex_BC=None, kex_AC=None):
839 """Calculate the chi-squared value for the 'NS MMQ 3-site' models.
840
841 @keyword R20A: The R2 value for state A in the absence of exchange.
842 @type R20A: list of float
843 @keyword R20B: The R2 value for state B in the absence of exchange.
844 @type R20B: list of float
845 @keyword R20C: The R2 value for state C in the absence of exchange.
846 @type R20C: list of float
847 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s.
848 @type dw_AB: float
849 @keyword dw_BC: The chemical exchange difference between states B and C in rad/s.
850 @type dw_BC: float
851 @keyword dwH_AB: The proton chemical exchange difference between states A and B in rad/s.
852 @type dwH_AB: float
853 @keyword dwH_BC: The proton chemical exchange difference between states B and C in rad/s.
854 @type dwH_BC: float
855 @keyword pA: The population of state A.
856 @type pA: float
857 @keyword kex_AB: The rate of exchange between states A and B.
858 @type kex_AB: float
859 @keyword kex_BC: The rate of exchange between states B and C.
860 @type kex_BC: float
861 @keyword kex_AC: The rate of exchange between states A and C.
862 @type kex_AC: float
863 @return: The chi-squared value.
864 @rtype: float
865 """
866
867
868 multiply( multiply.outer( dw_AB.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_AB_struct )
869 multiply( multiply.outer( dw_BC.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_BC_struct )
870 multiply( multiply.outer( dwH_AB.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_H, out=self.dwH_AB_struct )
871 multiply( multiply.outer( dwH_BC.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_H, out=self.dwH_BC_struct )
872
873
874 self.r20a_struct[:] = multiply.outer( R20A.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
875 self.r20b_struct[:] = multiply.outer( R20B.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
876 self.r20c_struct[:] = multiply.outer( R20C.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
877
878
879 for ei in range(self.NE):
880 r20a = self.r20a_struct[ei]
881 r20b = self.r20b_struct[ei]
882 r20c = self.r20b_struct[ei]
883 dw_AB_frq = self.dw_AB_struct[ei]
884 dw_BC_frq = self.dw_BC_struct[ei]
885 dwH_AB_frq = self.dwH_AB_struct[ei]
886 dwH_BC_frq = self.dwH_BC_struct[ei]
887
888
889 aliased_dwH_AB = 0.0 * self.dwH_AB_struct[ei]
890 aliased_dwH_BC = 0.0 * self.dwH_BC_struct[ei]
891 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ:
892 aliased_dw_AB = dw_AB_frq
893 aliased_dw_BC = dw_BC_frq
894 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ:
895 aliased_dw_AB = dwH_AB_frq
896 aliased_dw_BC = dwH_BC_frq
897 elif self.exp_types[ei] == EXP_TYPE_CPMG_DQ:
898 aliased_dw_AB = dw_AB_frq + dwH_AB_frq
899 aliased_dw_BC = dw_BC_frq + dwH_BC_frq
900 elif self.exp_types[ei] == EXP_TYPE_CPMG_ZQ:
901 aliased_dw_AB = dw_AB_frq - dwH_AB_frq
902 aliased_dw_BC = dw_BC_frq - dwH_BC_frq
903 elif self.exp_types[ei] == EXP_TYPE_CPMG_MQ:
904 aliased_dw_AB = dw_AB_frq
905 aliased_dw_BC = dw_BC_frq
906 aliased_dwH_AB = dwH_AB_frq
907 aliased_dwH_BC = dwH_BC_frq
908 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_MQ:
909 aliased_dw_AB = dwH_AB_frq
910 aliased_dw_BC = dwH_BC_frq
911 aliased_dwH_AB = dw_AB_frq
912 aliased_dwH_BC = dw_BC_frq
913
914
915 self.r2eff_ns_mmq[ei](M0=self.M0, R20A=r20a, R20B=r20b, R20C=r20c, pA=pA, pB=pB, dw_AB=aliased_dw_AB, dw_BC=aliased_dw_BC, dwH_AB=aliased_dwH_AB, dwH_BC=aliased_dwH_BC, kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=kex_AC, inv_tcpmg=self.inv_relax_times[ei], tcp=self.tau_cpmg[ei], back_calc=self.back_calc[ei], num_points=self.num_disp_points[ei], power=self.power[ei])
916
917
918 self.back_calc = self.back_calc*self.disp_struct
919
920
921 if self.has_missing:
922
923 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
924
925
926 return chi2_rankN(self.values, self.back_calc, self.errors)
927
928
930 """Calculation function for the reduced numerical solution for the 2-site Bloch-McConnell equations for R1rho data.
931
932 @keyword R1: The R1 value.
933 @type R1: list of float
934 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
935 @type r1rho_prime: list of float
936 @keyword dw: The chemical shift differences in ppm for each spin.
937 @type dw: list of float
938 @keyword pA: The population of state A.
939 @type pA: float
940 @keyword kex: The rate of exchange.
941 @type kex: float
942 @return: The chi-squared value.
943 @rtype: float
944 """
945
946
947 self.r1rho_prime_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
948
949
950 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
951
952
953 ns_r1rho_2site(M0=self.M0, M0_T=self.M0_T, r1rho_prime=self.r1rho_prime_struct, omega=self.chemical_shifts, offset=self.offset, r1=R1, pA=pA, dw=self.dw_struct, kex=kex, spin_lock_fields=self.spin_lock_omega1, relax_time=self.relax_times, inv_relax_time=self.inv_relax_times, back_calc=self.back_calc)
954
955
956 self.back_calc = self.back_calc*self.disp_struct
957
958
959 if self.has_missing:
960
961 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
962
963
964 return chi2_rankN(self.values, self.back_calc, self.errors)
965
966
967 - def calc_ns_r1rho_3site_chi2(self, r1rho_prime=None, dw_AB=None, dw_BC=None, pA=None, pB=None, kex_AB=None, kex_BC=None, kex_AC=None):
968 """Calculate the chi-squared value for the 'NS MMQ 3-site' models.
969
970 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
971 @type r1rho_prime: list of float
972 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s.
973 @type dw_AB: float
974 @keyword dw_BC: The chemical exchange difference between states B and C in rad/s.
975 @type dw_BC: float
976 @keyword pA: The population of state A.
977 @type pA: float
978 @keyword kex_AB: The rate of exchange between states A and B.
979 @type kex_AB: float
980 @keyword kex_BC: The rate of exchange between states B and C.
981 @type kex_BC: float
982 @keyword kex_AC: The rate of exchange between states A and C.
983 @type kex_AC: float
984 @return: The chi-squared value.
985 @rtype: float
986 """
987
988
989 multiply( multiply.outer( dw_AB.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_AB_struct )
990 multiply( multiply.outer( dw_BC.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_BC_struct )
991
992
993 self.r20_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
994
995
996 ns_r1rho_3site(M0=self.M0, M0_T=self.M0_T, r1rho_prime=self.r20_struct, omega=self.chemical_shifts, offset=self.offset, r1=self.r1, pA=pA, pB=pB, dw_AB=self.dw_AB_struct, dw_BC=self.dw_BC_struct, kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=kex_AC, spin_lock_fields=self.spin_lock_omega1, relax_time=self.relax_times, inv_relax_time=self.inv_relax_times, back_calc=self.back_calc, num_points=self.num_disp_points)
997
998
999 self.back_calc = self.back_calc*self.disp_struct
1000
1001
1002 if self.has_missing:
1003
1004 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1005
1006
1007 return chi2_rankN(self.values, self.back_calc, self.errors)
1008
1009
1010 - def calc_TAP03(self, R1=None, r1rho_prime=None, dw=None, pA=None, kex=None):
1011 """Calculation function for the Trott, Abergel and Palmer (2003) R1rho off-resonance 2-site model.
1012
1013 @keyword R1: The R1 value.
1014 @type R1: list of float
1015 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
1016 @type r1rho_prime: list of float
1017 @keyword dw: The chemical shift differences in ppm for each spin.
1018 @type dw: list of float
1019 @keyword pA: The population of state A.
1020 @type pA: float
1021 @keyword kex: The rate of exchange.
1022 @type kex: float
1023 @return: The chi-squared value.
1024 @rtype: float
1025 """
1026
1027
1028 self.r1rho_prime_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1029
1030
1031 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1032
1033
1034 r1rho_TAP03(r1rho_prime=self.r1rho_prime_struct, omega=self.chemical_shifts, offset=self.offset, pA=pA, dw=self.dw_struct, kex=kex, R1=R1, spin_lock_fields=self.spin_lock_omega1, spin_lock_fields2=self.spin_lock_omega1_squared, back_calc=self.back_calc)
1035
1036
1037 self.back_calc = self.back_calc*self.disp_struct
1038
1039
1040 if self.has_missing:
1041
1042 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1043
1044
1045 return chi2_rankN(self.values, self.back_calc, self.errors)
1046
1047
1048 - def calc_TP02(self, R1=None, r1rho_prime=None, dw=None, pA=None, kex=None):
1049 """Calculation function for the Trott and Palmer (2002) R1rho off-resonance 2-site model.
1050
1051 @keyword R1: The R1 value.
1052 @type R1: list of float
1053 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange.
1054 @type r1rho_prime: list of float
1055 @keyword dw: The chemical shift differences in ppm for each spin.
1056 @type dw: list of float
1057 @keyword pA: The population of state A.
1058 @type pA: float
1059 @keyword kex: The rate of exchange.
1060 @type kex: float
1061 @return: The chi-squared value.
1062 @rtype: float
1063 """
1064
1065
1066 self.r1rho_prime_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1067
1068
1069 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1070
1071
1072 r1rho_TP02(r1rho_prime=self.r1rho_prime_struct, omega=self.chemical_shifts, offset=self.offset, pA=pA, dw=self.dw_struct, kex=kex, R1=R1, spin_lock_fields=self.spin_lock_omega1, spin_lock_fields2=self.spin_lock_omega1_squared, back_calc=self.back_calc)
1073
1074
1075 self.back_calc = self.back_calc*self.disp_struct
1076
1077
1078 if self.has_missing:
1079
1080 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1081
1082
1083 return chi2_rankN(self.values, self.back_calc, self.errors)
1084
1085
1087 """Check the experiment types and simplify data structures.
1088
1089 For the single experiment type models, the first dimension of the values, errors, and missing data structures will be removed to simplify the target functions.
1090 """
1091
1092
1093 if self.model in MODEL_LIST_MMQ:
1094
1095 self.r2eff_ns_mmq = []
1096
1097
1098 for ei in range(self.NE):
1099
1100 if self.exp_types[ei] in [EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_DQ, EXP_TYPE_CPMG_ZQ]:
1101 if self.model == MODEL_NS_MMQ_2SITE:
1102 self.r2eff_ns_mmq.append(r2eff_ns_mmq_2site_sq_dq_zq)
1103 else:
1104 self.r2eff_ns_mmq.append(r2eff_ns_mmq_3site_sq_dq_zq)
1105
1106
1107 elif self.exp_types[ei] in [EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ]:
1108 if self.model == MODEL_NS_MMQ_2SITE:
1109 self.r2eff_ns_mmq.append(r2eff_ns_mmq_2site_mq)
1110 else:
1111 self.r2eff_ns_mmq.append(r2eff_ns_mmq_3site_mq)
1112
1113
1114 else:
1115
1116 if self.model != MODEL_NOREX and self.model in MODEL_LIST_CPMG and self.exp_types[0] != EXP_TYPE_CPMG_SQ:
1117 raise RelaxError("The '%s' CPMG model is not compatible with the '%s' experiment type." % (self.model, self.exp_types[0]))
1118 if self.model != MODEL_NOREX and self.model in MODEL_LIST_R1RHO and self.exp_types[0] != EXP_TYPE_R1RHO:
1119 raise RelaxError("The '%s' R1rho model is not compatible with the '%s' experiment type." % (self.model, self.exp_types[0]))
1120 if self.model != MODEL_NOREX and self.model in MODEL_LIST_MQ_CPMG and self.exp_types[0] != EXP_TYPE_CPMG_MQ:
1121 raise RelaxError("The '%s' CPMG model is not compatible with the '%s' experiment type." % (self.model, self.exp_types[0]))
1122
1123
1125 """Target function for the Baldwin (2014) 2-site exact solution model for all time scales, whereby the simplification R20A = R20B is assumed.
1126
1127 This assumes that pA > pB, and hence this must be implemented as a constraint.
1128
1129
1130 @param params: The vector of parameter values.
1131 @type params: numpy rank-1 float array
1132 @return: The chi-squared value.
1133 @rtype: float
1134 """
1135
1136
1137 if self.scaling_flag:
1138 params = dot(params, self.scaling_matrix)
1139
1140
1141 R20 = params[:self.end_index[0]]
1142 dw = params[self.end_index[0]:self.end_index[1]]
1143 pA = params[self.end_index[1]]
1144 kex = params[self.end_index[1]+1]
1145
1146
1147 return self.calc_B14_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)
1148
1149
1151 """Target function for the Baldwin (2014) 2-site exact solution model for all time scales.
1152
1153 This assumes that pA > pB, and hence this must be implemented as a constraint.
1154
1155
1156 @param params: The vector of parameter values.
1157 @type params: numpy rank-1 float array
1158 @return: The chi-squared value.
1159 @rtype: float
1160 """
1161
1162
1163 if self.scaling_flag:
1164 params = dot(params, self.scaling_matrix)
1165
1166
1167 R20 = params[:self.end_index[1]].reshape(self.NS*2, self.NM)
1168 R20A = R20[::2].flatten()
1169 R20B = R20[1::2].flatten()
1170 dw = params[self.end_index[1]:self.end_index[2]]
1171 pA = params[self.end_index[2]]
1172 kex = params[self.end_index[2]+1]
1173
1174
1175 return self.calc_B14_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)
1176
1177
1179 """Target function for the reduced Carver and Richards (1972) 2-site exchange model on all time scales.
1180
1181 This assumes that pA > pB, and hence this must be implemented as a constraint. For this model, the simplification R20A = R20B is assumed.
1182
1183
1184 @param params: The vector of parameter values.
1185 @type params: numpy rank-1 float array
1186 @return: The chi-squared value.
1187 @rtype: float
1188 """
1189
1190
1191 if self.scaling_flag:
1192 params = dot(params, self.scaling_matrix)
1193
1194
1195 R20 = params[:self.end_index[0]]
1196 dw = params[self.end_index[0]:self.end_index[1]]
1197 pA = params[self.end_index[1]]
1198 kex = params[self.end_index[1]+1]
1199
1200
1201 return self.calc_CR72_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)
1202
1203
1205 """Target function for the full Carver and Richards (1972) 2-site exchange model on all time scales.
1206
1207 This assumes that pA > pB, and hence this must be implemented as a constraint.
1208
1209
1210 @param params: The vector of parameter values.
1211 @type params: numpy rank-1 float array
1212 @return: The chi-squared value.
1213 @rtype: float
1214 """
1215
1216
1217 if self.scaling_flag:
1218 params = dot(params, self.scaling_matrix)
1219
1220
1221 R20 = params[:self.end_index[1]].reshape(self.NS*2, self.NM)
1222 R20A = R20[::2].flatten()
1223 R20B = R20[1::2].flatten()
1224 dw = params[self.end_index[1]:self.end_index[2]]
1225 pA = params[self.end_index[2]]
1226 kex = params[self.end_index[2]+1]
1227
1228
1229 return self.calc_CR72_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)
1230
1231
1233 """Target function for the Davis, Perlman and London (1994) fast 2-site off-resonance exchange model for R1rho-type experiments.
1234
1235 @param params: The vector of parameter values.
1236 @type params: numpy rank-1 float array
1237 @return: The chi-squared value.
1238 @rtype: float
1239 """
1240
1241
1242 if self.scaling_flag:
1243 params = dot(params, self.scaling_matrix)
1244
1245
1246 r1rho_prime = params[:self.end_index[0]]
1247 phi_ex = params[self.end_index[0]:self.end_index[1]]
1248 kex = params[self.end_index[1]]
1249
1250
1251 return self.calc_DPL94(R1=self.r1, r1rho_prime=r1rho_prime, phi_ex=phi_ex, kex=kex)
1252
1253
1255 """Target function for the Davis, Perlman and London (1994) fast 2-site off-resonance exchange model for R1rho-type experiments, whereby R1 is fitted.
1256
1257 @param params: The vector of parameter values.
1258 @type params: numpy rank-1 float array
1259 @return: The chi-squared value.
1260 @rtype: float
1261 """
1262
1263
1264 if self.scaling_flag:
1265 params = dot(params, self.scaling_matrix)
1266
1267
1268 r1 = params[:self.end_index[0]]
1269 r1rho_prime = params[self.end_index[0]:self.end_index[1]]
1270 phi_ex = params[self.end_index[1]:self.end_index[2]]
1271 kex = params[self.end_index[2]]
1272
1273
1274 self.r1_struct[:] = multiply.outer( r1.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1275
1276
1277 return self.calc_DPL94(R1=self.r1_struct, r1rho_prime=r1rho_prime, phi_ex=phi_ex, kex=kex)
1278
1279
1281 """Target function for the Ishima and Torchia (1999) 2-site model for all timescales with pA >> pB.
1282
1283 @param params: The vector of parameter values.
1284 @type params: numpy rank-1 float array
1285 @return: The chi-squared value.
1286 @rtype: float
1287 """
1288
1289
1290 if self.scaling_flag:
1291 params = dot(params, self.scaling_matrix)
1292
1293
1294 R20 = params[:self.end_index[0]]
1295 dw = params[self.end_index[0]:self.end_index[1]]
1296 pA = params[self.end_index[1]]
1297 tex = params[self.end_index[1]+1]
1298
1299
1300 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1301
1302
1303 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1304
1305
1306 r2eff_IT99(r20=self.r20_struct, pA=pA, dw=self.dw_struct, dw_orig=dw, tex=tex, cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc)
1307
1308
1309 self.back_calc = self.back_calc*self.disp_struct
1310
1311
1312 if self.has_missing:
1313
1314 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1315
1316
1317 return chi2_rankN(self.values, self.back_calc, self.errors)
1318
1319
1321 """Target function for the Luz and Meiboom (1963) fast 3-site exchange model.
1322
1323 @param params: The vector of parameter values.
1324 @type params: numpy rank-1 float array
1325 @return: The chi-squared value.
1326 @rtype: float
1327 """
1328
1329
1330 if self.scaling_flag:
1331 params = dot(params, self.scaling_matrix)
1332
1333
1334 R20 = params[:self.end_index[0]]
1335 phi_ex_B = params[self.end_index[0]:self.end_index[1]]
1336 phi_ex_C = params[self.end_index[1]:self.end_index[2]]
1337 kB = params[self.end_index[2]]
1338 kC = params[self.end_index[2]+1]
1339
1340
1341 multiply( multiply.outer( phi_ex_B.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_squared, out=self.phi_ex_B_struct )
1342 multiply( multiply.outer( phi_ex_C.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_squared, out=self.phi_ex_C_struct )
1343
1344
1345 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1346
1347
1348 r2eff_LM63_3site(r20=self.r20_struct, phi_ex_B=self.phi_ex_B_struct, phi_ex_C=self.phi_ex_C_struct, kB=kB, kC=kC, cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc)
1349
1350
1351 self.back_calc = self.back_calc*self.disp_struct
1352
1353
1354 if self.has_missing:
1355
1356 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1357
1358
1359 return chi2_rankN(self.values, self.back_calc, self.errors)
1360
1361
1363 """Target function for the Luz and Meiboom (1963) fast 2-site exchange model.
1364
1365 @param params: The vector of parameter values.
1366 @type params: numpy rank-1 float array
1367 @return: The chi-squared value.
1368 @rtype: float
1369 """
1370
1371
1372 if self.scaling_flag:
1373 params = dot(params, self.scaling_matrix)
1374
1375
1376 R20 = params[:self.end_index[0]]
1377 phi_ex = params[self.end_index[0]:self.end_index[1]]
1378 kex = params[self.end_index[1]]
1379
1380
1381 multiply( multiply.outer( phi_ex.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_squared, out=self.phi_ex_struct )
1382
1383
1384 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1385
1386
1387 r2eff_LM63(r20=self.r20_struct, phi_ex=self.phi_ex_struct, kex=kex, cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc)
1388
1389
1390 self.back_calc = self.back_calc*self.disp_struct
1391
1392
1393 if self.has_missing:
1394
1395 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1396
1397
1398 return chi2_rankN(self.values, self.back_calc, self.errors)
1399
1400
1402 """Target function for the Meiboom (1961) fast 2-site exchange model for R1rho-type experiments.
1403
1404 @param params: The vector of parameter values.
1405 @type params: numpy rank-1 float array
1406 @return: The chi-squared value.
1407 @rtype: float
1408 """
1409
1410
1411 if self.scaling_flag:
1412 params = dot(params, self.scaling_matrix)
1413
1414
1415 R20 = params[:self.end_index[0]]
1416 phi_ex = params[self.end_index[0]:self.end_index[1]]
1417 kex = params[self.end_index[1]]
1418
1419
1420 multiply( multiply.outer( phi_ex.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_squared, out=self.phi_ex_struct )
1421
1422
1423 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1424
1425
1426 r1rho_M61(r1rho_prime=self.r20_struct, phi_ex=self.phi_ex_struct, kex=kex, spin_lock_fields2=self.spin_lock_omega1_squared, back_calc=self.back_calc)
1427
1428
1429 self.back_calc = self.back_calc*self.disp_struct
1430
1431
1432 if self.has_missing:
1433
1434 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1435
1436
1437 return chi2_rankN(self.values, self.back_calc, self.errors)
1438
1439
1441 """Target function for the Meiboom (1961) R1rho on-resonance 2-site model for skewed populations (pA >> pB).
1442
1443 @param params: The vector of parameter values.
1444 @type params: numpy rank-1 float array
1445 @return: The chi-squared value.
1446 @rtype: float
1447 """
1448
1449
1450 if self.scaling_flag:
1451 params = dot(params, self.scaling_matrix)
1452
1453
1454 R20 = params[:self.end_index[0]]
1455 dw = params[self.end_index[0]:self.end_index[1]]
1456 pA = params[self.end_index[1]]
1457 kex = params[self.end_index[1]+1]
1458
1459
1460 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1461
1462
1463 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1464
1465
1466 r1rho_M61b(r1rho_prime=self.r20_struct, pA=pA, dw=self.dw_struct, kex=kex, spin_lock_fields2=self.spin_lock_omega1_squared, back_calc=self.back_calc)
1467
1468
1469 self.back_calc = self.back_calc*self.disp_struct
1470
1471
1472 if self.has_missing:
1473
1474 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1475
1476
1477 return chi2_rankN(self.values, self.back_calc, self.errors)
1478
1479
1481 """Target function for the Miloushev and Palmer (2005) R1rho off-resonance 2-site model.
1482
1483 @param params: The vector of parameter values.
1484 @type params: numpy rank-1 float array
1485 @return: The chi-squared value.
1486 @rtype: float
1487 """
1488
1489
1490 if self.scaling_flag:
1491 params = dot(params, self.scaling_matrix)
1492
1493
1494 r1rho_prime = params[:self.end_index[0]]
1495 dw = params[self.end_index[0]:self.end_index[1]]
1496 pA = params[self.end_index[1]]
1497 kex = params[self.end_index[1]+1]
1498
1499
1500 return self.calc_MP05(R1=self.r1, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
1501
1502
1504 """Target function for the Miloushev and Palmer (2005) R1rho off-resonance 2-site model, whereby R1 is fitted.
1505
1506 @param params: The vector of parameter values.
1507 @type params: numpy rank-1 float array
1508 @return: The chi-squared value.
1509 @rtype: float
1510 """
1511
1512
1513 if self.scaling_flag:
1514 params = dot(params, self.scaling_matrix)
1515
1516
1517 r1 = params[:self.end_index[0]]
1518 r1rho_prime = params[self.end_index[0]:self.end_index[1]]
1519 dw = params[self.end_index[1]:self.end_index[2]]
1520 pA = params[self.end_index[2]]
1521 kex = params[self.end_index[2]+1]
1522
1523
1524 self.r1_struct[:] = multiply.outer( r1.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1525
1526
1527 return self.calc_MP05(R1=self.r1_struct, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
1528
1529
1531 """Target function for the CR72 model extended for MQ CPMG data.
1532
1533 @param params: The vector of parameter values.
1534 @type params: numpy rank-1 float array
1535 @return: The chi-squared value.
1536 @rtype: float
1537 """
1538
1539
1540 if self.scaling_flag:
1541 params = dot(params, self.scaling_matrix)
1542
1543
1544 R20 = params[:self.end_index[0]]
1545 dw = params[self.end_index[0]:self.end_index[1]]
1546 dwH = params[self.end_index[1]:self.end_index[2]]
1547 pA = params[self.end_index[2]]
1548 kex = params[self.end_index[2]+1]
1549
1550
1551 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1552 multiply( multiply.outer( dwH.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_H, out=self.dwH_struct )
1553
1554
1555 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1556
1557
1558 for ei in range(self.NE):
1559 r20 = self.r20_struct[ei]
1560 dw_frq = self.dw_struct[ei]
1561 dwH_frq = self.dwH_struct[ei]
1562
1563
1564 aliased_dwH = 0.0
1565 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ:
1566 aliased_dw = dw_frq
1567 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ:
1568 aliased_dw = dwH_frq
1569 elif self.exp_types[ei] == EXP_TYPE_CPMG_DQ:
1570 aliased_dw = dw_frq + dwH_frq
1571 elif self.exp_types[ei] == EXP_TYPE_CPMG_ZQ:
1572 aliased_dw = dw_frq - dwH_frq
1573 elif self.exp_types[ei] == EXP_TYPE_CPMG_MQ:
1574 aliased_dw = dw_frq
1575 aliased_dwH = dwH_frq
1576 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_MQ:
1577 aliased_dw = dwH_frq
1578 aliased_dwH = dw_frq
1579
1580
1581 r2eff_mmq_cr72(r20=r20, pA=pA, dw=aliased_dw, dwH=aliased_dwH, kex=kex, cpmg_frqs=self.cpmg_frqs[ei], inv_tcpmg=self.inv_relax_times[ei], tcp=self.tau_cpmg[ei], back_calc=self.back_calc[ei])
1582
1583
1584 self.back_calc = self.back_calc*self.disp_struct
1585
1586
1587 if self.has_missing:
1588
1589 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1590
1591
1592 return chi2_rankN(self.values, self.back_calc, self.errors)
1593
1594
1596 """Target function for no exchange.
1597
1598 @param params: The vector of parameter values.
1599 @type params: numpy rank-1 float array
1600 @return: The chi-squared value.
1601 @rtype: float
1602 """
1603
1604
1605 if self.scaling_flag:
1606 params = dot(params, self.scaling_matrix)
1607
1608
1609 R20 = params
1610
1611
1612 self.back_calc[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1613
1614
1615 self.back_calc = self.back_calc*self.disp_struct
1616
1617
1618 if self.has_missing:
1619
1620 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1621
1622
1623 return chi2_rankN(self.values, self.back_calc, self.errors)
1624
1625
1627 """Target function for no exchange, for R1rho off resonance models.
1628
1629 @param params: The vector of parameter values.
1630 @type params: numpy rank-1 float array
1631 @return: The chi-squared value.
1632 @rtype: float
1633 """
1634
1635
1636 if self.scaling_flag:
1637 params = dot(params, self.scaling_matrix)
1638
1639
1640 r1rho_prime = params
1641
1642
1643 return self.calc_NOREX_R1RHO(R1=self.r1, r1rho_prime=r1rho_prime)
1644
1645
1647 """Target function for no exchange, for R1rho off resonance models, whereby R1 is fitted.
1648
1649 @param params: The vector of parameter values.
1650 @type params: numpy rank-1 float array
1651 @return: The chi-squared value.
1652 @rtype: float
1653 """
1654
1655
1656 if self.scaling_flag:
1657 params = dot(params, self.scaling_matrix)
1658
1659
1660 r1 = params[:self.end_index[0]]
1661 r1rho_prime = params[self.end_index[0]:self.end_index[1]]
1662
1663
1664 self.r1_struct[:] = multiply.outer( r1.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1665
1666
1667 return self.calc_NOREX_R1RHO(R1=self.r1_struct, r1rho_prime=r1rho_prime)
1668
1669
1671 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations.
1672
1673 @param params: The vector of parameter values.
1674 @type params: numpy rank-1 float array
1675 @return: The chi-squared value.
1676 @rtype: float
1677 """
1678
1679
1680 if self.scaling_flag:
1681 params = dot(params, self.scaling_matrix)
1682
1683
1684 R20 = params[:self.end_index[0]]
1685 dw = params[self.end_index[0]:self.end_index[1]]
1686 pA = params[self.end_index[1]]
1687 kex = params[self.end_index[1]+1]
1688
1689
1690 return self.calc_ns_cpmg_2site_3D_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)
1691
1692
1694 """Target function for the full numerical solution for the 2-site Bloch-McConnell equations.
1695
1696 @param params: The vector of parameter values.
1697 @type params: numpy rank-1 float array
1698 @return: The chi-squared value.
1699 @rtype: float
1700 """
1701
1702
1703 if self.scaling_flag:
1704 params = dot(params, self.scaling_matrix)
1705
1706
1707 R20 = params[:self.end_index[1]].reshape(self.NS*2, self.NM)
1708 R20A = R20[::2].flatten()
1709 R20B = R20[1::2].flatten()
1710 dw = params[self.end_index[1]:self.end_index[2]]
1711 pA = params[self.end_index[2]]
1712 kex = params[self.end_index[2]+1]
1713
1714
1715 return self.calc_ns_cpmg_2site_3D_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)
1716
1717
1719 """Target function for the numerical solution for the 2-site Bloch-McConnell equations using the expanded notation.
1720
1721 @param params: The vector of parameter values.
1722 @type params: numpy rank-1 float array
1723 @return: The chi-squared value.
1724 @rtype: float
1725 """
1726
1727
1728 if self.scaling_flag:
1729 params = dot(params, self.scaling_matrix)
1730
1731
1732 R20 = params[:self.end_index[0]]
1733 dw = params[self.end_index[0]:self.end_index[1]]
1734 pA = params[self.end_index[1]]
1735 kex = params[self.end_index[1]+1]
1736
1737
1738 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1739
1740
1741 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1742
1743
1744 r2eff_ns_cpmg_2site_expanded(r20=self.r20_struct, pA=pA, dw=self.dw_struct, dw_orig=dw, kex=kex, relax_time=self.relax_times, inv_relax_time=self.inv_relax_times, tcp=self.tau_cpmg, back_calc=self.back_calc, num_cpmg=self.power)
1745
1746
1747 self.back_calc = self.back_calc*self.disp_struct
1748
1749
1750 if self.has_missing:
1751
1752 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1753
1754
1755 return chi2_rankN(self.values, self.back_calc, self.errors)
1756
1757
1759 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations using complex conjugate matrices.
1760
1761 This is the model whereby the simplification R20A = R20B is assumed.
1762
1763
1764 @param params: The vector of parameter values.
1765 @type params: numpy rank-1 float array
1766 @return: The chi-squared value.
1767 @rtype: float
1768 """
1769
1770
1771 if self.scaling_flag:
1772 params = dot(params, self.scaling_matrix)
1773
1774
1775 R20 = params[:self.end_index[0]]
1776 dw = params[self.end_index[0]:self.end_index[1]]
1777 pA = params[self.end_index[1]]
1778 kex = params[self.end_index[1]+1]
1779
1780
1781 return self.calc_ns_cpmg_2site_star_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)
1782
1783
1785 """Target function for the full numerical solution for the 2-site Bloch-McConnell equations using complex conjugate matrices.
1786
1787 @param params: The vector of parameter values.
1788 @type params: numpy rank-1 float array
1789 @return: The chi-squared value.
1790 @rtype: float
1791 """
1792
1793
1794 if self.scaling_flag:
1795 params = dot(params, self.scaling_matrix)
1796
1797
1798 R20 = params[:self.end_index[1]].reshape(self.NS*2, self.NM)
1799 R20A = R20[::2].flatten()
1800 R20B = R20[1::2].flatten()
1801 dw = params[self.end_index[1]:self.end_index[2]]
1802 pA = params[self.end_index[2]]
1803 kex = params[self.end_index[2]+1]
1804
1805
1806 return self.calc_ns_cpmg_2site_star_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)
1807
1808
1810 """Target function for the combined SQ, ZQ, DQ and MQ CPMG numeric solution.
1811
1812 @param params: The vector of parameter values.
1813 @type params: numpy rank-1 float array
1814 @return: The chi-squared value.
1815 @rtype: float
1816 """
1817
1818
1819 if self.scaling_flag:
1820 params = dot(params, self.scaling_matrix)
1821
1822
1823 R20 = params[:self.end_index[0]]
1824 dw = params[self.end_index[0]:self.end_index[1]]
1825 dwH = params[self.end_index[1]:self.end_index[2]]
1826 pA = params[self.end_index[2]]
1827 kex = params[self.end_index[2]+1]
1828
1829 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
1830 multiply( multiply.outer( dwH.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs_H, out=self.dwH_struct )
1831
1832
1833 self.r20_struct[:] = multiply.outer( R20.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1834
1835
1836 for ei in range(self.NE):
1837 r20 = self.r20_struct[ei]
1838 dw_frq = self.dw_struct[ei]
1839 dwH_frq = self.dwH_struct[ei]
1840
1841
1842 aliased_dwH = 0.0
1843 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ:
1844 aliased_dw = dw_frq
1845 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ:
1846 aliased_dw = dwH_frq
1847 elif self.exp_types[ei] == EXP_TYPE_CPMG_DQ:
1848 aliased_dw = dw_frq + dwH_frq
1849 elif self.exp_types[ei] == EXP_TYPE_CPMG_ZQ:
1850 aliased_dw = dw_frq - dwH_frq
1851 elif self.exp_types[ei] == EXP_TYPE_CPMG_MQ:
1852 aliased_dw = dw_frq
1853 aliased_dwH = dwH_frq
1854 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_MQ:
1855 aliased_dw = dwH_frq
1856 aliased_dwH = dw_frq
1857
1858
1859 self.r2eff_ns_mmq[ei](M0=self.M0, R20A=r20, R20B=r20, pA=pA, dw=aliased_dw, dwH=aliased_dwH, kex=kex, inv_tcpmg=self.inv_relax_times[ei], tcp=self.tau_cpmg[ei], back_calc=self.back_calc[ei], num_points=self.num_disp_points[ei], power=self.power[ei])
1860
1861
1862 self.back_calc = self.back_calc*self.disp_struct
1863
1864
1865 if self.has_missing:
1866
1867 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
1868
1869
1870 return chi2_rankN(self.values, self.back_calc, self.errors)
1871
1872
1874 """Target function for the combined SQ, ZQ, DQ and MQ 3-site MMQ CPMG numeric solution.
1875
1876 @param params: The vector of parameter values.
1877 @type params: numpy rank-1 float array
1878 @return: The chi-squared value.
1879 @rtype: float
1880 """
1881
1882
1883 if self.scaling_flag:
1884 params = dot(params, self.scaling_matrix)
1885
1886
1887 R20 = params[:self.end_index[0]]
1888 dw_AB = params[self.end_index[0]:self.end_index[1]]
1889 dw_BC = params[self.end_index[1]:self.end_index[2]]
1890 dwH_AB = params[self.end_index[2]:self.end_index[3]]
1891 dwH_BC = params[self.end_index[3]:self.end_index[4]]
1892 pA = params[self.end_index[4]]
1893 kex_AB = params[self.end_index[4]+1]
1894 pB = params[self.end_index[4]+2]
1895 kex_BC = params[self.end_index[4]+3]
1896 kex_AC = params[self.end_index[4]+4]
1897
1898
1899 return self.calc_ns_mmq_3site_chi2(R20A=R20, R20B=R20, R20C=R20, dw_AB=dw_AB, dw_BC=dw_BC, dwH_AB=dwH_AB, dwH_BC=dwH_BC, pA=pA, pB=pB, kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=kex_AC)
1900
1901
1903 """Target function for the combined SQ, ZQ, DQ and MQ 3-site linearised MMQ CPMG numeric solution.
1904
1905 @param params: The vector of parameter values.
1906 @type params: numpy rank-1 float array
1907 @return: The chi-squared value.
1908 @rtype: float
1909 """
1910
1911
1912 if self.scaling_flag:
1913 params = dot(params, self.scaling_matrix)
1914
1915
1916 R20 = params[:self.end_index[0]]
1917 dw_AB = params[self.end_index[0]:self.end_index[1]]
1918 dw_BC = params[self.end_index[1]:self.end_index[2]]
1919 dwH_AB = params[self.end_index[2]:self.end_index[3]]
1920 dwH_BC = params[self.end_index[3]:self.end_index[4]]
1921 pA = params[self.end_index[4]]
1922 kex_AB = params[self.end_index[4]+1]
1923 pB = params[self.end_index[4]+2]
1924 kex_BC = params[self.end_index[4]+3]
1925
1926
1927 return self.calc_ns_mmq_3site_chi2(R20A=R20, R20B=R20, R20C=R20, dw_AB=dw_AB, dw_BC=dw_BC, dwH_AB=dwH_AB, dwH_BC=dwH_BC, pA=pA, pB=pB, kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=0.0)
1928
1929
1931 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations for R1rho data.
1932
1933 @param params: The vector of parameter values.
1934 @type params: numpy rank-1 float array
1935 @return: The chi-squared value.
1936 @rtype: float
1937 """
1938
1939
1940 if self.scaling_flag:
1941 params = dot(params, self.scaling_matrix)
1942
1943
1944 r1rho_prime = params[:self.end_index[0]]
1945 dw = params[self.end_index[0]:self.end_index[1]]
1946 pA = params[self.end_index[1]]
1947 kex = params[self.end_index[1]+1]
1948
1949
1950 return self.calc_ns_r1rho_2site(R1=self.r1, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
1951
1952
1954 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations for R1rho data, whereby R1 is fitted.
1955
1956 @param params: The vector of parameter values.
1957 @type params: numpy rank-1 float array
1958 @return: The chi-squared value.
1959 @rtype: float
1960 """
1961
1962
1963 if self.scaling_flag:
1964 params = dot(params, self.scaling_matrix)
1965
1966
1967 r1 = params[:self.end_index[0]]
1968 r1rho_prime = params[self.end_index[0]:self.end_index[1]]
1969 dw = params[self.end_index[1]:self.end_index[2]]
1970 pA = params[self.end_index[2]]
1971 kex = params[self.end_index[2]+1]
1972
1973
1974 self.r1_struct[:] = multiply.outer( r1.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
1975
1976
1977 return self.calc_ns_r1rho_2site(R1=self.r1_struct, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
1978
1979
1981 """Target function for the R1rho 3-site numeric solution.
1982
1983 @param params: The vector of parameter values.
1984 @type params: numpy rank-1 float array
1985 @return: The chi-squared value.
1986 @rtype: float
1987 """
1988
1989
1990 if self.scaling_flag:
1991 params = dot(params, self.scaling_matrix)
1992
1993
1994 r1rho_prime = params[:self.end_index[0]]
1995 dw_AB = params[self.end_index[0]:self.end_index[1]]
1996 dw_BC = params[self.end_index[1]:self.end_index[2]]
1997 pA = params[self.end_index[2]]
1998 kex_AB = params[self.end_index[2]+1]
1999 pB = params[self.end_index[2]+2]
2000 kex_BC = params[self.end_index[2]+3]
2001 kex_AC = params[self.end_index[2]+4]
2002
2003
2004 return self.calc_ns_r1rho_3site_chi2(r1rho_prime=r1rho_prime, dw_AB=dw_AB, dw_BC=dw_BC, pA=pA, pB=pB, kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=kex_AC)
2005
2006
2008 """Target function for the R1rho 3-site numeric solution linearised with kAC = kCA = 0.
2009
2010 @param params: The vector of parameter values.
2011 @type params: numpy rank-1 float array
2012 @return: The chi-squared value.
2013 @rtype: float
2014 """
2015
2016
2017 if self.scaling_flag:
2018 params = dot(params, self.scaling_matrix)
2019
2020
2021 r1rho_prime = params[:self.end_index[0]]
2022 dw_AB = params[self.end_index[0]:self.end_index[1]]
2023 dw_BC = params[self.end_index[1]:self.end_index[2]]
2024 pA = params[self.end_index[2]]
2025 kex_AB = params[self.end_index[2]+1]
2026 pB = params[self.end_index[2]+2]
2027 kex_BC = params[self.end_index[2]+3]
2028
2029
2030 return self.calc_ns_r1rho_3site_chi2(r1rho_prime=r1rho_prime, dw_AB=dw_AB, dw_BC=dw_BC, pA=pA, pB=pB, kex_AB=kex_AB, kex_BC=kex_BC, kex_AC=0.0)
2031
2032
2034 """Target function for the Trott, Abergel and Palmer (2003) R1rho off-resonance 2-site model.
2035
2036 @param params: The vector of parameter values.
2037 @type params: numpy rank-1 float array
2038 @return: The chi-squared value.
2039 @rtype: float
2040 """
2041
2042
2043 if self.scaling_flag:
2044 params = dot(params, self.scaling_matrix)
2045
2046
2047 r1rho_prime = params[:self.end_index[0]]
2048 dw = params[self.end_index[0]:self.end_index[1]]
2049 pA = params[self.end_index[1]]
2050 kex = params[self.end_index[1]+1]
2051
2052
2053 return self.calc_TAP03(R1=self.r1, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
2054
2055
2057 """Target function for the Trott, Abergel and Palmer (2003) R1rho off-resonance 2-site model, whereby R1 is fitted.
2058
2059 @param params: The vector of parameter values.
2060 @type params: numpy rank-1 float array
2061 @return: The chi-squared value.
2062 @rtype: float
2063 """
2064
2065
2066 if self.scaling_flag:
2067 params = dot(params, self.scaling_matrix)
2068
2069
2070 r1 = params[:self.end_index[0]]
2071 r1rho_prime = params[self.end_index[0]:self.end_index[1]]
2072 dw = params[self.end_index[1]:self.end_index[2]]
2073 pA = params[self.end_index[2]]
2074 kex = params[self.end_index[2]+1]
2075
2076
2077 self.r1_struct[:] = multiply.outer( r1.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
2078
2079
2080 return self.calc_TAP03(R1=self.r1_struct, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
2081
2082
2084 """Target function for the Trott and Palmer (2002) R1rho off-resonance 2-site model.
2085
2086 @param params: The vector of parameter values.
2087 @type params: numpy rank-1 float array
2088 @return: The chi-squared value.
2089 @rtype: float
2090 """
2091
2092
2093 if self.scaling_flag:
2094 params = dot(params, self.scaling_matrix)
2095
2096
2097 r1rho_prime = params[:self.end_index[0]]
2098 dw = params[self.end_index[0]:self.end_index[1]]
2099 pA = params[self.end_index[1]]
2100 kex = params[self.end_index[1]+1]
2101
2102
2103 return self.calc_TP02(R1=self.r1, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
2104
2105
2107 """Target function for the Trott and Palmer (2002) R1rho off-resonance 2-site model, whereby R1 is fitted.
2108
2109 @param params: The vector of parameter values.
2110 @type params: numpy rank-1 float array
2111 @return: The chi-squared value.
2112 @rtype: float
2113 """
2114
2115
2116 if self.scaling_flag:
2117 params = dot(params, self.scaling_matrix)
2118
2119
2120 r1 = params[:self.end_index[0]]
2121 r1rho_prime = params[self.end_index[0]:self.end_index[1]]
2122 dw = params[self.end_index[1]:self.end_index[2]]
2123 pA = params[self.end_index[2]]
2124 kex = params[self.end_index[2]+1]
2125
2126
2127 self.r1_struct[:] = multiply.outer( r1.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
2128
2129
2130 return self.calc_TP02(R1=self.r1_struct, r1rho_prime=r1rho_prime, dw=dw, pA=pA, kex=kex)
2131
2132
2134 """Target function for the the Tollinger et al. (2001) 2-site very-slow exchange model, range of microsecond to second time scale.
2135
2136 @param params: The vector of parameter values.
2137 @type params: numpy rank-1 float array
2138 @return: The chi-squared value.
2139 @rtype: float
2140 """
2141
2142
2143 if self.scaling_flag:
2144 params = dot(params, self.scaling_matrix)
2145
2146
2147 R20A = params[:self.end_index[0]]
2148 dw = params[self.end_index[0]:self.end_index[1]]
2149 k_AB = params[self.end_index[1]]
2150
2151
2152 multiply( multiply.outer( dw.reshape(1, self.NS), self.nm_no_nd_ones ), self.frqs, out=self.dw_struct )
2153
2154
2155 self.r20a_struct[:] = multiply.outer( R20A.reshape(self.NE, self.NS, self.NM), self.no_nd_ones )
2156
2157
2158 r2eff_TSMFK01(r20a=self.r20a_struct, dw=self.dw_struct, dw_orig=dw, k_AB=k_AB, tcp=self.tau_cpmg, back_calc=self.back_calc)
2159
2160
2161 self.back_calc = self.back_calc*self.disp_struct
2162
2163
2164 if self.has_missing:
2165
2166 self.back_calc[self.mask_replace_blank.mask] = self.values[self.mask_replace_blank.mask]
2167
2168
2169 return chi2_rankN(self.values, self.back_calc, self.errors)
2170
2171
2173 """Class function to return back_calc as lists of lists. Number of values in should match number of dispersion points or spin_lock.
2174
2175 @return: back calculation of the R2eff/R1rho values in structure of list of lists. The dimensions are {Ei, Si, Mi, Oi, Di}.
2176 @rtype: rank-4 list of numpy rank-1 float arrays
2177 """
2178
2179 back_calc_return = deepcopy(self.values_orig)
2180
2181
2182 for ei in range(self.NE):
2183 exp_type = self.exp_types[ei]
2184 for si in range(self.NS):
2185 for mi in range(self.NM):
2186 for oi in range(self.NO):
2187 if exp_type in EXP_TYPE_LIST_CPMG:
2188 num = len(self.cpmg_frqs_orig[ei][mi][oi])
2189 else:
2190 num = len(self.spin_lock_nu1_orig[ei][mi][oi])
2191 back_calc_return[ei][si][mi][oi][:] = self.back_calc[ei, si, mi, oi, :num]
2192
2193 return back_calc_return
2194