specific_analyses.relax

63 """Class containing functions for relaxation dispersion curve fitting.""" 64

65 - def __init__(self):

66 """Initialise the class by placing API_common methods into the API.""" 67 68 # Execute the base class __init__ method. 69 super(Relax_disp, self).__init__() 70 71 # Place methods into the API. 72 self.data_init = self._data_init_spin 73 self.model_type = self._model_type_local 74 self.return_conversion_factor = self._return_no_conversion_factor 75 self.return_value = self.return_value 76 self.set_param_values = self._set_param_values_spin 77 78 # Set up the spin parameters. 79 self.PARAMS.add('intensities', scope='spin', py_type=dict, grace_string='\\qPeak intensities\\Q') 80 self.PARAMS.add('relax_times', scope='spin', py_type=dict, grace_string='\\qRelaxation time period (s)\\Q') 81 self.PARAMS.add('cpmg_frqs', scope='spin', py_type=dict, grace_string='\\qCPMG pulse train frequency (Hz)\\Q') 82 self.PARAMS.add('spin_lock_nu1', scope='spin', py_type=dict, grace_string='\\qSpin-lock field strength (Hz)\\Q') 83 self.PARAMS.add('r2eff', scope='spin', default=15.0, desc='The effective transversal relaxation rate', set='params', py_type=dict, grace_string='\\qR\\s2,eff\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 84 self.PARAMS.add('i0', scope='spin', default=10000.0, desc='The initial intensity', py_type=dict, set='params', grace_string='\\qI\\s0\\Q', err=True, sim=True) 85 self.PARAMS.add('r2', scope='spin', default=15.0, desc='The transversal relaxation rate', set='params', py_type=dict, grace_string='\\qR\\s2\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 86 self.PARAMS.add('r2a', scope='spin', default=15.0, desc='The transversal relaxation rate for state A in the absence of exchange', set='params', py_type=dict, grace_string='\\qR\\s2,A\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 87 self.PARAMS.add('r2b', scope='spin', default=15.0, desc='The transversal relaxation rate for state B in the absence of exchange', set='params', py_type=dict, grace_string='\\qR\\s2,B\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 88 self.PARAMS.add('pA', scope='spin', default=0.5, desc='The population for state A', set='params', py_type=float, grace_string='\\qp\\sA\\N\\Q', err=True, sim=True) 89 self.PARAMS.add('pB', scope='spin', default=0.5, desc='The population for state B', set='params', py_type=float, grace_string='\\qp\\sB\\N\\Q', err=True, sim=True) 90 self.PARAMS.add('pC', scope='spin', default=0.5, desc='The population for state C', set='params', py_type=float, grace_string='\\qp\\sC\\N\\Q', err=True, sim=True) 91 self.PARAMS.add('phi_ex', scope='spin', default=5.0, desc='The phi_ex = pA.pB.dw**2 value (ppm^2)', set='params', py_type=float, grace_string='\\xF\\B\\sex\\N = \\q p\\sA\\N.p\\sB\\N.\\xDw\\B\\S2\\N\\Q (ppm\\S2\\N)', err=True, sim=True) 92 self.PARAMS.add('phi_ex_B', scope='spin', default=5.0, desc='The fast exchange factor between sites A and B (ppm^2)', set='params', py_type=float, grace_string='\\xF\\B\\sex,B\\N (ppm\\S2\\N)', err=True, sim=True) 93 self.PARAMS.add('phi_ex_C', scope='spin', default=5.0, desc='The fast exchange factor between sites A and C (ppm^2)', set='params', py_type=float, grace_string='\\xF\\B\\sex,C\\N (ppm\\S2\\N)', err=True, sim=True) 94 self.PARAMS.add('padw2', scope='spin', default=1.0, desc='The pA.dw**2 value (ppm^2)', set='params', py_type=float, grace_string='\\qp\\sA\\N.\\xDw\\B\\S2\\N\\Q (ppm\\S2\\N)', err=True, sim=True) 95 self.PARAMS.add('dw', scope='spin', default=0.0, desc='The chemical shift difference between states A and B (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\Q (ppm)', err=True, sim=True) 96 self.PARAMS.add('dw_AB', scope='spin', default=0.0, desc='The chemical shift difference between states A and B for 3-site exchange (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\Q\\SAB\\N (ppm)', err=True, sim=True) 97 self.PARAMS.add('dw_AC', scope='spin', default=0.0, desc='The chemical shift difference between states A and C for 3-site exchange (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\Q\\SAC\\N (ppm)', err=True, sim=True) 98 self.PARAMS.add('dw_BC', scope='spin', default=0.0, desc='The chemical shift difference between states B and C for 3-site exchange (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\Q\\SBC\\N (ppm)', err=True, sim=True) 99 self.PARAMS.add('dwH', scope='spin', default=0.0, desc='The proton chemical shift difference between states A and B (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\sH\\N\\Q (ppm)', err=True, sim=True) 100 self.PARAMS.add('dwH_AB', scope='spin', default=0.0, desc='The proton chemical shift difference between states A and B for 3-site exchange (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\sH\\N\\Q\\SAB\\N (ppm)', err=True, sim=True) 101 self.PARAMS.add('dwH_AC', scope='spin', default=0.0, desc='The proton chemical shift difference between states A and C for 3-site exchange (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\sH\\N\\Q\\SAC\\N (ppm)', err=True, sim=True) 102 self.PARAMS.add('dwH_BC', scope='spin', default=0.0, desc='The proton chemical shift difference between states B and C for 3-site exchange (in ppm)', set='params', py_type=float, grace_string='\\q\\xDw\\B\\sH\\N\\Q\\SBC\\N (ppm)', err=True, sim=True) 103 self.PARAMS.add('kex', scope='spin', default=10000.0, desc='The exchange rate', set='params', py_type=float, grace_string='\\qk\\sex\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 104 self.PARAMS.add('kex_AB', scope='spin', default=10000.0, desc='The exchange rate between sites A and B for 3-site exchange with kex_AB = k_AB + k_BA (rad.s^-1)', set='params', py_type=float, grace_string='\\qk\\sex\\N\\Q\\SAB\\N (rad.s\\S-1\\N)', err=True, sim=True) 105 self.PARAMS.add('kex_AC', scope='spin', default=10000.0, desc='The exchange rate between sites A and C for 3-site exchange with kex_AC = k_AC + k_CA (rad.s^-1)', set='params', py_type=float, grace_string='\\qk\\sex\\N\\Q\\SAC\\N (rad.s\\S-1\\N)', err=True, sim=True) 106 self.PARAMS.add('kex_BC', scope='spin', default=10000.0, desc='The exchange rate between sites B and C for 3-site exchange with kex_BC = k_BC + k_CB (rad.s^-1)', set='params', py_type=float, grace_string='\\qk\\sex\\N\\Q\\SBC\\N (rad.s\\S-1\\N)', err=True, sim=True) 107 self.PARAMS.add('kB', scope='spin', default=10000.0, desc='Approximate chemical exchange rate constant between sites A and B (rad.s^-1)', set='params', py_type=float, grace_string='\\qk\\sB\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 108 self.PARAMS.add('kC', scope='spin', default=10000.0, desc='Approximate chemical exchange rate constant between sites A and C (rad.s^-1)', set='params', py_type=float, grace_string='\\qk\\sC\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 109 self.PARAMS.add('tex', scope='spin', default=1.0/10000.0, desc='The time of exchange (tex = 1/kex)', set='params', py_type=float, grace_string='\\q\\xt\\B\\sex\\N\\Q (s.rad\\S-1\\N)', err=True, sim=True) 110 self.PARAMS.add('theta', scope='spin', desc='Rotating frame tilt angle : ( theta = arctan(w_1 / Omega) ) (rad)', set='params', grace_string='Rotating frame tilt angle (rad)', py_type=dict, err=False, sim=False) 111 self.PARAMS.add('w_eff', scope='spin', desc='Effective field in rotating frame : ( w_eff = sqrt(Omega^2 + w_1^2) ) (rad.s^-1)', grace_string='Effective field in rotating frame (rad.s\\S-1\\N)', set='params', py_type=dict, err=False, sim=False) 112 self.PARAMS.add('k_AB', scope='spin', default=10000.0, desc='The exchange rate from state A to state B', set='params', py_type=float, grace_string='\\qk\\sAB\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 113 self.PARAMS.add('k_BA', scope='spin', default=10000.0, desc='The exchange rate from state B to state A', set='params', py_type=float, grace_string='\\qk\\sBA\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True) 114 self.PARAMS.add('params', scope='spin', desc='The model parameters', py_type=list) 115 self.PARAMS.add('model', scope='spin', desc='The dispersion model', py_type=str) 116 117 # Add the minimisation data. 118 self.PARAMS.add_min_data(min_stats_global=False, min_stats_spin=True)

119 120

121 - def _back_calc_peak_intensities(self, spin=None, exp_type=None, frq=None, offset=None, point=None):

122 """Back-calculation of peak intensity for the given relaxation time. 123 124 @keyword spin: The specific spin data container. 125 @type spin: SpinContainer instance 126 @keyword exp_type: The experiment type. 127 @type exp_type: str 128 @keyword frq: The spectrometer frequency. 129 @type frq: float 130 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm. 131 @type offset: None or float 132 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 133 @type point: float 134 @return: The back-calculated peak intensities for the given exponential curve. 135 @rtype: numpy rank-1 float array 136 """ 137 138 # Check. 139 if not has_exponential_exp_type(): 140 raise RelaxError("Back-calculation is not allowed for the fixed time experiment types.") 141 142 # The key. 143 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 144 145 # Create the initial parameter vector. 146 param_vector = assemble_param_vector(spins=[spin], key=param_key) 147 148 # Create a scaling matrix. 149 scaling_matrix = assemble_scaling_matrix(spins=[spin], key=param_key, scaling=False) 150 151 # The peak intensities and times. 152 values = [] 153 errors = [] 154 times = [] 155 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 156 # The data. 157 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) 158 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 159 times.append(time) 160 161 # The scaling matrix in a diagonalised list form. 162 scaling_list = [] 163 for i in range(len(scaling_matrix)): 164 scaling_list.append(scaling_matrix[i, i]) 165 166 # Initialise the relaxation fit functions. 167 setup(num_params=len(param_vector), num_times=len(times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list) 168 169 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 170 func(param_vector) 171 172 # Get the data back. 173 results = back_calc_I() 174 175 # Return the correct peak height. 176 return results

177 178

179 - def _calculate_r2eff(self):

180 """Calculate the R2eff values for fixed relaxation time period data.""" 181 182 # Data checks. 183 check_exp_type() 184 check_disp_points() 185 check_exp_type_fixed_time() 186 187 # Printouts. 188 print("Calculating the R2eff/R1rho values for fixed relaxation time period data.") 189 190 # Loop over the spins. 191 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 192 # Spin ID printout. 193 print("Spin '%s'." % spin_id) 194 195 # Skip spins which have no data. 196 if not hasattr(spin, 'intensities'): 197 continue 198 199 # Initialise the data structures. 200 if not hasattr(spin, 'r2eff'): 201 spin.r2eff = {} 202 if not hasattr(spin, 'r2eff_err'): 203 spin.r2eff_err = {} 204 205 # Loop over all the data. 206 for exp_type, frq, offset, point, time in loop_exp_frq_offset_point_time(): 207 # The three keys. 208 ref_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=None, time=time) 209 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 210 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 211 212 # Check for missing data. 213 missing = False 214 for i in range(len(ref_keys)): 215 if ref_keys[i] not in spin.intensities: 216 missing = True 217 for i in range(len(int_keys)): 218 if int_keys[i] not in spin.intensities: 219 missing = True 220 if missing: 221 continue 222 223 # Average the reference intensity data and errors. 224 ref_intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time) 225 ref_intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time, error=True) 226 227 # Average the intensity data and errors. 228 intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 229 intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True) 230 231 # Calculate the R2eff value. 232 spin.r2eff[param_key] = calc_two_point_r2eff(relax_time=time, I_ref=ref_intensity, I=intensity) 233 234 # Calculate the R2eff error. 235 spin.r2eff_err[param_key] = calc_two_point_r2eff_err(relax_time=time, I_ref=ref_intensity, I=intensity, I_ref_err=ref_intensity_err, I_err=intensity_err)

236 237

238 - def _cluster(self, cluster_id=None, spin_id=None):

239 """Define spin clustering. 240 241 @keyword cluster_id: The cluster ID string. 242 @type cluster_id: str 243 @keyword spin_id: The spin ID string for the spin or group of spins to add to the cluster. 244 @type spin_id: str 245 """ 246 247 # Initialise. 248 if not hasattr(cdp, 'clustering'): 249 # Create the dictionary. 250 cdp.clustering = {} 251 cdp.clustering['free spins'] = [] 252 253 # Add all spin IDs to the cluster. 254 for spin, id in spin_loop(return_id=True): 255 cdp.clustering['free spins'].append(id) 256 257 # Add the key. 258 if cluster_id not in cdp.clustering: 259 cdp.clustering[cluster_id] = [] 260 261 # Loop over the spins to add to the cluster. 262 for spin, id in spin_loop(selection=spin_id, return_id=True): 263 # First remove the ID from all clusters. 264 for key in cdp.clustering.keys(): 265 if id in cdp.clustering[key]: 266 cdp.clustering[key].pop(cdp.clustering[key].index(id)) 267 268 # Then add the ID to the cluster. 269 cdp.clustering[cluster_id].append(id) 270 271 # Clean up - delete any empty clusters (except the free spins). 272 clean = [] 273 for key in cdp.clustering.keys(): 274 if key == 'free spins': 275 continue 276 if cdp.clustering[key] == []: 277 clean.append(key) 278 for key in clean: 279 cdp.clustering.pop(key)

280 281

282 - def _cluster_ids(self):

283 """Return the current list of cluster ID strings. 284 285 @return: The list of cluster IDs. 286 @rtype: list of str 287 """ 288 289 # Initialise. 290 ids = ['free spins'] 291 292 # Add the defined IDs. 293 if hasattr(cdp, 'clustering'): 294 for key in list(cdp.clustering.keys()): 295 if key not in ids: 296 ids.append(key) 297 298 # Return the IDs. 299 return ids

300 301

302 - def _minimise_r2eff(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

303 """Optimise the R2eff model by fitting the 2-parameter exponential curves. 304 305 This mimics the R1 and R2 relax_fit analysis. 306 307 308 @keyword min_algor: The minimisation algorithm to use. 309 @type min_algor: str 310 @keyword min_options: An array of options to be used by the minimisation algorithm. 311 @type min_options: array of str 312 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 313 @type func_tol: None or float 314 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 315 @type grad_tol: None or float 316 @keyword max_iterations: The maximum number of iterations for the algorithm. 317 @type max_iterations: int 318 @keyword constraints: If True, constraints are used during optimisation. 319 @type constraints: bool 320 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 321 @type scaling: bool 322 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 323 @type verbosity: int 324 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 325 @type sim_index: None or int 326 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 327 @type lower: array of numbers 328 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 329 @type upper: array of numbers 330 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search. 331 @type inc: array of int 332 """ 333 334 # Check that the C modules have been compiled. 335 if not C_module_exp_fn: 336 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 337 338 # Loop over the spins. 339 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 340 # Skip spins which have no data. 341 if not hasattr(spin, 'intensities'): 342 continue 343 344 # Loop over each spectrometer frequency and dispersion point. 345 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 346 # The parameter key. 347 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 348 349 # The initial parameter vector. 350 param_vector = assemble_param_vector(spins=[spin], key=param_key, sim_index=sim_index) 351 352 # Diagonal scaling. 353 scaling_matrix = assemble_scaling_matrix(spins=[spin], key=param_key, scaling=scaling) 354 if len(scaling_matrix): 355 param_vector = dot(inv(scaling_matrix), param_vector) 356 357 # Get the grid search minimisation options. 358 lower_new, upper_new = None, None 359 if match('^[Gg]rid', min_algor): 360 grid_size, inc_new, lower_new, upper_new = grid_search_setup(spins=[spin], spin_ids=[spin_id], param_vector=param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix) 361 362 # Linear constraints. 363 A, b = None, None 364 if constraints: 365 A, b = linear_constraints(spins=[spin], scaling_matrix=scaling_matrix) 366 367 # Print out. 368 if verbosity >= 1: 369 # Individual spin section. 370 top = 2 371 if verbosity >= 2: 372 top += 2 373 text = "Fitting to spin %s, frequency %s and dispersion point %s" % (spin_id, frq, point) 374 subsection(file=sys.stdout, text=text, prespace=top) 375 376 # Grid search printout. 377 if match('^[Gg]rid', min_algor): 378 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % grid_size) 379 380 # The peak intensities, errors and times. 381 values = [] 382 errors = [] 383 times = [] 384 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 385 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, sim_index=sim_index)) 386 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 387 times.append(time) 388 389 # The scaling matrix in a diagonalised list form. 390 scaling_list = [] 391 for i in range(len(scaling_matrix)): 392 scaling_list.append(scaling_matrix[i, i]) 393 394 # Initialise the function to minimise. 395 setup(num_params=len(param_vector), num_times=len(times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list) 396 397 # Grid search. 398 if search('^[Gg]rid', min_algor): 399 results = grid(func=func, args=(), num_incs=inc_new, lower=lower_new, upper=upper_new, A=A, b=b, verbosity=verbosity) 400 401 # Unpack the results. 402 param_vector, chi2, iter_count, warning = results 403 f_count = iter_count 404 g_count = 0.0 405 h_count = 0.0 406 407 # Minimisation. 408 else: 409 results = generic_minimise(func=func, dfunc=dfunc, d2func=d2func, args=(), x0=param_vector, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, maxiter=max_iterations, A=A, b=b, full_output=True, print_flag=verbosity) 410 411 # Unpack the results. 412 if results == None: 413 return 414 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 415 416 # Scaling. 417 if scaling: 418 param_vector = dot(scaling_matrix, param_vector) 419 420 # Disassemble the parameter vector. 421 disassemble_param_vector(param_vector=param_vector, spins=[spin], key=param_key, sim_index=sim_index) 422 423 # Monte Carlo minimisation statistics. 424 if sim_index != None: 425 # Chi-squared statistic. 426 spin.chi2_sim[sim_index] = chi2 427 428 # Iterations. 429 spin.iter_sim[sim_index] = iter_count 430 431 # Function evaluations. 432 spin.f_count_sim[sim_index] = f_count 433 434 # Gradient evaluations. 435 spin.g_count_sim[sim_index] = g_count 436 437 # Hessian evaluations. 438 spin.h_count_sim[sim_index] = h_count 439 440 # Warning. 441 spin.warning_sim[sim_index] = warning 442 443 # Normal statistics. 444 else: 445 # Chi-squared statistic. 446 spin.chi2 = chi2 447 448 # Iterations. 449 spin.iter = iter_count 450 451 # Function evaluations. 452 spin.f_count = f_count 453 454 # Gradient evaluations. 455 spin.g_count = g_count 456 457 # Hessian evaluations. 458 spin.h_count = h_count 459 460 # Warning. 461 spin.warning = warning

462 463

464 - def _model_setup(self, model, params):

465 """Update various model specific data structures. 466 467 @param model: The relaxation dispersion curve type. 468 @type model: str 469 @param params: A list consisting of the model parameters. 470 @type params: list of str 471 """ 472 473 # The model group. 474 if model == MODEL_R2EFF: 475 cdp.model_type = 'R2eff' 476 else: 477 cdp.model_type = 'disp' 478 479 # Loop over the sequence. 480 for spin in spin_loop(skip_desel=True): 481 # The model and parameter names. 482 spin.model = model 483 spin.params = params 484 485 # Initialise the data structures (if needed). 486 self.data_init(spin)

487 488

489 - def _select_model(self, model=MODEL_R2EFF):

490 """Set up the model for the relaxation dispersion analysis. 491 492 @keyword model: The relaxation dispersion analysis type. 493 @type model: str 494 """ 495 496 # Data checks. 497 pipes.test() 498 check_pipe_type() 499 check_mol_res_spin_data() 500 check_exp_type() 501 if model == MODEL_R2EFF: 502 check_c_modules() 503 504 # The curve type. 505 curve_type = get_curve_type() 506 507 # R2eff/R1rho model. 508 if model == MODEL_R2EFF: 509 print("R2eff/R1rho value and error determination.") 510 if curve_type == 'exponential': 511 params = ['r2eff', 'i0'] 512 else: 513 params = ['r2eff'] 514 515 # The model for no chemical exchange relaxation. 516 elif model == MODEL_NOREX: 517 print("The model for no chemical exchange relaxation.") 518 params = ['r2'] 519 520 # LM63 model. 521 elif model == MODEL_LM63: 522 print("The Luz and Meiboom (1963) 2-site fast exchange model.") 523 params = ['r2', 'phi_ex', 'kex'] 524 525 # LM63 3-site model. 526 elif model == MODEL_LM63_3SITE: 527 print("The Luz and Meiboom (1963) 3-site fast exchange model.") 528 params = ['r2', 'phi_ex_B', 'phi_ex_C', 'kB', 'kC'] 529 530 # Full CR72 model. 531 elif model == MODEL_CR72_FULL: 532 print("The full Carver and Richards (1972) 2-site model for all time scales.") 533 params = ['r2a', 'r2b', 'pA', 'dw', 'kex'] 534 535 # Reduced CR72 model. 536 elif model == MODEL_CR72: 537 print("The reduced Carver and Richards (1972) 2-site model for all time scales, whereby the simplification R20A = R20B is assumed.") 538 params = ['r2', 'pA', 'dw', 'kex'] 539 540 # IT99 model. 541 elif model == MODEL_IT99: 542 print("The Ishima and Torchia (1999) CPMG 2-site model for all time scales with pA >> pB.") 543 params = ['r2', 'pA', 'dw', 'tex'] 544 545 # TSMFK01 model. 546 elif model == MODEL_TSMFK01: 547 print("The Tollinger et al. (2001) 2-site very-slow exchange model, range of microsecond to second time scale.") 548 params = ['r2a', 'dw', 'k_AB'] 549 550 # Full NS CPMG 2-site 3D model. 551 elif model == MODEL_NS_CPMG_2SITE_3D_FULL: 552 print("The full numerical solution for the 2-site Bloch-McConnell equations for CPMG data using 3D magnetisation vectors.") 553 params = ['r2a', 'r2b', 'pA', 'dw', 'kex'] 554 555 # Reduced NS CPMG 2-site 3D model. 556 elif model == MODEL_NS_CPMG_2SITE_3D: 557 print("The reduced numerical solution for the 2-site Bloch-McConnell equations for CPMG data using 3D magnetisation vectors, whereby the simplification R20A = R20B is assumed.") 558 params = ['r2', 'pA', 'dw', 'kex'] 559 560 # NS CPMG 2-site expanded model. 561 elif model == MODEL_NS_CPMG_2SITE_EXPANDED: 562 print("The numerical solution for the 2-site Bloch-McConnell equations for CPMG data expanded using Maple by Nikolai Skrynnikov.") 563 params = ['r2', 'pA', 'dw', 'kex'] 564 565 # Full NS CPMG 2-site star model. 566 elif model == MODEL_NS_CPMG_2SITE_STAR_FULL: 567 print("The full numerical solution for the 2-site Bloch-McConnell equations for CPMG data using complex conjugate matrices.") 568 params = ['r2a', 'r2b', 'pA', 'dw', 'kex'] 569 570 # Reduced NS CPMG 2-site star model. 571 elif model == MODEL_NS_CPMG_2SITE_STAR: 572 print("The numerical reduced solution for the 2-site Bloch-McConnell equations for CPMG data using complex conjugate matrices, whereby the simplification R20A = R20B is assumed.") 573 params = ['r2', 'pA', 'dw', 'kex'] 574 575 # M61 model. 576 elif model == MODEL_M61: 577 print("The Meiboom (1961) 2-site fast exchange model for R1rho-type experiments.") 578 params = ['r2', 'phi_ex', 'kex'] 579 580 # M61 skew model. 581 elif model == MODEL_M61B: 582 print("The Meiboom (1961) on-resonance 2-site model with skewed populations (pA >> pB) for R1rho-type experiments.") 583 params = ['r2', 'pA', 'dw', 'kex'] 584 585 # DPL94 model. 586 elif model == MODEL_DPL94: 587 print("The Davis, Perlman and London (1994) 2-site fast exchange model for R1rho-type experiments.") 588 params = ['r2', 'phi_ex', 'kex'] 589 590 # TP02 model. 591 elif model == MODEL_TP02: 592 print("The Trott and Palmer (2002) off-resonance 2-site model for R1rho-type experiments.") 593 params = ['r2', 'pA', 'dw', 'kex'] 594 595 # TAP03 model. 596 elif model == MODEL_TAP03: 597 print("The Trott, Abergel and Palmer (2003) off-resonance 2-site model for R1rho-type experiments.") 598 params = ['r2', 'pA', 'dw', 'kex'] 599 600 # MP05 model. 601 elif model == MODEL_MP05: 602 print("The Miloushev and Palmer (2005) off-resonance 2-site model for R1rho-type experiments.") 603 params = ['r2', 'pA', 'dw', 'kex'] 604 605 # Reduced NS R1rho 2-site model. 606 elif model == MODEL_NS_R1RHO_2SITE: 607 print("The reduced numerical solution for the 2-site Bloch-McConnell equations for R1rho data using 3D magnetisation vectors, whereby the simplification R20A = R20B is assumed.") 608 params = ['r2', 'pA', 'dw', 'kex'] 609 610 # NS R1rho CPMG 3-site model. 611 elif model == MODEL_NS_R1RHO_3SITE: 612 print("The numerical solution for the 3-site Bloch-McConnell equations for R1rho data using 3D magnetisation vectors whereby the simplification R20A = R20B = R20C is assumed.") 613 params = ['r2', 'pA', 'dw_AB', 'kex_AB', 'pB', 'dw_BC', 'kex_BC', 'kex_AC'] 614 615 # NS R1rho CPMG 3-site linearised model. 616 elif model == MODEL_NS_R1RHO_3SITE_LINEAR: 617 print("The numerical solution for the 3-site Bloch-McConnell equations for R1rho data using 3D magnetisation vectors linearised with kAC = kCA = 0 whereby the simplification R20A = R20B = R20C is assumed.") 618 params = ['r2', 'pA', 'dw_AB', 'kex_AB', 'pB', 'dw_BC', 'kex_BC'] 619 620 # MMQ CR72 model. 621 elif model == MODEL_MMQ_CR72: 622 print("The Carver and Richards (1972) 2-site model for all time scales expanded for MMQ CPMG data by Korzhnev et al., 2004.") 623 params = ['r2', 'pA', 'dw', 'dwH', 'kex'] 624 625 # NS MQ CPMG 2-site model. 626 elif model == MODEL_NS_MMQ_2SITE: 627 print("The reduced numerical solution for the 2-site Bloch-McConnell equations for MQ CPMG data using 3D magnetisation vectors, whereby the simplification R20A = R20B is assumed.") 628 params = ['r2', 'pA', 'dw', 'dwH', 'kex'] 629 630 # NS MMQ CPMG 3-site model. 631 elif model == MODEL_NS_MMQ_3SITE: 632 print("The numerical solution for the 3-site Bloch-McConnell equations for MMQ CPMG data whereby the simplification R20A = R20B = R20C is assumed.") 633 params = ['r2', 'pA', 'dw_AB', 'dwH_AB', 'kex_AB', 'pB', 'dw_BC', 'dwH_BC', 'kex_BC', 'kex_AC'] 634 635 # NS MMQ CPMG 3-site linearised model. 636 elif model == MODEL_NS_MMQ_3SITE_LINEAR: 637 print("The numerical solution for the 3-site Bloch-McConnell equations for MMQ CPMG data linearised with kAC = kCA = 0 whereby the simplification R20A = R20B = R20C is assumed.") 638 params = ['r2', 'pA', 'dw_AB', 'dwH_AB', 'kex_AB', 'pB', 'dw_BC', 'dwH_BC', 'kex_BC'] 639 640 # Invalid model. 641 else: 642 raise RelaxError("The model '%s' must be one of %s." % (model, MODEL_LIST_FULL)) 643 644 # Set up the model. 645 self._model_setup(model, params)

646 647

648 - def base_data_loop(self):

649 """Custom generator method for looping over the base data. 650 651 For the R2eff model, the base data is the peak intensity data defining a single exponential curve. For all other models, the base data is the R2eff/R1rho values for individual spins. 652 653 654 @return: For the R2eff model, a tuple of the spin container and the exponential curve identifying key (the CPMG frequency or R1rho spin-lock field strength). For all other models, the spin container and ID from the spin loop. 655 @rtype: (tuple of SpinContainer instance and float) or (SpinContainer instance and str) 656 """ 657 658 # The R2eff model data (the base data is peak intensities). 659 if cdp.model_type == 'R2eff': 660 # Loop over the sequence. 661 for spin in spin_loop(): 662 # Skip deselected spins. 663 if not spin.select: 664 continue 665 666 # Skip spins with no peak intensity data. 667 if not hasattr(spin, 'intensities'): 668 continue 669 670 # Loop over each spectrometer frequency and dispersion point. 671 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 672 yield spin, exp_type, frq, offset, point 673 674 # All other models (the base data is the R2eff/R1rho values). 675 else: 676 # 1H MMQ flag. 677 proton_mmq_flag = has_proton_mmq_cpmg() 678 679 # Loop over the sequence. 680 for spin, spin_id in spin_loop(return_id=True): 681 # Skip deselected spins. 682 if not spin.select: 683 continue 684 685 # Skip protons for MMQ data. 686 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 687 continue 688 689 # Get the attached proton. 690 proton = None 691 if proton_mmq_flag: 692 proton = return_attached_protons(spin_id)[0] 693 694 # Skip spins with no R2eff/R1rho values. 695 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'): 696 continue 697 698 # Yield the spin container and ID. 699 yield spin, spin_id

700 701

702 - def calculate(self, spin_id=None, verbosity=1, sim_index=None):

703 """Calculate the model chi-squared value or the R2eff values for fixed time period data. 704 705 @keyword spin_id: The spin identification string. 706 @type spin_id: None or str 707 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 708 @type verbosity: int 709 @keyword sim_index: The MC simulation index (unused). 710 @type sim_index: None 711 """ 712 713 # Data checks. 714 pipes.test() 715 check_mol_res_spin_data() 716 check_model_type() 717 718 # Special exponential curve-fitting for the 'R2eff' model. 719 if cdp.model_type == 'R2eff': 720 self._calculate_r2eff() 721 722 # Calculate the chi-squared value. 723 else: 724 # 1H MMQ flag. 725 proton_mmq_flag = has_proton_mmq_cpmg() 726 727 # Loop over all spins. 728 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 729 # Skip protons for MMQ data. 730 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 731 continue 732 733 # Get the attached proton. 734 proton = None 735 if proton_mmq_flag: 736 proton = return_attached_protons(spin_id)[0] 737 738 # The back calculated values. 739 back_calc = back_calc_r2eff(spin=spin, spin_id=spin_id, store_chi2=True) 740 741 # Pack the data. 742 pack_back_calc_r2eff(spin=spin, spin_id=spin_id, si=0, back_calc=back_calc, proton_mmq_flag=proton_mmq_flag)

743 744

745 - def constraint_algorithm(self):

746 """Return the 'Log barrier' optimisation constraint algorithm. 747 748 @return: The 'Log barrier' constraint algorithm. 749 @rtype: str 750 """ 751 752 # The log barrier algorithm, as required by minfx. 753 return 'Log barrier'

754 755

756 - def create_mc_data(self, data_id):

757 """Create the Monte Carlo peak intensity data. 758 759 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 760 @type data_id: SpinContainer instance and float 761 @return: The Monte Carlo simulation data. 762 @rtype: list of floats 763 """ 764 765 # The R2eff model (with peak intensity base data). 766 if cdp.model_type == 'R2eff': 767 # Unpack the data. 768 spin, exp_type, frq, offset, point = data_id 769 770 # Back calculate the peak intensities. 771 values = self._back_calc_peak_intensities(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point) 772 773 # All other models (with R2eff/R1rho base data). 774 else: 775 # 1H MMQ flag. 776 proton_mmq_flag = has_proton_mmq_cpmg() 777 778 # Unpack the data. 779 spin, spin_id = data_id 780 781 # Back calculate the R2eff/R1rho data. 782 back_calc = back_calc_r2eff(spin=spin, spin_id=spin_id) 783 784 # Get the attached proton data. 785 if proton_mmq_flag: 786 proton = return_attached_protons(spin_id)[0] 787 proton_back_calc = back_calc_r2eff(spin=proton, spin_id=spin_id) 788 789 # Convert to a dictionary matching the R2eff data structure. 790 values = {} 791 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 792 # Alias the correct data. 793 current_bc = back_calc 794 current_spin = spin 795 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 796 current_spin = proton 797 current_bc = proton_back_calc 798 799 # The parameter key. 800 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 801 802 # Skip missing data. 803 if not hasattr(current_spin, 'r2eff') or param_key not in current_spin.r2eff.keys(): 804 continue 805 806 # Store the result. 807 values[param_key] = back_calc[ei][0][mi][oi][di] 808 809 # Return the MC data. 810 return values

811 812 813 default_value_doc = Desc_container("Relaxation dispersion default values") 814 _table = uf_tables.add_table(label="table: dispersion default values", caption="Relaxation dispersion default values.") 815 _table.add_headings(["Data type", "Object name", "Value"]) 816 _table.add_row(["Transversal relaxation rate (rad/s)", "'r2'", "15.0"]) 817 _table.add_row(["Transversal relaxation rate for state A (rad/s)", "'r2a'", "15.0"]) 818 _table.add_row(["Transversal relaxation rate for state B (rad/s)", "'r2b'", "15.0"]) 819 _table.add_row(["Population of state A", "'pA'", "0.5"]) 820 _table.add_row(["Population of state B", "'pB'", "0.5"]) 821 _table.add_row(["Population of state C", "'pC'", "0.5"]) 822 _table.add_row(["The pA.pB.dw**2 parameter (ppm^2)", "'phi_ex'", "5.0"]) 823 _table.add_row(["The pA.pB.dw**2 parameter of state B (ppm^2)", "'phi_ex_B'", "5.0"]) 824 _table.add_row(["The pA.pB.dw**2 parameter of state C (ppm^2)", "'phi_ex_C'", "5.0"]) 825 _table.add_row(["The pA.dw**2 parameter (ppm^2)", "'padw2'", "1.0"]) 826 _table.add_row(["Chemical shift difference between states A and B (ppm)", "'dw'", "0.0"]) 827 _table.add_row(["Chemical shift difference between states A and B for 3-site exchange (ppm)", "'dw_AB'", "0.0"]) 828 _table.add_row(["Chemical shift difference between states A and C for 3-site exchange (ppm)", "'dw_AC'", "0.0"]) 829 _table.add_row(["Chemical shift difference between states B and C for 3-site exchange (ppm)", "'dw_BC'", "0.0"]) 830 _table.add_row(["Proton chemical shift difference between states A and B (ppm)", "'dwH'", "0.0"]) 831 _table.add_row(["Proton chemical shift difference between states A and B for 3-site exchange (ppm)", "'dwH_AB'", "0.0"]) 832 _table.add_row(["Proton chemical shift difference between states A and C for 3-site exchange (ppm)", "'dwH_AC'", "0.0"]) 833 _table.add_row(["Proton chemical shift difference between states B and C for 3-site exchange (ppm)", "'dwH_BC'", "0.0"]) 834 _table.add_row(["Exchange rate (rad/s)", "'kex'", "10000.0"]) 835 _table.add_row(["Exchange rate between sites A and B for 3-site exchange (rad/s)", "'kex_AB'", "10000.0"]) 836 _table.add_row(["Exchange rate between sites A and C for 3-site exchange (rad/s)", "'kex_AC'", "10000.0"]) 837 _table.add_row(["Exchange rate between sites B and C for 3-site exchange (rad/s)", "'kex_BC'", "10000.0"]) 838 _table.add_row(["Exchange rate between sites A and B (rad/s)", "'kB'", "10000.0"]) 839 _table.add_row(["Exchange rate between sites A and C (rad/s)", "'kC'", "10000.0"]) 840 _table.add_row(["Time of exchange (s/rad)", "'tex'", "1.0/10000.0"]) 841 _table.add_row(["Rotating frame tilt angle", "'theta'", "0.0"]) 842 _table.add_row(["Effective field in rotating frame", "'w_eff'", "0.0"]) 843 _table.add_row(["Exchange rate from state A to state B (rad/s)", "'k_AB'", "10000.0"]) 844 _table.add_row(["Exchange rate from state B to state A (rad/s)", "'k_BA'", "10000.0"]) 845 default_value_doc.add_table(_table.label) 846 847

848 - def deselect(self, model_info, sim_index=None):

849 """Deselect models or simulations. 850 851 @param model_info: The spin ID list from the model_loop() API method. 852 @type model_info: int 853 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will. 854 @type sim_index: None or int 855 """ 856 857 # Loop over all the spins and deselect them. 858 for spin_id in model_info: 859 # Get the spin. 860 spin = return_spin(spin_id) 861 862 # Spin deselection. 863 if sim_index == None: 864 spin.select = False 865 866 # Simulation deselection. 867 else: 868 spin.select_sim[sim_index] = False

869 870

871 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):

872 """Duplicate the data specific to a single model. 873 874 @keyword pipe_from: The data pipe to copy the data from. 875 @type pipe_from: str 876 @keyword pipe_to: The data pipe to copy the data to. 877 @type pipe_to: str 878 @keyword model_info: The model index from model_info(). 879 @type model_info: int 880 @keyword global_stats: The global statistics flag. 881 @type global_stats: bool 882 @keyword verbose: A flag which if True will cause info to be printed out. 883 @type verbose: bool 884 """ 885 886 # First create the pipe_to data pipe, if it doesn't exist, but don't switch to it. 887 if not pipes.has_pipe(pipe_to): 888 pipes.create(pipe_to, pipe_type='relax_disp', switch=False) 889 890 # Get the data pipes. 891 dp_from = pipes.get_pipe(pipe_from) 892 dp_to = pipes.get_pipe(pipe_to) 893 894 # Duplicate all non-sequence specific data. 895 for data_name in dir(dp_from): 896 # Skip the container objects. 897 if data_name in ['mol', 'interatomic', 'structure', 'exp_info', 'result_files']: 898 continue 899 900 # Skip dispersion specific parameters. 901 if data_name in ['model']: 902 continue 903 904 # Skip special objects. 905 if search('^_', data_name) or data_name in list(dp_from.__class__.__dict__.keys()): 906 continue 907 908 # Get the original object. 909 data_from = getattr(dp_from, data_name) 910 911 # The data already exists. 912 if hasattr(dp_to, data_name): 913 # Get the object in the target pipe. 914 data_to = getattr(dp_to, data_name) 915 916 # The data must match! 917 if data_from != data_to: 918 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 919 920 # Skip the data. 921 continue 922 923 # Duplicate the data. 924 setattr(dp_to, data_name, deepcopy(data_from)) 925 926 # No sequence data, so skip the rest. 927 if dp_from.mol.is_empty(): 928 return 929 930 # Duplicate the sequence data if it doesn't exist. 931 if dp_to.mol.is_empty(): 932 sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose) 933 934 # Loop over the cluster. 935 for id in model_info: 936 # The original spin container. 937 spin = return_spin(id, pipe=pipe_from) 938 939 # Duplicate the spin specific data. 940 for name in dir(spin): 941 # Skip special objects. 942 if search('^__', name): 943 continue 944 945 # Get the object. 946 obj = getattr(spin, name) 947 948 # Skip methods. 949 if isinstance(obj, MethodType): 950 continue 951 952 # Duplicate the object. 953 new_obj = deepcopy(getattr(spin, name)) 954 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj)

955 956

957 - def eliminate(self, name, value, model_info, args, sim=None):

958 """Relaxation dispersion model elimination, parameter by parameter. 959 960 @param name: The parameter name. 961 @type name: str 962 @param value: The parameter value. 963 @type value: float 964 @param model_info: The list of spin IDs from the model_loop() API method. 965 @type model_info: int 966 @param args: The c1 and c2 elimination constant overrides. 967 @type args: None or tuple of float 968 @keyword sim: The Monte Carlo simulation index. 969 @type sim: int 970 @return: True if the model is to be eliminated, False otherwise. 971 @rtype: bool 972 """ 973 974 # Skip the R2eff model parameters. 975 if name in ['r2eff', 'i0']: 976 return False 977 978 # Default limits. 979 c1 = 0.501 980 c2 = 0.999 981 c3 = 1.0 982 983 # Depack the arguments. 984 if args != None: 985 c1, c2, c3 = args 986 987 # Elimination text. 988 elim_text = "Data pipe '%s': The %s parameter of %.5f is %s, eliminating " % (pipes.cdp_name(), name, value, "%s") 989 if sim == None: 990 elim_text += "the spin cluster %s." % model_info 991 else: 992 elim_text += "simulation %i of the spin cluster %s." % (sim, model_info) 993 994 # The pA parameter. 995 if name == 'pA': 996 if value < c1: 997 print(elim_text % ("less than %.5f" % c1)) 998 return True 999 if value > c2: 1000 print(elim_text % ("greater than %.5f" % c2)) 1001 return True 1002 1003 # The tex parameter. 1004 if name == 'tex': 1005 if value > c3: 1006 print(elim_text % ("greater than %.5f" % c3)) 1007 return True 1008 1009 # Accept model. 1010 return False

1011 1012

1013 - def get_param_names(self, model_info=None):

1014 """Return a vector of parameter names. 1015 1016 @keyword model_info: The list spin ID strings from the model_loop() API method. 1017 @type model_info: list of str 1018 @return: The vector of parameter names. 1019 @rtype: list of str 1020 """ 1021 1022 # Get the spin containers. 1023 spins = [] 1024 for spin_id in model_info: 1025 # Get the spin. 1026 spin = return_spin(spin_id) 1027 1028 # Skip deselected spins. 1029 if not spin.select: 1030 continue 1031 1032 # Add the spin. 1033 spins.append(spin) 1034 1035 # No spins. 1036 if not len(spins): 1037 return None 1038 1039 # Call the get_param_names() function. 1040 return get_param_names(spins)

1041 1042

1043 - def get_param_values(self, model_info=None, sim_index=None):

1044 """Return a vector of parameter values. 1045 1046 @keyword model_info: The model index from model_info(). This is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1047 @type model_info: int 1048 @keyword sim_index: The Monte Carlo simulation index. 1049 @type sim_index: int 1050 @return: The vector of parameter values. 1051 @rtype: list of str 1052 """ 1053 1054 # Get the spin containers. 1055 spins = [] 1056 for spin_id in model_info: 1057 # Get the spin. 1058 spin = return_spin(spin_id) 1059 1060 # Skip deselected spins. 1061 if not spin.select: 1062 continue 1063 1064 # Add the spin. 1065 spins.append(spin) 1066 1067 # No spins. 1068 if not len(spins): 1069 return None 1070 1071 # Loop over the parameters of the cluster, fetching their values. 1072 values = [] 1073 for param_name, param_index, si, r20_key in loop_parameters(spins=spins): 1074 values.append(get_value(spins=spins, sim_index=sim_index, param_name=param_name, spin_index=si, r20_key=r20_key)) 1075 1076 # Return all values. 1077 return values

1078 1079

1080 - def grid_search(self, lower=None, upper=None, inc=None, constraints=True, verbosity=1, sim_index=None):

1081 """The relaxation dispersion curve fitting grid search function. 1082 1083 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 1084 @type lower: array of numbers 1085 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 1086 @type upper: array of numbers 1087 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. 1088 @type inc: array of int 1089 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 1090 @type constraints: bool 1091 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 1092 @type verbosity: int 1093 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 1094 @type sim_index: int 1095 """ 1096 1097 # Minimisation. 1098 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

1099 1100

1101 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

1102 """Relaxation dispersion curve fitting function. 1103 1104 @keyword min_algor: The minimisation algorithm to use. 1105 @type min_algor: str 1106 @keyword min_options: An array of options to be used by the minimisation algorithm. 1107 @type min_options: array of str 1108 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1109 @type func_tol: None or float 1110 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1111 @type grad_tol: None or float 1112 @keyword max_iterations: The maximum number of iterations for the algorithm. 1113 @type max_iterations: int 1114 @keyword constraints: If True, constraints are used during optimisation. 1115 @type constraints: bool 1116 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 1117 @type scaling: bool 1118 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 1119 @type verbosity: int 1120 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 1121 @type sim_index: None or int 1122 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 1123 @type lower: array of numbers 1124 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 1125 @type upper: array of numbers 1126 @keyword inc: The increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search. 1127 @type inc: array of int 1128 """ 1129 1130 # Data checks. 1131 check_mol_res_spin_data() 1132 check_model_type() 1133 1134 # Initialise some empty data pipe structures so that the target function set up does not fail. 1135 if not hasattr(cdp, 'cpmg_frqs_list'): 1136 cdp.cpmg_frqs_list = [] 1137 if not hasattr(cdp, 'spin_lock_nu1_list'): 1138 cdp.spin_lock_nu1_list = [] 1139 1140 # Special exponential curve-fitting for the 'R2eff' model. 1141 if cdp.model_type == 'R2eff': 1142 # Sanity checks. 1143 if not has_exponential_exp_type(): 1144 raise RelaxError("The R2eff model with the fixed time period dispersion experiments cannot be optimised.") 1145 1146 # Optimisation. 1147 self._minimise_r2eff(min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, scaling=scaling, verbosity=verbosity, sim_index=sim_index, lower=lower, upper=upper, inc=inc) 1148 1149 # Exit the method. 1150 return 1151 1152 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 1153 processor_box = Processor_box() 1154 processor = processor_box.processor 1155 1156 # The number of time points for the exponential curves (if present). 1157 num_time_pts = 1 1158 if hasattr(cdp, 'num_time_pts'): 1159 num_time_pts = cdp.num_time_pts 1160 1161 # Number of spectrometer fields. 1162 fields = [None] 1163 field_count = 1 1164 if hasattr(cdp, 'spectrometer_frq'): 1165 fields = cdp.spectrometer_frq_list 1166 field_count = cdp.spectrometer_frq_count 1167 1168 # Loop over the spin blocks. 1169 for spin_ids in self.model_loop(): 1170 # The spin containers. 1171 spins = spin_ids_to_containers(spin_ids) 1172 1173 # Skip deselected clusters. 1174 skip = True 1175 for spin in spins: 1176 if spin.select: 1177 skip = False 1178 if skip: 1179 continue 1180 1181 # Diagonal scaling. 1182 scaling_matrix = assemble_scaling_matrix(spins=spins, scaling=scaling) 1183 1184 # Set up the slave command object. 1185 command = Disp_minimise_command(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, verbosity=verbosity, lower=lower, upper=upper, inc=inc, fields=fields, param_names=get_param_names(spins)) 1186 1187 # Set up the memo. 1188 memo = Disp_memo(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix, verbosity=verbosity) 1189 1190 # Add the slave command and memo to the processor queue. 1191 processor.add_to_queue(command, memo)

1192 1193

1194 - def model_desc(self, model_info):

1195 """Return a description of the model. 1196 1197 @param model_info: The model index from model_info(). 1198 @type model_info: int 1199 @return: The model description. 1200 @rtype: str 1201 """ 1202 1203 # The model loop is over the spin clusters, so return a description of the cluster. 1204 return "The spin cluster %s." % model_info

1205 1206

1207 - def model_loop(self):

1208 """Loop over the spin groupings for one model applied to multiple spins. 1209 1210 @return: The list of spins per block will be yielded, as well as the list of spin IDs. 1211 @rtype: tuple of list of SpinContainer instances and list of spin IDs 1212 """ 1213 1214 # The cluster loop. 1215 for spin_ids in loop_cluster(skip_desel=False): 1216 yield spin_ids

1217 1218

1219 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):

1220 """Return the k, n, and chi2 model statistics. 1221 1222 k - number of parameters. 1223 n - number of data points. 1224 chi2 - the chi-squared value. 1225 1226 1227 @keyword model_info: The model index from model_info(). 1228 @type model_info: None or int 1229 @keyword spin_id: The spin identification string. This is ignored in the N-state model. 1230 @type spin_id: None or str 1231 @keyword global_stats: A parameter which determines if global or local statistics are returned. For the N-state model, this argument is ignored. 1232 @type global_stats: None or bool 1233 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2). 1234 @rtype: tuple of (int, int, float) 1235 """ 1236 1237 # Unpack the data. 1238 spin_ids = model_info 1239 spins = spin_ids_to_containers(spin_ids) 1240 1241 # The number of parameters for the cluster. 1242 k = param_num(spins=spins) 1243 1244 # The number of points from all spins. 1245 n = 0 1246 for spin in spins: 1247 # Skip deselected spins. 1248 if not spin.select: 1249 continue 1250 1251 n += len(spin.r2eff) 1252 1253 # Take the chi-squared from the first spin of the cluster (which has a value). 1254 chi2 = None 1255 for spin in spins: 1256 # Skip deselected spins. 1257 if not spin.select: 1258 continue 1259 1260 if hasattr(spin, 'chi2'): 1261 chi2 = spin.chi2 1262 break 1263 1264 # Return the values. 1265 return k, n, chi2

1266 1267

1268 - def overfit_deselect(self, data_check=True, verbose=True):

1269 """Deselect spins which have insufficient data to support minimisation. 1270 1271 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 1272 @type data_check: bool 1273 @keyword verbose: A flag which if True will allow printouts. 1274 @type verbose: bool 1275 """ 1276 1277 # Test the sequence data exists. 1278 check_mol_res_spin_data() 1279 1280 # 1H MMQ flag. 1281 proton_mmq_flag = has_proton_mmq_cpmg() 1282 1283 # Loop over spin data. 1284 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 1285 # Skip protons for MMQ data. 1286 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 1287 continue 1288 1289 # Get the attached proton. 1290 proton = None 1291 if proton_mmq_flag: 1292 # Get all protons. 1293 proton_spins = return_attached_protons(spin_id) 1294 1295 # Only one allowed. 1296 if len(proton_spins) > 1: 1297 print("Multiple protons attached to the spin '%s', but one one attached proton is supported for the MMQ-type models." % spin_id) 1298 spin.select = False 1299 continue 1300 1301 # Alias the single proton. 1302 if len(proton_spins): 1303 proton = proton_spins[0] 1304 1305 # Check if data exists. 1306 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'): 1307 print("No R2eff data could be found, deselecting the '%s' spin." % spin_id) 1308 spin.select = False 1309 continue

1310 1311

1312 - def return_data(self, data_id=None):

1313 """Return the peak intensity data structure. 1314 1315 @param data_id: The spin ID string, as yielded by the base_data_loop() generator method. 1316 @type data_id: str 1317 @return: The peak intensity data structure. 1318 @rtype: list of float 1319 """ 1320 1321 # The R2eff model. 1322 if cdp.model_type == 'R2eff': 1323 # Unpack the data. 1324 spin, exp_type, frq, offset, point = data_id 1325 1326 # Return the data. 1327 return spin.intensities 1328 1329 # All other models. 1330 else: 1331 raise RelaxImplementError

1332 1333 1334 return_data_name_doc = Desc_container("Relaxation dispersion curve fitting data type string matching patterns") 1335 _table = uf_tables.add_table(label="table: dispersion curve-fit data type patterns", caption="Relaxation dispersion curve fitting data type string matching patterns.") 1336 _table.add_headings(["Data type", "Object name"]) 1337 _table.add_row(["Transversal relaxation rate (rad/s)", "'r2'"]) 1338 _table.add_row(["Transversal relaxation rate for state A (rad/s)", "'r2a'"]) 1339 _table.add_row(["Transversal relaxation rate for state B (rad/s)", "'r2b'"]) 1340 _table.add_row(["Population of state A", "'pA'"]) 1341 _table.add_row(["Population of state B", "'pB'"]) 1342 _table.add_row(["Population of state C", "'pC'"]) 1343 _table.add_row(["The pA.pB.dw**2 parameter (ppm^2)", "'phi_ex'"]) 1344 _table.add_row(["The pA.dw**2 parameter (ppm^2)", "'padw2'"]) 1345 _table.add_row(["Chemical shift difference between states A and B (ppm)", "'dw'"]) 1346 _table.add_row(["Chemical shift difference between states A and B for 3-site exchange (ppm)", "'dw_AB'"]) 1347 _table.add_row(["Chemical shift difference between states A and C for 3-site exchange (ppm)", "'dw_AC'"]) 1348 _table.add_row(["Chemical shift difference between states B and C for 3-site exchange (ppm)", "'dw_BC'"]) 1349 _table.add_row(["Proton chemical shift difference between states A and B (ppm)", "'dwH'"]) 1350 _table.add_row(["Proton chemical shift difference between states A and B for 3-site exchange (ppm)", "'dwH_AB'"]) 1351 _table.add_row(["Proton chemical shift difference between states A and C for 3-site exchange (ppm)", "'dwH_AC'"]) 1352 _table.add_row(["Proton chemical shift difference between states B and C for 3-site exchange (ppm)", "'dwH_BC'"]) 1353 _table.add_row(["Exchange rate (rad/s)", "'kex'"]) 1354 _table.add_row(["Exchange rate between sites A and B for 3-site exchange (rad/s)", "'kex_AB'"]) 1355 _table.add_row(["Exchange rate between sites A and C for 3-site exchange (rad/s)", "'kex_AC'"]) 1356 _table.add_row(["Exchange rate between sites B and C for 3-site exchange (rad/s)", "'kex_BC'"]) 1357 _table.add_row(["Exchange rate from state A to state B (rad/s)", "'k_AB'"]) 1358 _table.add_row(["Exchange rate from state B to state A (rad/s)", "'k_BA'"]) 1359 _table.add_row(["Time of exchange (s/rad)", "'tex'"]) 1360 _table.add_row(["Rotating frame tilt angle", "'theta'"]) 1361 _table.add_row(["Effective field in rotating frame", "'w_eff'"]) 1362 _table.add_row(["Peak intensities (series)", "'intensities'"]) 1363 _table.add_row(["CPMG pulse train frequency (series, Hz)", "'cpmg_frqs'"]) 1364 return_data_name_doc.add_table(_table.label) 1365 1366

1367 - def return_error(self, data_id=None):

1368 """Return the standard deviation data structure. 1369 1370 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 1371 @type data_id: SpinContainer instance and float 1372 @return: The standard deviation data structure. 1373 @rtype: list of float 1374 """ 1375 1376 # The R2eff model. 1377 if cdp.model_type == 'R2eff': 1378 # Unpack the data. 1379 spin, exp_type, frq, offset, point = data_id 1380 1381 # Generate the data structure to return. 1382 errors = [] 1383 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 1384 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 1385 1386 # All other models. 1387 else: 1388 # Unpack the data. 1389 spin, spin_id = data_id 1390 1391 # 1H MMQ flag. 1392 proton_mmq_flag = has_proton_mmq_cpmg() 1393 1394 # Get the attached proton. 1395 proton = None 1396 if proton_mmq_flag: 1397 proton = return_attached_protons(spin_id)[0] 1398 1399 # The errors. 1400 r2eff_err = {} 1401 if hasattr(spin, 'r2eff_err'): 1402 r2eff_err.update(spin.r2eff_err) 1403 if hasattr(proton, 'r2eff_err'): 1404 r2eff_err.update(proton.r2eff_err) 1405 return r2eff_err 1406 1407 # Return the error list. 1408 return errors

1409 1410 1411 set_doc = Desc_container("Relaxation dispersion curve fitting set details") 1412 set_doc.add_paragraph("Only three parameters can be set for either the slow- or the fast-exchange regime. For the slow-exchange regime, these parameters include the transversal relaxation rate for state A (R2A), the exchange rate from state A to state (k_AB) and the chemical shift difference between states A and B (dw). For the fast-exchange regime, these include the transversal relaxation rate (R2), the chemical exchange contribution to R2 (Rex) and the exchange rate (kex). Setting parameters for a non selected model has no effect.") 1413 1414

1415 - def return_value(self, spin, param, sim=None, bc=False):

1416 """Return the value and error corresponding to the parameter. 1417 1418 If sim is set to an integer, return the value of the simulation and None. 1419 1420 1421 @param spin: The SpinContainer object. 1422 @type spin: SpinContainer 1423 @param param: The name of the parameter to return values for. 1424 @type param: str 1425 @keyword sim: The Monte Carlo simulation index. 1426 @type sim: None or int 1427 @keyword bc: The back-calculated data flag. If True, then the back-calculated data will be returned rather than the actual data. 1428 @type bc: bool 1429 @return: The value and error corresponding to 1430 @rtype: tuple of length 2 of floats or None 1431 """ 1432 1433 # Define list of special parameters. 1434 special_parameters = ['theta', 'w_eff'] 1435 1436 # Use api_common function for paramets not defined in special_parameters. 1437 if param not in special_parameters: 1438 returnval = self._return_value_general(spin=spin, param=param, sim=sim, bc=bc) 1439 return returnval 1440 1441 # If parameter in special_parameters, the do the following. 1442 1443 # Initial values. 1444 value = None 1445 error = None 1446 1447 # Return the data 1448 theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=spin) 1449 1450 # Return for parameter theta 1451 if param == "theta": 1452 value = theta_spin_dic 1453 1454 # Return for parameter theta 1455 elif param == "w_eff": 1456 value = w_eff_spin_dic 1457 1458 # Return the data. 1459 return value, error

1460 1461

1462 - def set_error(self, model_info, index, error):

1463 """Set the parameter errors. 1464 1465 @param model_info: The spin container originating from model_loop(). 1466 @type model_info: unknown 1467 @param index: The index of the parameter to set the errors for. 1468 @type index: int 1469 @param error: The error value. 1470 @type error: float 1471 """ 1472 1473 # Unpack the data. 1474 spin_ids = model_info 1475 spins = spin_ids_to_containers(spin_ids) 1476 1477 # The number of parameters. 1478 total_param_num = param_num(spins=spins) 1479 1480 # No more model parameters. 1481 model_param = True 1482 if index >= total_param_num: 1483 model_param = False 1484 1485 # The auxiliary cluster parameters. 1486 aux_params = [] 1487 if 'pA' in spins[0].params: 1488 aux_params.append('pB') 1489 if 'pB' in spins[0].params: 1490 aux_params.append('pA') 1491 if 'kex' in spins[0].params: 1492 aux_params.append('tex') 1493 if 'tex' in spins[0].params: 1494 aux_params.append('kex') 1495 1496 # Convert the parameter index. 1497 if model_param: 1498 param_name, si, mi = param_index_to_param_info(index=index, spins=spins) 1499 else: 1500 param_name = aux_params[index - total_param_num] 1501 si = 0 1502 mi = 0 1503 1504 # The parameter error name. 1505 err_name = param_name + "_err" 1506 1507 # The exponential curve parameters. 1508 if param_name in ['r2eff', 'i0']: 1509 # Initialise if needed. 1510 if not hasattr(spins[si], err_name): 1511 setattr(spins[si], err_name, {}) 1512 1513 # Set the value. 1514 setattr(spins[si], err_name, error) 1515 1516 # Model and auxiliary parameters. 1517 else: 1518 for spin in spins: 1519 setattr(spin, err_name, error)

1520 1521

1522 - def set_selected_sim(self, model_info, select_sim):

1523 """Set the simulation selection flag. 1524 1525 @param model_info: The list of spins and spin IDs per cluster originating from model_loop(). 1526 @type model_info: tuple of list of SpinContainer instances and list of spin IDs 1527 @param select_sim: The selection flag for the simulations. 1528 @type select_sim: bool 1529 """ 1530 1531 # Unpack the data. 1532 spin_ids = model_info 1533 spins = spin_ids_to_containers(spin_ids) 1534 1535 # Loop over the spins, storing the structure for each spin. 1536 for spin in spins: 1537 spin.select_sim = deepcopy(select_sim)

1538 1539

1540 - def sim_init_values(self):

1541 """Initialise the Monte Carlo parameter values.""" 1542 1543 # Get the parameter object names. 1544 param_names = self.data_names(set='params') 1545 1546 # Add the names of kex-tex pair. 1547 pairs = {} 1548 pairs['kex'] = 'tex' 1549 pairs['tex'] = 'kex' 1550 1551 # Add the names of pA-pB pair. 1552 pairs['pA'] = 'pB' 1553 pairs['pB'] = 'pA' 1554 1555 # Add the names of kex-k_AB pair and kex-k_BA pair. 1556 pairs['k_AB'] = 'kex' 1557 pairs['k_BA'] = 'kex' 1558 1559 # Get the minimisation statistic object names. 1560 min_names = self.data_names(set='min') 1561 1562 1563 # Test if Monte Carlo parameter values have already been set. 1564 ############################################################# 1565 1566 # Loop over the spins. 1567 for spin in spin_loop(): 1568 # Skip deselected spins. 1569 if not spin.select: 1570 continue 1571 1572 # Loop over all the parameter names. 1573 for object_name in param_names: 1574 # Name for the simulation object. 1575 sim_object_name = object_name + '_sim' 1576 1577 1578 # Set the Monte Carlo parameter values. 1579 ####################################### 1580 1581 # Loop over the spins. 1582 for spin in spin_loop(): 1583 # Skip deselected spins. 1584 if not spin.select: 1585 continue 1586 1587 # Skip protons for MMQ data. 1588 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 1589 continue 1590 1591 # Loop over all the data names. 1592 for object_name in param_names: 1593 # Not a parameter of the model. 1594 if not (object_name in spin.params or (object_name in pairs and pairs[object_name] in spin.params)): 1595 continue 1596 1597 # Name for the simulation object. 1598 sim_object_name = object_name + '_sim' 1599 1600 # Create the simulation object. 1601 setattr(spin, sim_object_name, []) 1602 1603 # Get the simulation object. 1604 sim_object = getattr(spin, sim_object_name) 1605 1606 # Loop over the simulations. 1607 for j in range(cdp.sim_number): 1608 # The non-simulation value. 1609 if object_name not in spin.params: 1610 value = deepcopy(getattr(spin, pairs[object_name])) 1611 else: 1612 value = deepcopy(getattr(spin, object_name)) 1613 1614 # Copy and append the data. 1615 sim_object.append(value) 1616 1617 # Loop over all the minimisation object names. 1618 for object_name in min_names: 1619 # Name for the simulation object. 1620 sim_object_name = object_name + '_sim' 1621 1622 # Create the simulation object. 1623 setattr(spin, sim_object_name, []) 1624 1625 # Get the simulation object. 1626 sim_object = getattr(spin, sim_object_name) 1627 1628 # Loop over the simulations. 1629 for j in range(cdp.sim_number): 1630 # Copy and append the data. 1631 sim_object.append(deepcopy(getattr(spin, object_name)))

1632 1633

1634 - def sim_pack_data(self, data_id, sim_data):

1635 """Pack the Monte Carlo simulation data. 1636 1637 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 1638 @type data_id: SpinContainer instance and float 1639 @param sim_data: The Monte Carlo simulation data. 1640 @type sim_data: list of float 1641 """ 1642 1643 # The R2eff model (with peak intensity base data). 1644 if cdp.model_type == 'R2eff': 1645 # Unpack the data. 1646 spin, exp_type, frq, offset, point = data_id 1647 1648 # Initialise the data structure if needed. 1649 if not hasattr(spin, 'intensity_sim'): 1650 spin.intensity_sim = {} 1651 1652 # Loop over each time point. 1653 ti = 0 1654 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 1655 # Get the intensity keys. 1656 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 1657 1658 # Loop over the intensity keys. 1659 for int_key in int_keys: 1660 # Test if the simulation data point already exists. 1661 if int_key in spin.intensity_sim: 1662 raise RelaxError("Monte Carlo simulation data for the key '%s' already exists." % int_key) 1663 1664 # Initialise the list. 1665 spin.intensity_sim[int_key] = [] 1666 1667 # Loop over the simulations, appending the corresponding data. 1668 for i in range(cdp.sim_number): 1669 spin.intensity_sim[int_key].append(sim_data[i][ti]) 1670 1671 # Increment the time index. 1672 ti += 1 1673 1674 # All other models (with R2eff/R1rho base data). 1675 else: 1676 # Unpack the data. 1677 spin, spin_id = data_id 1678 1679 # 1H MMQ flag. 1680 proton_mmq_flag = has_proton_mmq_cpmg() 1681 1682 # Get the attached proton. 1683 proton = None 1684 if proton_mmq_flag: 1685 proton = return_attached_protons(spin_id)[0] 1686 1687 # Pack the data. 1688 spin.r2eff_sim = sim_data 1689 if proton != None: 1690 proton.r2eff_sim = sim_data

1691 1692

1693 - def sim_return_param(self, model_info, index):

1694 """Return the array of simulation parameter values. 1695 1696 @param model_info: The model information originating from model_loop(). 1697 @type model_info: unknown 1698 @param index: The index of the parameter to return the array of values for. 1699 @type index: int 1700 @return: The array of simulation parameter values. 1701 @rtype: list of float 1702 """ 1703 1704 # Unpack the data. 1705 spin_ids = model_info 1706 spins = spin_ids_to_containers(spin_ids) 1707 1708 # The number of parameters. 1709 total_param_num = param_num(spins=spins) 1710 1711 # No more model parameters. 1712 model_param = True 1713 if index >= total_param_num: 1714 model_param = False 1715 1716 # The auxiliary cluster parameters. 1717 aux_params = [] 1718 for spin in spins: 1719 if not spin.select: 1720 continue 1721 if 'pA' in spin.params: 1722 aux_params.append('pB') 1723 if 'pB' in spin.params: 1724 aux_params.append('pA') 1725 if 'kex' in spin.params: 1726 aux_params.append('tex') 1727 if 'tex' in spin.params: 1728 aux_params.append('kex') 1729 break 1730 1731 # No more auxiliary parameters. 1732 total_aux_num = total_param_num + len(aux_params) 1733 if index >= total_aux_num: 1734 return 1735 1736 # Convert the parameter index. 1737 if model_param: 1738 param_name, si, mi = param_index_to_param_info(index=index, spins=spins) 1739 if not param_name in ['r2eff', 'i0']: 1740 si = 0 1741 else: 1742 param_name = aux_params[index - total_param_num] 1743 si = 0 1744 mi = 0 1745 1746 # The exponential curve parameters. 1747 sim_data = [] 1748 if param_name == 'r2eff': 1749 for i in range(cdp.sim_number): 1750 sim_data.append(spins[si].r2eff_sim[i]) 1751 elif param_name == 'i0': 1752 for i in range(cdp.sim_number): 1753 sim_data.append(spins[si].i0_sim[i]) 1754 1755 # Model and auxiliary parameters. 1756 else: 1757 sim_data = getattr(spins[si], param_name + "_sim") 1758 1759 # Set the sim data to None if empty. 1760 if sim_data == []: 1761 sim_data = None 1762 1763 # Return the object. 1764 return sim_data

1765 1766

1767 - def sim_return_selected(self, model_info):

1768 """Return the array of selected simulation flags. 1769 1770 @param model_info: The list of spins and spin IDs per cluster originating from model_loop(). 1771 @type model_info: tuple of list of SpinContainer instances and list of spin IDs 1772 @return: The array of selected simulation flags. 1773 @rtype: list of int 1774 """ 1775 1776 # Unpack the data. 1777 spin_ids = model_info 1778 spins = spin_ids_to_containers(spin_ids) 1779 1780 # Return the array from the first spin, as this array will be identical for all spins in the cluster. 1781 return spins[0].select_sim

Source Code for Module specific_analyses.relax_disp.api