target_functions.relax

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):

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 always zero. 109 - Di: The index for each dispersion point (either the spin-lock field strength or the nu_CPMG frequency). 110 111 112 @keyword model: The relaxation dispersion model to fit. 113 @type model: str 114 @keyword num_param: The number of parameters in the model. 115 @type num_param: int 116 @keyword num_spins: The number of spins in the cluster. 117 @type num_spins: int 118 @keyword num_frq: The number of spectrometer field strengths. 119 @type num_frq: int 120 @keyword exp_types: The list of experiment types. The dimensions are {Ei}. 121 @type exp_types: list of str 122 @keyword values: The R2eff/R1rho values. The dimensions are {Ei, Si, Mi, Oi, Di}. 123 @type values: rank-4 list of numpy rank-1 float arrays 124 @keyword errors: The R2eff/R1rho errors. The dimensions are {Ei, Si, Mi, Oi, Di}. 125 @type errors: rank-4 list of numpy rank-1 float arrays 126 @keyword missing: The data structure indicating missing R2eff/R1rho data. The dimensions are {Ei, Si, Mi, Oi, Di}. 127 @type missing: rank-4 list of numpy rank-1 int arrays 128 @keyword frqs: The spin Larmor frequencies (in MHz*2pi to speed up the ppm to rad/s conversion). The dimensions are {Ei, Si, Mi}. 129 @type frqs: rank-3 list of floats 130 @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}. 131 @type frqs_H: rank-3 list of floats 132 @keyword cpmg_frqs: The CPMG frequencies in Hertz. This will be ignored for R1rho experiments. The dimensions are {Ei, Mi}. 133 @type cpmg_frqs: rank-2 list of floats 134 @keyword spin_lock_nu1: The spin-lock field strengths in Hertz. This will be ignored for CPMG experiments. The dimensions are {Ei, Mi}. 135 @type spin_lock_nu1: rank-2 list of floats 136 @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}. 137 @type chemical_shifts: rank-3 list of floats 138 @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}. 139 @type offset: rank-4 list of floats 140 @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}. 141 @type tilt_angles: rank-5 list of floats 142 @keyword r1: The R1 relaxation rates. This is only used for off-resonance R1rho models. The dimensions are {Ei, Si, Mi}. 143 @type r1: rank-3 list of floats 144 @keyword relax_times: The experiment specific fixed time period for relaxation (in seconds). The dimensions are {Ei, Mi}. 145 @type relax_times: rank-2 list of floats 146 @keyword scaling_matrix: The square and diagonal scaling matrix. 147 @type scaling_matrix: numpy rank-2 float array 148 @keyword recalc_tau: A flag which if True will cause tau_CPMG to be recalculated to remove user input truncation. 149 @type recalc_tau: bool 150 """ 151 152 # Check the args. 153 if model not in MODEL_LIST_FULL: 154 raise RelaxError("The model '%s' is unknown." % model) 155 if values == None: 156 raise RelaxError("No values have been supplied to the target function.") 157 if errors == None: 158 raise RelaxError("No errors have been supplied to the target function.") 159 if missing == None: 160 raise RelaxError("No missing data information has been supplied to the target function.") 161 if model in [MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05]: 162 if chemical_shifts == None: 163 raise RelaxError("Chemical shifts must be supplied for the '%s' R1rho off-resonance dispersion model." % model) 164 if r1 == None: 165 raise RelaxError("R1 relaxation rates must be supplied for the '%s' R1rho off-resonance dispersion model." % model) 166 167 # Store the arguments. 168 self.model = model 169 self.num_params = num_params 170 self.num_spins = num_spins 171 self.num_frq = num_frq 172 self.exp_types = exp_types 173 self.values = values 174 self.errors = errors 175 self.missing = missing 176 self.frqs = frqs 177 self.frqs_H = frqs_H 178 self.cpmg_frqs = cpmg_frqs 179 self.spin_lock_nu1 = spin_lock_nu1 180 self.chemical_shifts = chemical_shifts 181 self.offset = offset 182 self.tilt_angles = tilt_angles 183 self.r1 = r1 184 self.relax_times = relax_times 185 self.scaling_matrix = scaling_matrix 186 187 # Create the structure for holding the back-calculated R2eff values (matching the dimensions of the values structure). 188 self.back_calc = deepcopy(values) 189 190 # Check the experiment types, simplifying the data structures as needed. 191 self.experiment_type_setup() 192 193 # Determine the number of dispersion points (for each experiment type, magnetic field strength, and offset) and offsets. 194 self.num_disp_points = [] 195 self.num_offsets = [] 196 for ei in range(len(values)): 197 self.num_disp_points.append([]) 198 self.num_offsets.append([]) 199 for si in range(len(values[ei])): 200 self.num_disp_points[ei].append([]) 201 self.num_offsets[ei].append([]) 202 for mi in range(len(values[ei][0])): 203 self.num_disp_points[ei][si].append([]) 204 self.num_offsets[ei][si].append([]) 205 206 # The offset data. 207 if len(offset[ei][si][mi]): 208 self.num_offsets[ei][si][mi] = len(self.offset[ei][si][mi]) 209 else: 210 self.num_offsets[ei][si][mi] = 0 211 212 # The dispersion points. 213 for oi in range(self.num_offsets[ei][si][mi]): 214 self.num_disp_points[ei][si][mi].append([]) 215 if cpmg_frqs != None and len(cpmg_frqs[ei][mi][oi]): 216 self.num_disp_points[ei][si][mi][oi] = len(self.cpmg_frqs[ei][mi][oi]) 217 elif spin_lock_nu1 != None and len(spin_lock_nu1[ei][mi][oi]): 218 self.num_disp_points[ei][si][mi][oi] = len(self.spin_lock_nu1[ei][mi][oi]) 219 else: 220 self.num_disp_points[ei][si][mi][oi] = 0 221 222 # Scaling initialisation. 223 self.scaling_flag = False 224 if self.scaling_matrix != None: 225 self.scaling_flag = True 226 227 # Initialise the post spin parameter indices. 228 self.end_index = [] 229 230 # The spin and frequency dependent R2 parameters. 231 self.end_index.append(self.num_exp * self.num_spins * self.num_frq) 232 if model in [MODEL_B14_FULL, MODEL_CR72_FULL, MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_STAR_FULL]: 233 self.end_index.append(2 * self.num_exp * self.num_spins * self.num_frq) 234 235 # The spin and dependent parameters (phi_ex, dw, padw2). 236 self.end_index.append(self.end_index[-1] + self.num_spins) 237 if model in [MODEL_IT99, MODEL_LM63_3SITE, MODEL_MMQ_CR72, MODEL_NS_MMQ_2SITE]: 238 self.end_index.append(self.end_index[-1] + self.num_spins) 239 elif model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 240 self.end_index.append(self.end_index[-1] + self.num_spins) 241 self.end_index.append(self.end_index[-1] + self.num_spins) 242 elif model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]: 243 self.end_index.append(self.end_index[-1] + self.num_spins) 244 self.end_index.append(self.end_index[-1] + self.num_spins) 245 self.end_index.append(self.end_index[-1] + self.num_spins) 246 247 # Set up the matrices for the numerical solutions. 248 if model in [MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL]: 249 # The matrix that contains only the R2 relaxation terms ("Redfield relaxation", i.e. non-exchange broadening). 250 self.Rr = zeros((2, 2), complex64) 251 252 # The matrix that contains the exchange terms between the two states A and B. 253 self.Rex = zeros((2, 2), complex64) 254 255 # The matrix that contains the chemical shift evolution. It works here only with X magnetization, and the complex notation allows to evolve in the transverse plane (x, y). 256 self.RCS = zeros((2, 2), complex64) 257 258 # The matrix that contains all the contributions to the evolution, i.e. relaxation, exchange and chemical shift evolution. 259 self.R = zeros((2, 2), complex64) 260 261 # Pi-pulse propagators. 262 if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL]: 263 self.r180x = r180x_3d() 264 265 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 266 if model in [MODEL_NS_CPMG_2SITE_STAR, MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_MMQ_2SITE]: 267 self.M0 = zeros(2, float64) 268 if model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]: 269 self.M0 = zeros(3, float64) 270 if model in [MODEL_NS_CPMG_2SITE_3D, MODEL_NS_CPMG_2SITE_3D_FULL]: 271 self.M0 = zeros(7, float64) 272 self.M0[0] = 0.5 273 if model in [MODEL_NS_R1RHO_2SITE]: 274 self.M0 = zeros(6, float64) 275 if model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 276 self.M0 = zeros(9, float64) 277 278 # Special CPMG-type data structures. 279 if model in [MODEL_B14, MODEL_B14_FULL, MODEL_MMQ_CR72, 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_TSMFK01]: 280 # The number of CPMG blocks. 281 self.power = [] 282 for ei in range(self.num_exp): 283 self.power.append([]) 284 for mi in range(self.num_frq): 285 self.power[ei].append(zeros(self.num_disp_points[ei][si][mi][0], int16)) 286 for di in range(self.num_disp_points[ei][si][mi][0]): 287 # Missing data for an entire field strength. 288 if isNaN(self.relax_times[ei][mi]): 289 self.power[ei][mi][di] = 0.0 290 291 # Normal value. 292 else: 293 self.power[ei][mi][di] = int(round(self.cpmg_frqs[ei][mi][0][di] * self.relax_times[ei][mi])) 294 295 # The tau_cpmg times. 296 self.tau_cpmg = [] 297 for ei in range(len(values)): 298 self.tau_cpmg.append([]) 299 for mi in range(len(values[ei][0])): 300 self.tau_cpmg[ei].append(zeros(self.num_disp_points[ei][si][mi][0], float64)) 301 for di in range(self.num_disp_points[ei][si][mi][0]): 302 # Recalculate the tau_cpmg times to avoid any user induced truncation in the input files. 303 if recalc_tau: 304 self.tau_cpmg[ei][mi][di] = 0.25 * self.relax_times[ei][mi] / self.power[ei][mi][di] 305 else: 306 self.tau_cpmg[ei][mi][di] = 0.25 / self.cpmg_frqs[ei][mi][0][di] 307 308 # Convert the spin-lock data to rad.s^-1. 309 if spin_lock_nu1 != None: 310 self.spin_lock_omega1 = [] 311 self.spin_lock_omega1_squared = [] 312 for ei in range(len(values)): 313 self.spin_lock_omega1.append([]) 314 self.spin_lock_omega1_squared.append([]) 315 for mi in range(len(values[ei][0])): 316 self.spin_lock_omega1[ei].append([]) 317 self.spin_lock_omega1_squared[ei].append([]) 318 for oi in range(len(values[ei][0][mi])): 319 self.spin_lock_omega1[ei][mi].append([]) 320 self.spin_lock_omega1_squared[ei][mi].append([]) 321 self.spin_lock_omega1[ei][mi][oi] = 2.0 * pi * self.spin_lock_nu1[ei][mi][oi] 322 self.spin_lock_omega1_squared[ei][mi][oi] = self.spin_lock_omega1[ei][mi][oi] ** 2 323 324 # The inverted relaxation delay. 325 if model in [MODEL_B14, MODEL_B14_FULL, MODEL_MMQ_CR72, 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]: 326 self.inv_relax_times = 1.0 / relax_times 327 328 # Special storage matrices for the multi-quantum CPMG N-site numerical models. 329 if model == MODEL_NS_MMQ_2SITE: 330 self.m1 = zeros((2, 2), complex64) 331 self.m2 = zeros((2, 2), complex64) 332 elif model in [MODEL_NS_MMQ_3SITE, MODEL_NS_MMQ_3SITE_LINEAR]: 333 self.m1 = zeros((3, 3), complex64) 334 self.m2 = zeros((3, 3), complex64) 335 elif model == MODEL_NS_R1RHO_2SITE: 336 self.matrix = zeros((6, 6), float64) 337 elif model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: 338 self.matrix = zeros((9, 9), float64) 339 340 # Set up the model. 341 if model == MODEL_NOREX: 342 self.func = self.func_NOREX 343 if model == MODEL_LM63: 344 self.func = self.func_LM63 345 if model == MODEL_LM63_3SITE: 346 self.func = self.func_LM63_3site 347 if model == MODEL_CR72_FULL: 348 self.func = self.func_CR72_full 349 if model == MODEL_CR72: 350 self.func = self.func_CR72 351 if model == MODEL_IT99: 352 self.func = self.func_IT99 353 if model == MODEL_TSMFK01: 354 self.func = self.func_TSMFK01 355 if model == MODEL_B14: 356 self.func = self.func_B14 357 if model == MODEL_B14_FULL: 358 self.func = self.func_B14_full 359 if model == MODEL_NS_CPMG_2SITE_3D_FULL: 360 self.func = self.func_ns_cpmg_2site_3D_full 361 if model == MODEL_NS_CPMG_2SITE_3D: 362 self.func = self.func_ns_cpmg_2site_3D 363 if model == MODEL_NS_CPMG_2SITE_EXPANDED: 364 self.func = self.func_ns_cpmg_2site_expanded 365 if model == MODEL_NS_CPMG_2SITE_STAR_FULL: 366 self.func = self.func_ns_cpmg_2site_star_full 367 if model == MODEL_NS_CPMG_2SITE_STAR: 368 self.func = self.func_ns_cpmg_2site_star 369 if model == MODEL_M61: 370 self.func = self.func_M61 371 if model == MODEL_M61B: 372 self.func = self.func_M61b 373 if model == MODEL_DPL94: 374 self.func = self.func_DPL94 375 if model == MODEL_TP02: 376 self.func = self.func_TP02 377 if model == MODEL_TAP03: 378 self.func = self.func_TAP03 379 if model == MODEL_MP05: 380 self.func = self.func_MP05 381 if model == MODEL_NS_R1RHO_2SITE: 382 self.func = self.func_ns_r1rho_2site 383 if model == MODEL_NS_R1RHO_3SITE: 384 self.func = self.func_ns_r1rho_3site 385 if model == MODEL_NS_R1RHO_3SITE_LINEAR: 386 self.func = self.func_ns_r1rho_3site_linear 387 if model == MODEL_MMQ_CR72: 388 self.func = self.func_mmq_CR72 389 if model == MODEL_NS_MMQ_2SITE: 390 self.func = self.func_ns_mmq_2site 391 if model == MODEL_NS_MMQ_3SITE: 392 self.func = self.func_ns_mmq_3site 393 if model == MODEL_NS_MMQ_3SITE_LINEAR: 394 self.func = self.func_ns_mmq_3site_linear

395 396

397 - def calc_B14_chi2(self, R20A=None, R20B=None, dw=None, pA=None, kex=None):

398 """Calculate the chi-squared value of the Baldwin (2014) 2-site exact solution model for all time scales. 399 400 401 @keyword R20A: The R2 value for state A in the absence of exchange. 402 @type R20A: list of float 403 @keyword R20B: The R2 value for state B in the absence of exchange. 404 @type R20B: list of float 405 @keyword dw: The chemical shift differences in ppm for each spin. 406 @type dw: list of float 407 @keyword pA: The population of state A. 408 @type pA: float 409 @keyword kex: The rate of exchange. 410 @type kex: float 411 @return: The chi-squared value. 412 @rtype: float 413 """ 414 415 # Once off parameter conversions. 416 pB = 1.0 - pA 417 k_BA = pA * kex 418 k_AB = pB * kex 419 420 # Initialise. 421 chi2_sum = 0.0 422 423 # Loop over the experiment types. 424 for ei in range(self.num_exp): 425 # Loop over the spins. 426 for si in range(self.num_spins): 427 # Loop over the spectrometer frequencies. 428 for mi in range(self.num_frq): 429 # The R20 index. 430 r20_index = mi + si*self.num_frq 431 432 # Convert dw from ppm to rad/s. 433 dw_frq = dw[si] * self.frqs[ei][si][mi] 434 435 # Alias the dw frequency combinations. 436 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ: 437 aliased_dw = dw_frq 438 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ: 439 aliased_dw = dw_frq 440 441 # Back calculate the R2eff values. 442 r2eff_B14(r20a=R20A[r20_index], r20b=R20B[r20_index], pA=pA, pB=pB, dw=aliased_dw, kex=kex, k_AB=k_AB, k_BA=k_BA, ncyc=self.power[ei][mi], inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi], back_calc=self.back_calc[ei][si][mi][0], num_points=self.num_disp_points[ei][si][mi][0]) 443 444 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 445 for di in range(self.num_disp_points[ei][si][mi][0]): 446 if self.missing[ei][si][mi][0][di]: 447 self.back_calc[ei][si][mi][0][di] = self.values[ei][si][mi][0][di] 448 449 # Calculate and return the chi-squared value. 450 chi2_sum += chi2(self.values[ei][si][mi][0], self.back_calc[ei][si][mi][0], self.errors[ei][si][mi][0]) 451 452 # Return the total chi-squared value. 453 return chi2_sum

454 455

456 - def calc_CR72_chi2(self, R20A=None, R20B=None, dw=None, pA=None, kex=None):

457 """Calculate the chi-squared value of the Carver and Richards (1972) 2-site exchange model on all time scales. 458 459 @keyword R20A: The R2 value for state A in the absence of exchange. 460 @type R20A: list of float 461 @keyword R20B: The R2 value for state B in the absence of exchange. 462 @type R20B: list of float 463 @keyword dw: The chemical shift differences in ppm for each spin. 464 @type dw: list of float 465 @keyword pA: The population of state A. 466 @type pA: float 467 @keyword kex: The rate of exchange. 468 @type kex: float 469 @return: The chi-squared value. 470 @rtype: float 471 """ 472 473 # Initialise. 474 chi2_sum = 0.0 475 476 # Loop over the spins. 477 for si in range(self.num_spins): 478 # Loop over the spectrometer frequencies. 479 for mi in range(self.num_frq): 480 # The R20 index. 481 r20_index = mi + si*self.num_frq 482 483 # Convert dw from ppm to rad/s. 484 dw_frq = dw[si] * self.frqs[0][si][mi] 485 486 # Back calculate the R2eff values. 487 r2eff_CR72(r20a=R20A[r20_index], r20b=R20B[r20_index], pA=pA, dw=dw_frq, kex=kex, cpmg_frqs=self.cpmg_frqs[0][mi][0], back_calc = self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 488 489 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 490 for di in range(self.num_disp_points[0][si][mi][0]): 491 if self.missing[0][si][mi][0][di]: 492 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 493 494 # Calculate and return the chi-squared value. 495 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 496 497 # Return the total chi-squared value. 498 return chi2_sum

499 500

501 - def calc_ns_cpmg_2site_3D_chi2(self, R20A=None, R20B=None, dw=None, pA=None, kex=None):

502 """Calculate the chi-squared value of the 'NS CPMG 2-site' models. 503 504 @keyword R20A: The R2 value for state A in the absence of exchange. 505 @type R20A: list of float 506 @keyword R20B: The R2 value for state B in the absence of exchange. 507 @type R20B: list of float 508 @keyword dw: The chemical shift differences in ppm for each spin. 509 @type dw: list of float 510 @keyword pA: The population of state A. 511 @type pA: float 512 @keyword kex: The rate of exchange. 513 @type kex: float 514 @return: The chi-squared value. 515 @rtype: float 516 """ 517 518 # Once off parameter conversions. 519 pB = 1.0 - pA 520 k_BA = pA * kex 521 k_AB = pB * kex 522 523 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 524 self.M0[1] = pA 525 self.M0[4] = pB 526 527 # Chi-squared initialisation. 528 chi2_sum = 0.0 529 530 # Loop over the spins. 531 for si in range(self.num_spins): 532 # Loop over the spectrometer frequencies. 533 for mi in range(self.num_frq): 534 # The R20 index. 535 r20_index = mi + si*self.num_frq 536 537 # Convert dw from ppm to rad/s. 538 dw_frq = dw[si] * self.frqs[0][si][mi] 539 540 # Back calculate the R2eff values. 541 r2eff_ns_cpmg_2site_3D(r180x=self.r180x, M0=self.M0, r20a=R20A[r20_index], r20b=R20B[r20_index], pA=pA, pB=pB, dw=dw_frq, k_AB=k_AB, k_BA=k_BA, inv_tcpmg=self.inv_relax_times[0][mi], tcp=self.tau_cpmg[0][mi], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0], power=self.power[0][mi]) 542 543 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 544 for di in range(self.num_disp_points[0][si][mi][0]): 545 if self.missing[0][si][mi][0][di]: 546 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 547 548 # Calculate and return the chi-squared value. 549 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 550 551 # Return the total chi-squared value. 552 return chi2_sum

553 554

555 - def calc_ns_cpmg_2site_star_chi2(self, R20A=None, R20B=None, dw=None, pA=None, kex=None):

556 """Calculate the chi-squared value of the 'NS CPMG 2-site star' models. 557 558 @keyword R20A: The R2 value for state A in the absence of exchange. 559 @type R20A: list of float 560 @keyword R20B: The R2 value for state B in the absence of exchange. 561 @type R20B: list of float 562 @keyword dw: The chemical shift differences in ppm for each spin. 563 @type dw: list of float 564 @keyword pA: The population of state A. 565 @type pA: float 566 @keyword kex: The rate of exchange. 567 @type kex: float 568 @return: The chi-squared value. 569 @rtype: float 570 """ 571 572 # Once off parameter conversions. 573 pB = 1.0 - pA 574 k_BA = pA * kex 575 k_AB = pB * kex 576 577 # Set up the matrix that contains the exchange terms between the two states A and B. 578 self.Rex[0, 0] = -k_AB 579 self.Rex[0, 1] = k_BA 580 self.Rex[1, 0] = k_AB 581 self.Rex[1, 1] = -k_BA 582 583 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 584 self.M0[0] = pA 585 self.M0[1] = pB 586 587 # Chi-squared initialisation. 588 chi2_sum = 0.0 589 590 # Loop over the spins. 591 for si in range(self.num_spins): 592 # Loop over the spectrometer frequencies. 593 for mi in range(self.num_frq): 594 # The R20 index. 595 r20_index = mi + si*self.num_frq 596 597 # Convert dw from ppm to rad/s. 598 dw_frq = dw[si] * self.frqs[0][si][mi] 599 600 # Back calculate the R2eff values. 601 r2eff_ns_cpmg_2site_star(Rr=self.Rr, Rex=self.Rex, RCS=self.RCS, R=self.R, M0=self.M0, r20a=R20A[r20_index], r20b=R20B[r20_index], dw=dw_frq, inv_tcpmg=self.inv_relax_times[0][mi], tcp=self.tau_cpmg[0][mi], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0], power=self.power[0][mi]) 602 603 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 604 for di in range(self.num_disp_points[0][si][mi][0]): 605 if self.missing[0][si][mi][0][di]: 606 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 607 608 # Calculate and return the chi-squared value. 609 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 610 611 # Return the total chi-squared value. 612 return chi2_sum

613 614

615 - 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):

616 """Calculate the chi-squared value for the 'NS MMQ 3-site' models. 617 618 @keyword R20A: The R2 value for state A in the absence of exchange. 619 @type R20A: list of float 620 @keyword R20B: The R2 value for state B in the absence of exchange. 621 @type R20B: list of float 622 @keyword R20C: The R2 value for state C in the absence of exchange. 623 @type R20C: list of float 624 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s. 625 @type dw_AB: float 626 @keyword dw_BC: The chemical exchange difference between states B and C in rad/s. 627 @type dw_BC: float 628 @keyword dwH_AB: The proton chemical exchange difference between states A and B in rad/s. 629 @type dwH_AB: float 630 @keyword dwH_BC: The proton chemical exchange difference between states B and C in rad/s. 631 @type dwH_BC: float 632 @keyword pA: The population of state A. 633 @type pA: float 634 @keyword kex_AB: The rate of exchange between states A and B. 635 @type kex_AB: float 636 @keyword kex_BC: The rate of exchange between states B and C. 637 @type kex_BC: float 638 @keyword kex_AC: The rate of exchange between states A and C. 639 @type kex_AC: float 640 @return: The chi-squared value. 641 @rtype: float 642 """ 643 644 # Once off parameter conversions. 645 pC = 1.0 - pA - pB 646 pA_pB = pA + pB 647 pA_pC = pA + pC 648 pB_pC = pB + pC 649 k_BA = pA * kex_AB / pA_pB 650 k_AB = pB * kex_AB / pA_pB 651 k_CB = pB * kex_BC / pB_pC 652 k_BC = pC * kex_BC / pB_pC 653 k_CA = pA * kex_AC / pA_pC 654 k_AC = pC * kex_AC / pA_pC 655 dw_AC = dw_AB + dw_BC 656 dwH_AC = dwH_AB + dwH_BC 657 658 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 659 self.M0[0] = pA 660 self.M0[1] = pB 661 self.M0[2] = pC 662 663 # Initialise. 664 chi2_sum = 0.0 665 666 # Loop over the experiment types. 667 for ei in range(self.num_exp): 668 # Loop over the spins. 669 for si in range(self.num_spins): 670 # Loop over the spectrometer frequencies. 671 for mi in range(self.num_frq): 672 # The R20 index. 673 r20_index = mi + ei*self.num_frq + si*self.num_frq*self.num_exp 674 675 # Convert dw from ppm to rad/s. 676 dw_AB_frq = dw_AB[si] * self.frqs[ei][si][mi] 677 dw_AC_frq = dw_AC[si] * self.frqs[ei][si][mi] 678 dwH_AB_frq = dwH_AB[si] * self.frqs_H[ei][si][mi] 679 dwH_AC_frq = dwH_AC[si] * self.frqs_H[ei][si][mi] 680 681 # Alias the dw frequency combinations. 682 aliased_dwH_AB = 0.0 683 aliased_dwH_AC = 0.0 684 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ: 685 aliased_dw_AB = dw_AB_frq 686 aliased_dw_AC = dw_AC_frq 687 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ: 688 aliased_dw_AB = dwH_AB_frq 689 aliased_dw_AC = dwH_AC_frq 690 elif self.exp_types[ei] == EXP_TYPE_CPMG_DQ: 691 aliased_dw_AB = dw_AB_frq + dwH_AB_frq 692 aliased_dw_AC = dw_AC_frq + dwH_AC_frq 693 elif self.exp_types[ei] == EXP_TYPE_CPMG_ZQ: 694 aliased_dw_AB = dw_AB_frq - dwH_AB_frq 695 aliased_dw_AC = dw_AC_frq - dwH_AC_frq 696 elif self.exp_types[ei] == EXP_TYPE_CPMG_MQ: 697 aliased_dw_AB = dw_AB_frq 698 aliased_dw_AC = dw_AC_frq 699 aliased_dwH_AB = dwH_AB_frq 700 aliased_dwH_AC = dwH_AC_frq 701 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_MQ: 702 aliased_dw_AB = dwH_AB_frq 703 aliased_dw_AC = dwH_AC_frq 704 aliased_dwH_AB = dw_AB_frq 705 aliased_dwH_AC = dw_AC_frq 706 707 # Back calculate the R2eff values for each experiment type. 708 self.r2eff_ns_mmq[ei](M0=self.M0, m1=self.m1, m2=self.m2, R20A=R20A[r20_index], R20B=R20B[r20_index], R20C=R20C[r20_index], pA=pA, pB=pB, pC=pC, dw_AB=aliased_dw_AB, dw_AC=aliased_dw_AC, dwH_AB=aliased_dwH_AB, dwH_AC=aliased_dwH_AC, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi], back_calc=self.back_calc[ei][si][mi][0], num_points=self.num_disp_points[ei][si][mi][0], power=self.power[ei][mi]) 709 710 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 711 for di in range(self.num_disp_points[ei][si][mi][0]): 712 if self.missing[ei][si][mi][0][di]: 713 self.back_calc[ei][si][mi][0][di] = self.values[ei][si][mi][0][di] 714 715 # Calculate and return the chi-squared value. 716 chi2_sum += chi2(self.values[ei][si][mi][0], self.back_calc[ei][si][mi][0], self.errors[ei][si][mi][0]) 717 718 # Return the total chi-squared value. 719 return chi2_sum

720 721

722 - 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):

723 """Calculate the chi-squared value for the 'NS MMQ 3-site' models. 724 725 @keyword r1rho_prime: The R1rho value for all states in the absence of exchange. 726 @type r1rho_prime: list of float 727 @keyword dw_AB: The chemical exchange difference between states A and B in rad/s. 728 @type dw_AB: float 729 @keyword dw_BC: The chemical exchange difference between states B and C in rad/s. 730 @type dw_BC: float 731 @keyword pA: The population of state A. 732 @type pA: float 733 @keyword kex_AB: The rate of exchange between states A and B. 734 @type kex_AB: float 735 @keyword kex_BC: The rate of exchange between states B and C. 736 @type kex_BC: float 737 @keyword kex_AC: The rate of exchange between states A and C. 738 @type kex_AC: float 739 @return: The chi-squared value. 740 @rtype: float 741 """ 742 743 # Once off parameter conversions. 744 pC = 1.0 - pA - pB 745 pA_pB = pA + pB 746 pA_pC = pA + pC 747 pB_pC = pB + pC 748 k_BA = pA * kex_AB / pA_pB 749 k_AB = pB * kex_AB / pA_pB 750 k_CB = pB * kex_BC / pB_pC 751 k_BC = pC * kex_BC / pB_pC 752 k_CA = pA * kex_AC / pA_pC 753 k_AC = pC * kex_AC / pA_pC 754 dw_AC = dw_AB + dw_BC 755 756 # Initialise. 757 chi2_sum = 0.0 758 759 # Loop over the spins. 760 for si in range(self.num_spins): 761 # Loop over the spectrometer frequencies. 762 for mi in range(self.num_frq): 763 # The R20 index. 764 r20_index = mi + si*self.num_frq 765 766 # Convert dw from ppm to rad/s. 767 dw_AB_frq = dw_AB[si] * self.frqs[0][si][mi] 768 dw_AC_frq = dw_AC[si] * self.frqs[0][si][mi] 769 770 # Loop over the offsets. 771 for oi in range(self.num_offsets[0][si][mi]): 772 # Back calculate the R2eff values for each experiment type. 773 ns_r1rho_3site(M0=self.M0, matrix=self.matrix, r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], r1=self.r1[si, mi], pA=pA, pB=pB, pC=pC, dw_AB=dw_AB_frq, dw_AC=dw_AC_frq, k_AB=k_AB, k_BA=k_BA, k_BC=k_BC, k_CB=k_CB, k_AC=k_AC, k_CA=k_CA, spin_lock_fields=self.spin_lock_omega1[0][mi][oi], relax_time=self.relax_times[0][mi], inv_relax_time=self.inv_relax_times[0][mi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) 774 775 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 776 for di in range(self.num_disp_points[0][si][mi][oi]): 777 if self.missing[0][si][mi][oi][di]: 778 self.back_calc[0][si][mi][oi][di] = self.values[0][si][mi][oi][di] 779 780 # Calculate and return the chi-squared value. 781 chi2_sum += chi2(self.values[0][si][mi][oi], self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi]) 782 783 # Return the total chi-squared value. 784 return chi2_sum

785 786

787 - def experiment_type_setup(self):

788 """Check the experiment types and simplify data structures. 789 790 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. 791 """ 792 793 # The number of experiments. 794 self.num_exp = len(self.exp_types) 795 796 # The MMQ combined data type models. 797 if self.model in MODEL_LIST_MMQ: 798 # Alias the r2eff functions. 799 self.r2eff_ns_mmq = [] 800 801 # Loop over the experiment types. 802 for ei in range(self.num_exp): 803 # SQ, DQ and ZQ data types. 804 if self.exp_types[ei] in [EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_DQ, EXP_TYPE_CPMG_ZQ]: 805 if self.model == MODEL_NS_MMQ_2SITE: 806 self.r2eff_ns_mmq.append(r2eff_ns_mmq_2site_sq_dq_zq) 807 else: 808 self.r2eff_ns_mmq.append(r2eff_ns_mmq_3site_sq_dq_zq) 809 810 # MQ data types. 811 elif self.exp_types[ei] in [EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ]: 812 if self.model == MODEL_NS_MMQ_2SITE: 813 self.r2eff_ns_mmq.append(r2eff_ns_mmq_2site_mq) 814 else: 815 self.r2eff_ns_mmq.append(r2eff_ns_mmq_3site_mq) 816 817 # The single data type models. 818 else: 819 # Check that the data is correct. 820 if self.model != MODEL_NOREX and self.model in MODEL_LIST_CPMG and self.exp_types[0] != EXP_TYPE_CPMG_SQ: 821 raise RelaxError("The '%s' CPMG model is not compatible with the '%s' experiment type." % (self.model, self.exp_types[0])) 822 if self.model != MODEL_NOREX and self.model in MODEL_LIST_R1RHO and self.exp_types[0] != EXP_TYPE_R1RHO: 823 raise RelaxError("The '%s' R1rho model is not compatible with the '%s' experiment type." % (self.model, self.exp_types[0])) 824 if self.model != MODEL_NOREX and self.model in MODEL_LIST_MQ_CPMG and self.exp_types[0] != EXP_TYPE_CPMG_MQ: 825 raise RelaxError("The '%s' CPMG model is not compatible with the '%s' experiment type." % (self.model, self.exp_types[0]))

826 827

828 - def func_B14(self, params):

829 """Target function for the Baldwin (2014) 2-site exact solution model for all time scales, whereby the simplification R20A = R20B is assumed. 830 831 This assumes that pA > pB, and hence this must be implemented as a constraint. 832 833 834 @param params: The vector of parameter values. 835 @type params: numpy rank-1 float array 836 @return: The chi-squared value. 837 @rtype: float 838 """ 839 840 # Scaling. 841 if self.scaling_flag: 842 params = dot(params, self.scaling_matrix) 843 844 # Unpack the parameter values. 845 R20 = params[:self.end_index[0]] 846 dw = params[self.end_index[0]:self.end_index[1]] 847 pA = params[self.end_index[1]] 848 kex = params[self.end_index[1]+1] 849 850 # Calculate and return the chi-squared value. 851 return self.calc_B14_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)

852 853

854 - def func_B14_full(self, params):

855 """Target function for the Baldwin (2014) 2-site exact solution model for all time scales. 856 857 This assumes that pA > pB, and hence this must be implemented as a constraint. 858 859 860 @param params: The vector of parameter values. 861 @type params: numpy rank-1 float array 862 @return: The chi-squared value. 863 @rtype: float 864 """ 865 866 # Scaling. 867 if self.scaling_flag: 868 params = dot(params, self.scaling_matrix) 869 870 # Unpack the parameter values. 871 R20 = params[:self.end_index[1]].reshape(self.num_spins*2, self.num_frq) 872 R20A = R20[::2].flatten() 873 R20B = R20[1::2].flatten() 874 dw = params[self.end_index[1]:self.end_index[2]] 875 pA = params[self.end_index[2]] 876 kex = params[self.end_index[2]+1] 877 878 # Calculate and return the chi-squared value. 879 return self.calc_B14_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)

880 881

882 - def func_CR72(self, params):

883 """Target function for the reduced Carver and Richards (1972) 2-site exchange model on all time scales. 884 885 This assumes that pA > pB, and hence this must be implemented as a constraint. For this model, the simplification R20A = R20B is assumed. 886 887 888 @param params: The vector of parameter values. 889 @type params: numpy rank-1 float array 890 @return: The chi-squared value. 891 @rtype: float 892 """ 893 894 # Scaling. 895 if self.scaling_flag: 896 params = dot(params, self.scaling_matrix) 897 898 # Unpack the parameter values. 899 R20 = params[:self.end_index[0]] 900 dw = params[self.end_index[0]:self.end_index[1]] 901 pA = params[self.end_index[1]] 902 kex = params[self.end_index[1]+1] 903 904 # Calculate and return the chi-squared value. 905 return self.calc_CR72_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)

906 907

908 - def func_CR72_full(self, params):

909 """Target function for the full Carver and Richards (1972) 2-site exchange model on all time scales. 910 911 This assumes that pA > pB, and hence this must be implemented as a constraint. 912 913 914 @param params: The vector of parameter values. 915 @type params: numpy rank-1 float array 916 @return: The chi-squared value. 917 @rtype: float 918 """ 919 920 # Scaling. 921 if self.scaling_flag: 922 params = dot(params, self.scaling_matrix) 923 924 # Unpack the parameter values. 925 R20 = params[:self.end_index[1]].reshape(self.num_spins*2, self.num_frq) 926 R20A = R20[::2].flatten() 927 R20B = R20[1::2].flatten() 928 dw = params[self.end_index[1]:self.end_index[2]] 929 pA = params[self.end_index[2]] 930 kex = params[self.end_index[2]+1] 931 932 # Calculate and return the chi-squared value. 933 return self.calc_CR72_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)

934 935

936 - def func_DPL94(self, params):

937 """Target function for the Davis, Perlman and London (1994) fast 2-site off-resonance exchange model for R1rho-type experiments. 938 939 @param params: The vector of parameter values. 940 @type params: numpy rank-1 float array 941 @return: The chi-squared value. 942 @rtype: float 943 """ 944 945 # Scaling. 946 if self.scaling_flag: 947 params = dot(params, self.scaling_matrix) 948 949 # Unpack the parameter values. 950 R20 = params[:self.end_index[0]] 951 phi_ex = params[self.end_index[0]:self.end_index[1]] 952 kex = params[self.end_index[1]] 953 954 # Initialise. 955 chi2_sum = 0.0 956 957 # Loop over the spins. 958 for si in range(self.num_spins): 959 # Loop over the spectrometer frequencies. 960 for mi in range(self.num_frq): 961 # The R20 index. 962 r20_index = mi + si*self.num_frq 963 964 # Convert phi_ex from ppm^2 to (rad/s)^2. 965 phi_ex_scaled = phi_ex[si] * self.frqs[0][si][mi]**2 966 967 # Loop over the offsets. 968 for oi in range(self.num_offsets[0][si][mi]): 969 # Back calculate the R2eff values. 970 r1rho_DPL94(r1rho_prime=R20[r20_index], phi_ex=phi_ex_scaled, kex=kex, theta=self.tilt_angles[0][si][mi][oi], R1=self.r1[si, mi], spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) 971 972 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 973 for di in range(self.num_disp_points[0][si][mi][oi]): 974 if self.missing[0][si][mi][oi][di]: 975 self.back_calc[0][si][mi][oi][di] = self.values[0][si][mi][oi][di] 976 977 # Calculate and return the chi-squared value. 978 chi2_sum += chi2(self.values[0][si][mi][oi], self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi]) 979 980 # Return the total chi-squared value. 981 return chi2_sum

982 983

984 - def func_IT99(self, params):

985 """Target function for the Ishima and Torchia (1999) 2-site model for all timescales with pA >> pB. 986 987 @param params: The vector of parameter values. 988 @type params: numpy rank-1 float array 989 @return: The chi-squared value. 990 @rtype: float 991 """ 992 993 # Scaling. 994 if self.scaling_flag: 995 params = dot(params, self.scaling_matrix) 996 997 # Unpack the parameter values. 998 R20 = params[:self.end_index[0]] 999 dw = params[self.end_index[0]:self.end_index[1]] 1000 pA = params[self.end_index[1]] 1001 tex = params[self.end_index[1]+1] 1002 1003 # Once off parameter conversions. 1004 pB = 1.0 - pA 1005 1006 # Initialise. 1007 chi2_sum = 0.0 1008 1009 # Loop over the spins. 1010 for si in range(self.num_spins): 1011 # Loop over the spectrometer frequencies. 1012 for mi in range(self.num_frq): 1013 # The R20 index. 1014 r20_index = mi + si*self.num_frq 1015 1016 # Convert dw from ppm to rad/s. 1017 dw_frq = dw[si] * self.frqs[0][si][mi] 1018 1019 # Back calculate the R2eff values. 1020 r2eff_IT99(r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq, tex=tex, cpmg_frqs=self.cpmg_frqs[0][mi][0], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 1021 1022 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1023 for di in range(self.num_disp_points[0][si][mi][0]): 1024 if self.missing[0][si][mi][0][di]: 1025 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1026 1027 # Calculate and return the chi-squared value. 1028 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1029 1030 # Return the total chi-squared value. 1031 return chi2_sum

1032 1033

1034 - def func_LM63_3site(self, params):

1035 """Target function for the Luz and Meiboom (1963) fast 3-site exchange model. 1036 1037 @param params: The vector of parameter values. 1038 @type params: numpy rank-1 float array 1039 @return: The chi-squared value. 1040 @rtype: float 1041 """ 1042 1043 # Scaling. 1044 if self.scaling_flag: 1045 params = dot(params, self.scaling_matrix) 1046 1047 # Unpack the parameter values. 1048 R20 = params[:self.end_index[0]] 1049 phi_ex_B = params[self.end_index[0]:self.end_index[1]] 1050 phi_ex_C = params[self.end_index[1]:self.end_index[2]] 1051 kB = params[self.end_index[2]] 1052 kC = params[self.end_index[2]+1] 1053 1054 # Once off parameter conversions. 1055 rex_B = phi_ex_B / kB 1056 rex_C = phi_ex_C / kC 1057 quart_kB = kB / 4.0 1058 quart_kC = kC / 4.0 1059 1060 # Initialise. 1061 chi2_sum = 0.0 1062 1063 # Loop over the spins. 1064 for si in range(self.num_spins): 1065 # Loop over the spectrometer frequencies. 1066 for mi in range(self.num_frq): 1067 # The R20 index. 1068 r20_index = mi + si*self.num_frq 1069 1070 # Convert phi_ex (or rex) from ppm^2 to (rad/s)^2. 1071 frq2 = self.frqs[0][si][mi]**2 1072 rex_B_scaled = rex_B[si] * frq2 1073 rex_C_scaled = rex_C[si] * frq2 1074 1075 # Back calculate the R2eff values. 1076 r2eff_LM63_3site(r20=R20[r20_index], rex_B=rex_B_scaled, rex_C=rex_C_scaled, quart_kB=quart_kB, quart_kC=quart_kC, cpmg_frqs=self.cpmg_frqs[0][mi][0], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 1077 1078 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1079 for di in range(self.num_disp_points[0][si][mi][0]): 1080 if self.missing[0][si][mi][0][di]: 1081 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1082 1083 # Calculate and return the chi-squared value. 1084 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1085 1086 # Return the total chi-squared value. 1087 return chi2_sum

1088 1089

1090 - def func_LM63(self, params):

1091 """Target function for the Luz and Meiboom (1963) fast 2-site exchange model. 1092 1093 @param params: The vector of parameter values. 1094 @type params: numpy rank-1 float array 1095 @return: The chi-squared value. 1096 @rtype: float 1097 """ 1098 1099 # Scaling. 1100 if self.scaling_flag: 1101 params = dot(params, self.scaling_matrix) 1102 1103 # Unpack the parameter values. 1104 R20 = params[:self.end_index[0]] 1105 phi_ex = params[self.end_index[0]:self.end_index[1]] 1106 kex = params[self.end_index[1]] 1107 1108 # Initialise. 1109 chi2_sum = 0.0 1110 1111 # Loop over the spins. 1112 for si in range(self.num_spins): 1113 # Loop over the spectrometer frequencies. 1114 for mi in range(self.num_frq): 1115 # The R20 index. 1116 r20_index = mi + si*self.num_frq 1117 1118 # Convert phi_ex from ppm^2 to (rad/s)^2. 1119 phi_ex_scaled = phi_ex[si] * self.frqs[0][si][mi]**2 1120 1121 # Back calculate the R2eff values. 1122 r2eff_LM63(r20=R20[r20_index], phi_ex=phi_ex_scaled, kex=kex, cpmg_frqs=self.cpmg_frqs[0][mi][0], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 1123 1124 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1125 for di in range(self.num_disp_points[0][si][mi][0]): 1126 if self.missing[0][si][mi][0][di]: 1127 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1128 1129 # Calculate and return the chi-squared value. 1130 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1131 1132 # Return the total chi-squared value. 1133 return chi2_sum

1134 1135

1136 - def func_M61(self, params):

1137 """Target function for the Meiboom (1961) fast 2-site exchange model for R1rho-type experiments. 1138 1139 @param params: The vector of parameter values. 1140 @type params: numpy rank-1 float array 1141 @return: The chi-squared value. 1142 @rtype: float 1143 """ 1144 1145 # Scaling. 1146 if self.scaling_flag: 1147 params = dot(params, self.scaling_matrix) 1148 1149 # Unpack the parameter values. 1150 R20 = params[:self.end_index[0]] 1151 phi_ex = params[self.end_index[0]:self.end_index[1]] 1152 kex = params[self.end_index[1]] 1153 1154 # Initialise. 1155 chi2_sum = 0.0 1156 1157 # Loop over the spins. 1158 for si in range(self.num_spins): 1159 # Loop over the spectrometer frequencies. 1160 for mi in range(self.num_frq): 1161 # The R20 index. 1162 r20_index = mi + si*self.num_frq 1163 1164 # Convert phi_ex from ppm^2 to (rad/s)^2. 1165 phi_ex_scaled = phi_ex[si] * self.frqs[0][si][mi]**2 1166 1167 # Back calculate the R2eff values. 1168 r1rho_M61(r1rho_prime=R20[r20_index], phi_ex=phi_ex_scaled, kex=kex, spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][0], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 1169 1170 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1171 for di in range(self.num_disp_points[0][si][mi][0]): 1172 if self.missing[0][si][mi][0][di]: 1173 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1174 1175 # Calculate and return the chi-squared value. 1176 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1177 1178 # Return the total chi-squared value. 1179 return chi2_sum

1180 1181

1182 - def func_M61b(self, params):

1183 """Target function for the Meiboom (1961) R1rho on-resonance 2-site model for skewed populations (pA >> pB). 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 # Scaling. 1192 if self.scaling_flag: 1193 params = dot(params, self.scaling_matrix) 1194 1195 # Unpack the parameter values. 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 # Initialise. 1202 chi2_sum = 0.0 1203 1204 # Loop over the spins. 1205 for si in range(self.num_spins): 1206 # Loop over the spectrometer frequencies. 1207 for mi in range(self.num_frq): 1208 # The R20 index. 1209 r20_index = mi + si*self.num_frq 1210 1211 # Convert dw from ppm to rad/s. 1212 dw_frq = dw[si] * self.frqs[0][si][mi] 1213 1214 # Back calculate the R1rho values. 1215 r1rho_M61b(r1rho_prime=R20[r20_index], pA=pA, dw=dw_frq, kex=kex, spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][0], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 1216 1217 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1218 for di in range(self.num_disp_points[0][si][mi][0]): 1219 if self.missing[0][si][mi][0][di]: 1220 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1221 1222 # Calculate and return the chi-squared value. 1223 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1224 1225 # Return the total chi-squared value. 1226 return chi2_sum

1227 1228

1229 - def func_MP05(self, params):

1230 """Target function for the Miloushev and Palmer (2005) R1rho off-resonance 2-site model. 1231 1232 @param params: The vector of parameter values. 1233 @type params: numpy rank-1 float array 1234 @return: The chi-squared value. 1235 @rtype: float 1236 """ 1237 1238 # Scaling. 1239 if self.scaling_flag: 1240 params = dot(params, self.scaling_matrix) 1241 1242 # Unpack the parameter values. 1243 R20 = params[:self.end_index[0]] 1244 dw = params[self.end_index[0]:self.end_index[1]] 1245 pA = params[self.end_index[1]] 1246 kex = params[self.end_index[1]+1] 1247 1248 # Once off parameter conversions. 1249 pB = 1.0 - pA 1250 1251 # Initialise. 1252 chi2_sum = 0.0 1253 1254 # Loop over the spins. 1255 for si in range(self.num_spins): 1256 # Loop over the spectrometer frequencies. 1257 for mi in range(self.num_frq): 1258 # The R20 index. 1259 r20_index = mi + si*self.num_frq 1260 1261 # Convert dw from ppm to rad/s. 1262 dw_frq = dw[si] * self.frqs[0][si][mi] 1263 1264 # Loop over the offsets. 1265 for oi in range(self.num_offsets[0][si][mi]): 1266 # Back calculate the R1rho values. 1267 r1rho_MP05(r1rho_prime=R20[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], spin_lock_fields=self.spin_lock_omega1[0][mi][oi], spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) 1268 1269 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1270 for di in range(self.num_disp_points[0][si][mi][oi]): 1271 if self.missing[0][si][mi][oi][di]: 1272 self.back_calc[0][si][mi][oi][di] = self.values[0][si][mi][oi][di] 1273 1274 # Calculate and return the chi-squared value. 1275 chi2_sum += chi2(self.values[0][si][mi][oi], self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi]) 1276 1277 # Return the total chi-squared value. 1278 return chi2_sum

1279 1280

1281 - def func_mmq_CR72(self, params):

1282 """Target function for the CR72 model extended for MQ CPMG data. 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 # Scaling. 1291 if self.scaling_flag: 1292 params = dot(params, self.scaling_matrix) 1293 1294 # Unpack the parameter values. 1295 R20 = params[:self.end_index[0]] 1296 dw = params[self.end_index[0]:self.end_index[1]] 1297 dwH = params[self.end_index[1]:self.end_index[2]] 1298 pA = params[self.end_index[2]] 1299 kex = params[self.end_index[2]+1] 1300 1301 # Once off parameter conversions. 1302 pB = 1.0 - pA 1303 k_BA = pA * kex 1304 k_AB = pB * kex 1305 1306 # Initialise. 1307 chi2_sum = 0.0 1308 1309 # Loop over the experiment types. 1310 for ei in range(self.num_exp): 1311 # Loop over the spins. 1312 for si in range(self.num_spins): 1313 # Loop over the spectrometer frequencies. 1314 for mi in range(self.num_frq): 1315 # The R20 index. 1316 r20_index = mi + ei*self.num_frq + si*self.num_frq*self.num_exp 1317 1318 # Convert dw from ppm to rad/s. 1319 dw_frq = dw[si] * self.frqs[ei][si][mi] 1320 dwH_frq = dwH[si] * self.frqs_H[ei][si][mi] 1321 1322 # Alias the dw frequency combinations. 1323 aliased_dwH = 0.0 1324 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ: 1325 aliased_dw = dw_frq 1326 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ: 1327 aliased_dw = dwH_frq 1328 elif self.exp_types[ei] == EXP_TYPE_CPMG_DQ: 1329 aliased_dw = dw_frq + dwH_frq 1330 elif self.exp_types[ei] == EXP_TYPE_CPMG_ZQ: 1331 aliased_dw = dw_frq - dwH_frq 1332 elif self.exp_types[ei] == EXP_TYPE_CPMG_MQ: 1333 aliased_dw = dw_frq 1334 aliased_dwH = dwH_frq 1335 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_MQ: 1336 aliased_dw = dwH_frq 1337 aliased_dwH = dw_frq 1338 1339 # Back calculate the R2eff values. 1340 r2eff_mmq_cr72(r20=R20[r20_index], pA=pA, pB=pB, dw=aliased_dw, dwH=aliased_dwH, kex=kex, k_AB=k_AB, k_BA=k_BA, cpmg_frqs=self.cpmg_frqs[ei][mi][0], inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi], back_calc=self.back_calc[ei][si][mi][0], num_points=self.num_disp_points[ei][si][mi][0], power=self.power[ei][mi]) 1341 1342 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1343 for di in range(self.num_disp_points[ei][si][mi][0]): 1344 if self.missing[ei][si][mi][0][di]: 1345 self.back_calc[ei][si][mi][0][di] = self.values[ei][si][mi][0][di] 1346 1347 # Calculate and return the chi-squared value. 1348 chi2_sum += chi2(self.values[ei][si][mi][0], self.back_calc[ei][si][mi][0], self.errors[ei][si][mi][0]) 1349 1350 # Return the total chi-squared value. 1351 return chi2_sum

1352 1353

1354 - def func_NOREX(self, params):

1355 """Target function for no exchange. 1356 1357 @param params: The vector of parameter values. 1358 @type params: numpy rank-1 float array 1359 @return: The chi-squared value. 1360 @rtype: float 1361 """ 1362 1363 # Scaling. 1364 if self.scaling_flag: 1365 params = dot(params, self.scaling_matrix) 1366 1367 # Unpack the parameter values. 1368 R20 = params 1369 1370 # Initialise. 1371 chi2_sum = 0.0 1372 1373 # Loop over the experiment types. 1374 for ei in range(self.num_exp): 1375 # Loop over the spins. 1376 for si in range(self.num_spins): 1377 # Loop over the spectrometer frequencies. 1378 for mi in range(self.num_frq): 1379 # The R20 index. 1380 r20_index = mi + si*self.num_frq 1381 1382 # Loop over the offsets. 1383 for oi in range(self.num_offsets[ei][si][mi]): 1384 # The R2eff values as R20 values. 1385 for di in range(self.num_disp_points[ei][si][mi][oi]): 1386 self.back_calc[ei][si][mi][oi][di] = R20[r20_index] 1387 1388 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1389 for di in range(self.num_disp_points[ei][si][mi][oi]): 1390 if self.missing[ei][si][mi][oi][di]: 1391 self.back_calc[ei][si][mi][oi][di] = self.values[ei][si][mi][oi][di] 1392 1393 # Calculate and return the chi-squared value. 1394 chi2_sum += chi2(self.values[ei][si][mi][oi], self.back_calc[ei][si][mi][oi], self.errors[ei][si][mi][oi]) 1395 1396 # Return the total chi-squared value. 1397 return chi2_sum

1398 1399

1400 - def func_ns_cpmg_2site_3D(self, params):

1401 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations. 1402 1403 @param params: The vector of parameter values. 1404 @type params: numpy rank-1 float array 1405 @return: The chi-squared value. 1406 @rtype: float 1407 """ 1408 1409 # Scaling. 1410 if self.scaling_flag: 1411 params = dot(params, self.scaling_matrix) 1412 1413 # Unpack the parameter values. 1414 R20 = params[:self.end_index[0]] 1415 dw = params[self.end_index[0]:self.end_index[1]] 1416 pA = params[self.end_index[1]] 1417 kex = params[self.end_index[1]+1] 1418 1419 # Calculate and return the chi-squared value. 1420 return self.calc_ns_cpmg_2site_3D_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)

1421 1422

1423 - def func_ns_cpmg_2site_3D_full(self, params):

1424 """Target function for the full numerical solution for the 2-site Bloch-McConnell equations. 1425 1426 @param params: The vector of parameter values. 1427 @type params: numpy rank-1 float array 1428 @return: The chi-squared value. 1429 @rtype: float 1430 """ 1431 1432 # Scaling. 1433 if self.scaling_flag: 1434 params = dot(params, self.scaling_matrix) 1435 1436 # Unpack the parameter values. 1437 R20 = params[:self.end_index[1]].reshape(self.num_spins*2, self.num_frq) 1438 R20A = R20[::2].flatten() 1439 R20B = R20[1::2].flatten() 1440 dw = params[self.end_index[1]:self.end_index[2]] 1441 pA = params[self.end_index[2]] 1442 kex = params[self.end_index[2]+1] 1443 1444 # Calculate and return the chi-squared value. 1445 return self.calc_ns_cpmg_2site_3D_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)

1446 1447

1448 - def func_ns_cpmg_2site_expanded(self, params):

1449 """Target function for the numerical solution for the 2-site Bloch-McConnell equations using the expanded notation. 1450 1451 @param params: The vector of parameter values. 1452 @type params: numpy rank-1 float array 1453 @return: The chi-squared value. 1454 @rtype: float 1455 """ 1456 1457 # Scaling. 1458 if self.scaling_flag: 1459 params = dot(params, self.scaling_matrix) 1460 1461 # Unpack the parameter values. 1462 R20 = params[:self.end_index[0]] 1463 dw = params[self.end_index[0]:self.end_index[1]] 1464 pA = params[self.end_index[1]] 1465 kex = params[self.end_index[1]+1] 1466 1467 # Once off parameter conversions. 1468 pB = 1.0 - pA 1469 k_BA = pA * kex 1470 k_AB = pB * kex 1471 1472 # Chi-squared initialisation. 1473 chi2_sum = 0.0 1474 1475 # Loop over the spins. 1476 for si in range(self.num_spins): 1477 # Loop over the spectrometer frequencies. 1478 for mi in range(self.num_frq): 1479 # The R20 index. 1480 r20_index = mi + si*self.num_frq 1481 1482 # Convert dw from ppm to rad/s. 1483 dw_frq = dw[si] * self.frqs[0][si][mi] 1484 1485 # Back calculate the R2eff values. 1486 r2eff_ns_cpmg_2site_expanded(r20=R20[r20_index], pA=pA, dw=dw_frq, k_AB=k_AB, k_BA=k_BA, relax_time=self.relax_times[0][mi], inv_relax_time=self.inv_relax_times[0][mi], tcp=self.tau_cpmg[0][mi], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0], num_cpmg=self.power[0][mi]) 1487 1488 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1489 for di in range(self.num_disp_points[0][si][mi][0]): 1490 if self.missing[0][si][mi][0][di]: 1491 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1492 1493 # Calculate and return the chi-squared value. 1494 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1495 1496 # Return the total chi-squared value. 1497 return chi2_sum

1498 1499

1500 - def func_ns_cpmg_2site_star(self, params):

1501 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations using complex conjugate matrices. 1502 1503 This is the model whereby the simplification R20A = R20B is assumed. 1504 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 # Scaling. 1513 if self.scaling_flag: 1514 params = dot(params, self.scaling_matrix) 1515 1516 # Unpack the parameter values. 1517 R20 = params[:self.end_index[0]] 1518 dw = params[self.end_index[0]:self.end_index[1]] 1519 pA = params[self.end_index[1]] 1520 kex = params[self.end_index[1]+1] 1521 1522 # Calculate and return the chi-squared value. 1523 return self.calc_ns_cpmg_2site_star_chi2(R20A=R20, R20B=R20, dw=dw, pA=pA, kex=kex)

1524 1525

1526 - def func_ns_cpmg_2site_star_full(self, params):

1527 """Target function for the full numerical solution for the 2-site Bloch-McConnell equations using complex conjugate matrices. 1528 1529 @param params: The vector of parameter values. 1530 @type params: numpy rank-1 float array 1531 @return: The chi-squared value. 1532 @rtype: float 1533 """ 1534 1535 # Scaling. 1536 if self.scaling_flag: 1537 params = dot(params, self.scaling_matrix) 1538 1539 # Unpack the parameter values. 1540 R20 = params[:self.end_index[1]].reshape(self.num_spins*2, self.num_frq) 1541 R20A = R20[::2].flatten() 1542 R20B = R20[1::2].flatten() 1543 dw = params[self.end_index[1]:self.end_index[2]] 1544 pA = params[self.end_index[2]] 1545 kex = params[self.end_index[2]+1] 1546 1547 # Calculate and return the chi-squared value. 1548 return self.calc_ns_cpmg_2site_star_chi2(R20A=R20A, R20B=R20B, dw=dw, pA=pA, kex=kex)

1549 1550

1551 - def func_ns_mmq_2site(self, params):

1552 """Target function for the combined SQ, ZQ, DQ and MQ CPMG numeric solution. 1553 1554 @param params: The vector of parameter values. 1555 @type params: numpy rank-1 float array 1556 @return: The chi-squared value. 1557 @rtype: float 1558 """ 1559 1560 # Scaling. 1561 if self.scaling_flag: 1562 params = dot(params, self.scaling_matrix) 1563 1564 # Unpack the parameter values. 1565 R20 = params[:self.end_index[0]] 1566 dw = params[self.end_index[0]:self.end_index[1]] 1567 dwH = params[self.end_index[1]:self.end_index[2]] 1568 pA = params[self.end_index[2]] 1569 kex = params[self.end_index[2]+1] 1570 1571 # Once off parameter conversions. 1572 pB = 1.0 - pA 1573 k_BA = pA * kex 1574 k_AB = pB * kex 1575 1576 # This is a vector that contains the initial magnetizations corresponding to the A and B state transverse magnetizations. 1577 self.M0[0] = pA 1578 self.M0[1] = pB 1579 1580 # Initialise. 1581 chi2_sum = 0.0 1582 1583 # Loop over the experiment types. 1584 for ei in range(self.num_exp): 1585 # Loop over the spins. 1586 for si in range(self.num_spins): 1587 # Loop over the spectrometer frequencies. 1588 for mi in range(self.num_frq): 1589 # The R20 index. 1590 r20_index = mi + ei*self.num_frq + si*self.num_frq*self.num_exp 1591 1592 # Convert dw from ppm to rad/s. 1593 dw_frq = dw[si] * self.frqs[ei][si][mi] 1594 dwH_frq = dwH[si] * self.frqs_H[ei][si][mi] 1595 1596 # Alias the dw frequency combinations. 1597 aliased_dwH = 0.0 1598 if self.exp_types[ei] == EXP_TYPE_CPMG_SQ: 1599 aliased_dw = dw_frq 1600 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_SQ: 1601 aliased_dw = dwH_frq 1602 elif self.exp_types[ei] == EXP_TYPE_CPMG_DQ: 1603 aliased_dw = dw_frq + dwH_frq 1604 elif self.exp_types[ei] == EXP_TYPE_CPMG_ZQ: 1605 aliased_dw = dw_frq - dwH_frq 1606 elif self.exp_types[ei] == EXP_TYPE_CPMG_MQ: 1607 aliased_dw = dw_frq 1608 aliased_dwH = dwH_frq 1609 elif self.exp_types[ei] == EXP_TYPE_CPMG_PROTON_MQ: 1610 aliased_dw = dwH_frq 1611 aliased_dwH = dw_frq 1612 1613 # Back calculate the R2eff values for each experiment type. 1614 self.r2eff_ns_mmq[ei](M0=self.M0, m1=self.m1, m2=self.m2, R20A=R20[r20_index], R20B=R20[r20_index], pA=pA, pB=pB, dw=aliased_dw, dwH=aliased_dwH, k_AB=k_AB, k_BA=k_BA, inv_tcpmg=self.inv_relax_times[ei][mi], tcp=self.tau_cpmg[ei][mi], back_calc=self.back_calc[ei][si][mi][0], num_points=self.num_disp_points[ei][si][mi][0], power=self.power[ei][mi]) 1615 1616 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1617 for di in range(self.num_disp_points[ei][si][mi][0]): 1618 if self.missing[ei][si][mi][0][di]: 1619 self.back_calc[ei][si][mi][0][di] = self.values[ei][si][mi][0][di] 1620 1621 # Calculate and return the chi-squared value. 1622 chi2_sum += chi2(self.values[ei][si][mi][0], self.back_calc[ei][si][mi][0], self.errors[ei][si][mi][0]) 1623 1624 # Return the total chi-squared value. 1625 return chi2_sum

1626 1627

1628 - def func_ns_mmq_3site(self, params):

1629 """Target function for the combined SQ, ZQ, DQ and MQ 3-site MMQ CPMG numeric solution. 1630 1631 @param params: The vector of parameter values. 1632 @type params: numpy rank-1 float array 1633 @return: The chi-squared value. 1634 @rtype: float 1635 """ 1636 1637 # Scaling. 1638 if self.scaling_flag: 1639 params = dot(params, self.scaling_matrix) 1640 1641 # Unpack the parameter values. 1642 R20 = params[:self.end_index[0]] 1643 dw_AB = params[self.end_index[0]:self.end_index[1]] 1644 dw_BC = params[self.end_index[1]:self.end_index[2]] 1645 dwH_AB = params[self.end_index[2]:self.end_index[3]] 1646 dwH_BC = params[self.end_index[3]:self.end_index[4]] 1647 pA = params[self.end_index[4]] 1648 kex_AB = params[self.end_index[4]+1] 1649 pB = params[self.end_index[4]+2] 1650 kex_BC = params[self.end_index[4]+3] 1651 kex_AC = params[self.end_index[4]+4] 1652 1653 # Calculate and return the chi-squared value. 1654 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)

1655 1656

1657 - def func_ns_mmq_3site_linear(self, params):

1658 """Target function for the combined SQ, ZQ, DQ and MQ 3-site linearised MMQ CPMG numeric solution. 1659 1660 @param params: The vector of parameter values. 1661 @type params: numpy rank-1 float array 1662 @return: The chi-squared value. 1663 @rtype: float 1664 """ 1665 1666 # Scaling. 1667 if self.scaling_flag: 1668 params = dot(params, self.scaling_matrix) 1669 1670 # Unpack the parameter values. 1671 R20 = params[:self.end_index[0]] 1672 dw_AB = params[self.end_index[0]:self.end_index[1]] 1673 dw_BC = params[self.end_index[1]:self.end_index[2]] 1674 dwH_AB = params[self.end_index[2]:self.end_index[3]] 1675 dwH_BC = params[self.end_index[3]:self.end_index[4]] 1676 pA = params[self.end_index[4]] 1677 kex_AB = params[self.end_index[4]+1] 1678 pB = params[self.end_index[4]+2] 1679 kex_BC = params[self.end_index[4]+3] 1680 1681 # Calculate and return the chi-squared value. 1682 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)

1683 1684

1685 - def func_ns_r1rho_2site(self, params):

1686 """Target function for the reduced numerical solution for the 2-site Bloch-McConnell equations for R1rho data. 1687 1688 @param params: The vector of parameter values. 1689 @type params: numpy rank-1 float array 1690 @return: The chi-squared value. 1691 @rtype: float 1692 """ 1693 1694 # Scaling. 1695 if self.scaling_flag: 1696 params = dot(params, self.scaling_matrix) 1697 1698 # Unpack the parameter values. 1699 r1rho_prime = params[:self.end_index[0]] 1700 dw = params[self.end_index[0]:self.end_index[1]] 1701 pA = params[self.end_index[1]] 1702 kex = params[self.end_index[1]+1] 1703 1704 # Once off parameter conversions. 1705 pB = 1.0 - pA 1706 k_BA = pA * kex 1707 k_AB = pB * kex 1708 1709 # Chi-squared initialisation. 1710 chi2_sum = 0.0 1711 1712 # Loop over the spins. 1713 for si in range(self.num_spins): 1714 # Loop over the spectrometer frequencies. 1715 for mi in range(self.num_frq): 1716 # The R20 index. 1717 r20_index = mi + si*self.num_frq 1718 1719 # Convert dw from ppm to rad/s. 1720 dw_frq = dw[si] * self.frqs[0][si][mi] 1721 1722 # Loop over the offsets. 1723 for oi in range(self.num_offsets[0][si][mi]): 1724 # Back calculate the R2eff values. 1725 ns_r1rho_2site(M0=self.M0, matrix=self.matrix, r1rho_prime=r1rho_prime[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], r1=self.r1[si, mi], pA=pA, pB=pB, dw=dw_frq, k_AB=k_AB, k_BA=k_BA, spin_lock_fields=self.spin_lock_omega1[0][mi][oi], relax_time=self.relax_times[0][mi], inv_relax_time=self.inv_relax_times[0][mi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) 1726 1727 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1728 for di in range(self.num_disp_points[0][si][mi][oi]): 1729 if self.missing[0][si][mi][oi][di]: 1730 self.back_calc[0][si][mi][oi][di] = self.values[0][si][mi][oi][di] 1731 1732 # Calculate and return the chi-squared value. 1733 chi2_sum += chi2(self.values[0][si][mi][oi], self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi]) 1734 1735 # Return the total chi-squared value. 1736 return chi2_sum

1737 1738

1739 - def func_ns_r1rho_3site(self, params):

1740 """Target function for the R1rho 3-site numeric solution. 1741 1742 @param params: The vector of parameter values. 1743 @type params: numpy rank-1 float array 1744 @return: The chi-squared value. 1745 @rtype: float 1746 """ 1747 1748 # Scaling. 1749 if self.scaling_flag: 1750 params = dot(params, self.scaling_matrix) 1751 1752 # Unpack the parameter values. 1753 r1rho_prime = params[:self.end_index[0]] 1754 dw_AB = params[self.end_index[0]:self.end_index[1]] 1755 dw_BC = params[self.end_index[1]:self.end_index[2]] 1756 pA = params[self.end_index[2]] 1757 kex_AB = params[self.end_index[2]+1] 1758 pB = params[self.end_index[2]+2] 1759 kex_BC = params[self.end_index[2]+3] 1760 kex_AC = params[self.end_index[2]+4] 1761 1762 # Calculate and return the chi-squared value. 1763 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)

1764 1765

1766 - def func_ns_r1rho_3site_linear(self, params):

1767 """Target function for the R1rho 3-site numeric solution linearised with kAC = kCA = 0. 1768 1769 @param params: The vector of parameter values. 1770 @type params: numpy rank-1 float array 1771 @return: The chi-squared value. 1772 @rtype: float 1773 """ 1774 1775 # Scaling. 1776 if self.scaling_flag: 1777 params = dot(params, self.scaling_matrix) 1778 1779 # Unpack the parameter values. 1780 r1rho_prime = params[:self.end_index[0]] 1781 dw_AB = params[self.end_index[0]:self.end_index[1]] 1782 dw_BC = params[self.end_index[1]:self.end_index[2]] 1783 pA = params[self.end_index[2]] 1784 kex_AB = params[self.end_index[2]+1] 1785 pB = params[self.end_index[2]+2] 1786 kex_BC = params[self.end_index[2]+3] 1787 1788 # Calculate and return the chi-squared value. 1789 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)

1790 1791

1792 - def func_TAP03(self, params):

1793 """Target function for the Trott, Abergel and Palmer (2003) R1rho off-resonance 2-site model. 1794 1795 @param params: The vector of parameter values. 1796 @type params: numpy rank-1 float array 1797 @return: The chi-squared value. 1798 @rtype: float 1799 """ 1800 1801 # Scaling. 1802 if self.scaling_flag: 1803 params = dot(params, self.scaling_matrix) 1804 1805 # Unpack the parameter values. 1806 R20 = params[:self.end_index[0]] 1807 dw = params[self.end_index[0]:self.end_index[1]] 1808 pA = params[self.end_index[1]] 1809 kex = params[self.end_index[1]+1] 1810 1811 # Once off parameter conversions. 1812 pB = 1.0 - pA 1813 1814 # Initialise. 1815 chi2_sum = 0.0 1816 1817 # Loop over the spins. 1818 for si in range(self.num_spins): 1819 # Loop over the spectrometer frequencies. 1820 for mi in range(self.num_frq): 1821 # The R20 index. 1822 r20_index = mi + si*self.num_frq 1823 1824 # Convert dw from ppm to rad/s. 1825 dw_frq = dw[si] * self.frqs[0][si][mi] 1826 1827 # Loop over the offsets. 1828 for oi in range(self.num_offsets[0][si][mi]): 1829 # Back calculate the R1rho values. 1830 r1rho_TAP03(r1rho_prime=R20[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], spin_lock_fields=self.spin_lock_omega1[0][mi][oi], spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) 1831 1832 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1833 for di in range(self.num_disp_points[0][si][mi][oi]): 1834 if self.missing[0][si][mi][oi][di]: 1835 self.back_calc[0][si][mi][oi][di] = self.values[0][si][mi][oi][di] 1836 1837 # Calculate and return the chi-squared value. 1838 chi2_sum += chi2(self.values[0][si][mi][oi], self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi]) 1839 1840 # Return the total chi-squared value. 1841 return chi2_sum

1842 1843

1844 - def func_TP02(self, params):

1845 """Target function for the Trott and Palmer (2002) R1rho off-resonance 2-site model. 1846 1847 @param params: The vector of parameter values. 1848 @type params: numpy rank-1 float array 1849 @return: The chi-squared value. 1850 @rtype: float 1851 """ 1852 1853 # Scaling. 1854 if self.scaling_flag: 1855 params = dot(params, self.scaling_matrix) 1856 1857 # Unpack the parameter values. 1858 R20 = params[:self.end_index[0]] 1859 dw = params[self.end_index[0]:self.end_index[1]] 1860 pA = params[self.end_index[1]] 1861 kex = params[self.end_index[1]+1] 1862 1863 # Once off parameter conversions. 1864 pB = 1.0 - pA 1865 1866 # Initialise. 1867 chi2_sum = 0.0 1868 1869 # Loop over the spins. 1870 for si in range(self.num_spins): 1871 # Loop over the spectrometer frequencies. 1872 for mi in range(self.num_frq): 1873 # The R20 index. 1874 r20_index = mi + si*self.num_frq 1875 1876 # Convert dw from ppm to rad/s. 1877 dw_frq = dw[si] * self.frqs[0][si][mi] 1878 1879 # Loop over the offsets. 1880 for oi in range(self.num_offsets[0][si][mi]): 1881 # Back calculate the R1rho values. 1882 r1rho_TP02(r1rho_prime=R20[r20_index], omega=self.chemical_shifts[0][si][mi], offset=self.offset[0][si][mi][oi], pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi], spin_lock_fields=self.spin_lock_omega1[0][mi][oi], spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi], back_calc=self.back_calc[0][si][mi][oi], num_points=self.num_disp_points[0][si][mi][oi]) 1883 1884 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1885 for di in range(self.num_disp_points[0][si][mi][oi]): 1886 if self.missing[0][si][mi][oi][di]: 1887 self.back_calc[0][si][mi][oi][di] = self.values[0][si][mi][oi][di] 1888 1889 # Calculate and return the chi-squared value. 1890 chi2_sum += chi2(self.values[0][si][mi][oi], self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi]) 1891 1892 # Return the total chi-squared value. 1893 return chi2_sum

1894 1895

1896 - def func_TSMFK01(self, params):

1897 """Target function for the the Tollinger et al. (2001) 2-site very-slow exchange model, range of microsecond to second time scale. 1898 1899 @param params: The vector of parameter values. 1900 @type params: numpy rank-1 float array 1901 @return: The chi-squared value. 1902 @rtype: float 1903 """ 1904 1905 # Scaling. 1906 if self.scaling_flag: 1907 params = dot(params, self.scaling_matrix) 1908 1909 # Unpack the parameter values. 1910 R20A = params[:self.end_index[0]] 1911 dw = params[self.end_index[0]:self.end_index[1]] 1912 k_AB = params[self.end_index[1]] 1913 1914 # Initialise. 1915 chi2_sum = 0.0 1916 1917 # Loop over the spins. 1918 for si in range(self.num_spins): 1919 # Loop over the spectrometer frequencies. 1920 for mi in range(self.num_frq): 1921 # The R20 index. 1922 r20a_index = mi + si*self.num_frq 1923 1924 # Convert dw from ppm to rad/s. 1925 dw_frq = dw[si] * self.frqs[0][si][mi] 1926 1927 # Back calculate the R2eff values. 1928 r2eff_TSMFK01(r20a=R20A[r20a_index], dw=dw_frq, k_AB=k_AB, tcp=self.tau_cpmg[0][mi], back_calc=self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) 1929 1930 # For all missing data points, set the back-calculated value to the measured values so that it has no effect on the chi-squared value. 1931 for di in range(self.num_disp_points[0][si][mi][0]): 1932 if self.missing[0][si][mi][0][di]: 1933 self.back_calc[0][si][mi][0][di] = self.values[0][si][mi][0][di] 1934 1935 # Calculate and return the chi-squared value. 1936 chi2_sum += chi2(self.values[0][si][mi][0], self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) 1937 1938 # Return the total chi-squared value. 1939 return chi2_sum

Source Code for Module target_functions.relax_disp