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