Package target_functions :: Module relax_disp
[hide private]
[frames] | no frames]

Source Code for Module target_functions.relax_disp

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