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