specific_analyses.relax_disp.optimisation

1 ############################################################################### 2 # # 3 # Copyright (C) 2013-2014 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for the optimisation of the relaxation dispersion models.""" 24 25 # Python module imports. 26 from minfx.generic import generic_minimise 27 from minfx.grid import grid 28 from numpy import dot, float64, int32, ones, zeros 29 from numpy.linalg import inv 30 from re import match, search 31 import sys 32 33 # relax module imports. 34 from dep_check import C_module_exp_fn 35 from lib.check_types import is_float 36 from lib.dispersion.two_point import calc_two_point_r2eff, calc_two_point_r2eff_err 37 from lib.errors import RelaxError 38 from lib.text.sectioning import subsection 39 from multi import Memo, Result_command, Slave_command 40 from pipe_control.mol_res_spin import spin_loop 41 from specific_analyses.relax_disp.checks import check_disp_points, check_exp_type, check_exp_type_fixed_time 42 from specific_analyses.relax_disp.data import average_intensity, count_spins, find_intensity_keys, has_exponential_exp_type, has_proton_mmq_cpmg, loop_exp, loop_exp_frq_offset_point, loop_exp_frq_offset_point_time, loop_frq, loop_offset, loop_time, pack_back_calc_r2eff, return_cpmg_frqs, return_offset_data, return_param_key_from_data, return_r1_data, return_r2eff_arrays, return_spin_lock_nu1 43 from specific_analyses.relax_disp.parameters import assemble_param_vector, assemble_scaling_matrix, disassemble_param_vector, linear_constraints, loop_parameters, param_conversion, param_num 44 from specific_analyses.relax_disp.variables import EXP_TYPE_LIST_CPMG, MODEL_CR72, MODEL_CR72_FULL, MODEL_LIST_MMQ, MODEL_LM63, MODEL_M61, MODEL_M61B, MODEL_MP05, MODEL_TAP03, MODEL_TP02 45 from target_functions.relax_disp import Dispersion 46 47 # C modules. 48 if C_module_exp_fn: 49 from target_functions.relax_fit import setup, func, dfunc, d2func, back_calc_I 50 51

52 -def back_calc_peak_intensities(spin=None, exp_type=None, frq=None, offset=None, point=None):

53 """Back-calculation of peak intensity for the given relaxation time. 54 55 @keyword spin: The specific spin data container. 56 @type spin: SpinContainer instance 57 @keyword exp_type: The experiment type. 58 @type exp_type: str 59 @keyword frq: The spectrometer frequency. 60 @type frq: float 61 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm. 62 @type offset: None or float 63 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 64 @type point: float 65 @return: The back-calculated peak intensities for the given exponential curve. 66 @rtype: numpy rank-1 float array 67 """ 68 69 # Check. 70 if not has_exponential_exp_type(): 71 raise RelaxError("Back-calculation is not allowed for the fixed time experiment types.") 72 73 # The key. 74 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 75 76 # Create the initial parameter vector. 77 param_vector = assemble_param_vector(spins=[spin], key=param_key) 78 79 # Create a scaling matrix. 80 scaling_matrix = assemble_scaling_matrix(spins=[spin], key=param_key, scaling=False) 81 82 # The peak intensities and times. 83 values = [] 84 errors = [] 85 times = [] 86 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 87 # The data. 88 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) 89 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 90 times.append(time) 91 92 # The scaling matrix in a diagonalised list form. 93 scaling_list = [] 94 for i in range(len(scaling_matrix)): 95 scaling_list.append(scaling_matrix[i, i]) 96 97 # Initialise the relaxation fit functions. 98 setup(num_params=len(param_vector), num_times=len(times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list) 99 100 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 101 func(param_vector) 102 103 # Get the data back. 104 results = back_calc_I() 105 106 # Return the correct peak height. 107 return results

108 109

110 -def back_calc_r2eff(spin=None, spin_id=None, cpmg_frqs=None, spin_lock_nu1=None, store_chi2=False):

111 """Back-calculation of R2eff/R1rho values for the given spin. 112 113 @keyword spin: The specific spin data container. 114 @type spin: SpinContainer instance 115 @keyword spin_id: The spin ID string for the spin container. 116 @type spin_id: str 117 @keyword cpmg_frqs: The CPMG frequencies to use instead of the user loaded values - to enable interpolation. 118 @type cpmg_frqs: list of lists of numpy rank-1 float arrays 119 @keyword spin_lock_nu1: The spin-lock field strengths to use instead of the user loaded values - to enable interpolation. 120 @type spin_lock_nu1: list of lists of numpy rank-1 float arrays 121 @keyword store_chi2: A flag which if True will cause the spin specific chi-squared value to be stored in the spin container. 122 @type store_chi2: bool 123 @return: The back-calculated R2eff/R1rho value for the given spin. 124 @rtype: numpy rank-1 float array 125 """ 126 127 # Skip protons for MMQ data. 128 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 129 return 130 131 # Create the initial parameter vector. 132 param_vector = assemble_param_vector(spins=[spin]) 133 134 # Create a scaling matrix. 135 scaling_matrix = assemble_scaling_matrix(spins=[spin], scaling=False) 136 137 # Number of spectrometer fields. 138 fields = [None] 139 field_count = 1 140 if hasattr(cdp, 'spectrometer_frq_count'): 141 fields = cdp.spectrometer_frq_list 142 field_count = cdp.spectrometer_frq_count 143 144 # Initialise the data structures for the target function. 145 values, errors, missing, frqs, frqs_H, exp_types, relax_times = return_r2eff_arrays(spins=[spin], spin_ids=[spin_id], fields=fields, field_count=field_count) 146 147 # The offset and R1 data. 148 chemical_shifts, offsets, tilt_angles, Delta_omega, w_eff = return_offset_data(spins=[spin], spin_ids=[spin_id], field_count=field_count, fields=spin_lock_nu1) 149 r1 = return_r1_data(spins=[spin], spin_ids=[spin_id], field_count=field_count) 150 151 # The dispersion data. 152 recalc_tau = True 153 if cpmg_frqs == None and spin_lock_nu1 == None: 154 cpmg_frqs = return_cpmg_frqs(ref_flag=False) 155 spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False) 156 157 # Reconstruct the structures for interpolation. 158 else: 159 recalc_tau = False 160 values = [] 161 errors = [] 162 missing = [] 163 for exp_type, ei in loop_exp(return_indices=True): 164 values.append([]) 165 errors.append([]) 166 missing.append([]) 167 for si in range(1): 168 values[ei].append([]) 169 errors[ei].append([]) 170 missing[ei].append([]) 171 for frq, mi in loop_frq(return_indices=True): 172 values[ei][si].append([]) 173 errors[ei][si].append([]) 174 missing[ei][si].append([]) 175 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 176 if exp_type in EXP_TYPE_LIST_CPMG: 177 num = len(cpmg_frqs[ei][mi][oi]) 178 else: 179 num = len(spin_lock_nu1[ei][mi][oi]) 180 values[ei][si][mi].append(zeros(num, float64)) 181 errors[ei][si][mi].append(ones(num, float64)) 182 missing[ei][si][mi].append(zeros(num, int32)) 183 184 # Initialise the relaxation dispersion fit functions. 185 model = Dispersion(model=spin.model, num_params=param_num(spins=[spin]), num_spins=1, num_frq=field_count, exp_types=exp_types, values=values, errors=errors, missing=missing, frqs=frqs, frqs_H=frqs_H, cpmg_frqs=cpmg_frqs, spin_lock_nu1=spin_lock_nu1, chemical_shifts=chemical_shifts, offset=offsets, tilt_angles=tilt_angles, r1=r1, relax_times=relax_times, scaling_matrix=scaling_matrix, recalc_tau=recalc_tau) 186 187 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 188 chi2 = model.func(param_vector) 189 190 # Store the chi-squared value. 191 if store_chi2: 192 spin.chi2 = chi2 193 194 # Return the structure. 195 return model.back_calc

196 197

198 -def calculate_r2eff():

199 """Calculate the R2eff values for fixed relaxation time period data.""" 200 201 # Data checks. 202 check_exp_type() 203 check_disp_points() 204 check_exp_type_fixed_time() 205 206 # Printouts. 207 print("Calculating the R2eff/R1rho values for fixed relaxation time period data.") 208 209 # Loop over the spins. 210 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 211 # Spin ID printout. 212 print("Spin '%s'." % spin_id) 213 214 # Skip spins which have no data. 215 if not hasattr(spin, 'peak_intensity'): 216 continue 217 218 # Initialise the data structures. 219 if not hasattr(spin, 'r2eff'): 220 spin.r2eff = {} 221 if not hasattr(spin, 'r2eff_err'): 222 spin.r2eff_err = {} 223 224 # Loop over all the data. 225 for exp_type, frq, offset, point, time in loop_exp_frq_offset_point_time(): 226 # The three keys. 227 ref_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=None, time=time) 228 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 229 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 230 231 # Check for missing data. 232 missing = False 233 for i in range(len(ref_keys)): 234 if ref_keys[i] not in spin.peak_intensity: 235 missing = True 236 for i in range(len(int_keys)): 237 if int_keys[i] not in spin.peak_intensity: 238 missing = True 239 if missing: 240 continue 241 242 # Average the reference intensity data and errors. 243 ref_intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time) 244 ref_intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time, error=True) 245 246 # Average the intensity data and errors. 247 intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 248 intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True) 249 250 # Calculate the R2eff value. 251 spin.r2eff[param_key] = calc_two_point_r2eff(relax_time=time, I_ref=ref_intensity, I=intensity) 252 253 # Calculate the R2eff error. 254 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)

255 256

257 -def grid_search_setup(spins=None, spin_ids=None, param_vector=None, lower=None, upper=None, inc=None, scaling_matrix=None):

258 """The grid search setup function. 259 260 @keyword spins: The list of spin data containers for the block. 261 @type spins: list of SpinContainer instances 262 @keyword spin_ids: The corresponding spin ID strings. 263 @type spin_ids: list of str 264 @keyword param_vector: The parameter vector. 265 @type param_vector: numpy array 266 @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. 267 @type lower: array of numbers 268 @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. 269 @type upper: array of numbers 270 @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. 271 @type inc: array of int 272 @keyword scaling_matrix: The scaling matrix. 273 @type scaling_matrix: numpy diagonal matrix 274 @return: A tuple of the grid size and the minimisation options. For the minimisation options, the first dimension corresponds to the model parameter. The second dimension is a list of the number of increments, the lower bound, and upper bound. 275 @rtype: (int, list of lists [int, float, float]) 276 """ 277 278 # The length of the parameter array. 279 n = len(param_vector) 280 281 # Make sure that the length of the parameter array is > 0. 282 if n == 0: 283 raise RelaxError("Cannot run a grid search on a model with zero parameters.") 284 285 # Lower bounds. 286 if lower != None and len(lower) != n: 287 raise RelaxLenError('lower bounds', n) 288 289 # Upper bounds. 290 if upper != None and len(upper) != n: 291 raise RelaxLenError('upper bounds', n) 292 293 # Increment. 294 if isinstance(inc, list) and len(inc) != n: 295 raise RelaxLenError('increment', n) 296 elif isinstance(inc, int): 297 inc = [inc]*n 298 299 # Set up the default bounds. 300 if not lower: 301 # Init. 302 lower = [] 303 upper = [] 304 305 # The R2eff model. 306 if cdp.model_type == 'R2eff': 307 # Loop over each experiment type, spectrometer frequency, offset and dispersion point. 308 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 309 # Loop over the parameters. 310 for param_name, param_index, si, r20_key in loop_parameters(spins=spins): 311 # R2eff relaxation rate (from 1 to 40 s^-1). 312 if param_name == 'r2eff': 313 lower.append(1.0) 314 upper.append(40.0) 315 316 # Intensity. 317 elif param_name == 'i0': 318 lower.append(0.0001) 319 upper.append(max(spins[si].peak_intensity.values())) 320 321 # All other models. 322 else: 323 # Loop over the parameters. 324 for param_name, param_index, si, r20_key in loop_parameters(spins=spins): 325 # Cluster specific parameter. 326 if si == None: 327 si = 0 328 329 # R2 relaxation rates (from 5 to 20 s^-1). 330 if param_name in ['r2', 'r2a', 'r2b']: 331 lower.append(5.0) 332 upper.append(20.0) 333 334 # The pA.pB.dw**2 and pA.dw**2 parameters. 335 elif param_name in ['phi_ex', 'phi_ex_B', 'phi_ex_C', 'padw2']: 336 lower.append(0.0) 337 upper.append(10.0) 338 339 # Chemical shift difference between states A and B (heteronucleus). 340 elif param_name in ['dw', 'dw_AB', 'dw_AC', 'dw_BC']: 341 if spins[si].model in MODEL_LIST_MMQ: 342 lower.append(-10.0) 343 else: 344 lower.append(0.0) 345 upper.append(10.0) 346 347 # Chemical shift difference between states A and B (proton). 348 elif param_name in ['dwH', 'dwH_AB', 'dwH_AC', 'dwH_BC']: 349 if spins[si].model in MODEL_LIST_MMQ: 350 lower.append(-3.0) 351 else: 352 lower.append(0.0) 353 upper.append(3.0) 354 355 # The population of state A. 356 elif param_name == 'pA': 357 if spins[si].model == MODEL_M61B: 358 lower.append(0.85) 359 else: 360 lower.append(0.5) 361 upper.append(1.0) 362 363 # The population of state B (for 3-site exchange). 364 elif param_name == 'pB': 365 lower.append(0.0) 366 upper.append(0.5) 367 368 # Exchange rates. 369 elif param_name in ['kex', 'kex_AB', 'kex_AC', 'kex_BC', 'k_AB', 'kB', 'kC']: 370 lower.append(1.0) 371 upper.append(10000.0) 372 373 # Time of exchange. 374 elif param_name in ['tex']: 375 lower.append(1/10000.0) 376 upper.append(1.0) 377 378 # Pre-set parameters. 379 for param_name, param_index, si, r20_key in loop_parameters(spins=spins): 380 # Cluster specific parameter. 381 if si == None: 382 si = 0 383 384 # Get the parameter. 385 if hasattr(spins[si], param_name): 386 val = getattr(spins[si], param_name) 387 388 # Value already set. 389 if is_float(val) and val != 0.0: 390 # Printout. 391 print("The spin '%s' parameter '%s' is pre-set to %s, skipping it in the grid search." % (spin_ids[si], param_name, val)) 392 393 # Turn of the grid search for this parameter. 394 inc[param_index] = 1 395 lower[param_index] = val 396 upper[param_index] = val 397 398 # Test if the value is a dict, for example for r2. 399 if isinstance(val, dict) and len(val) > 0: 400 # Test if r20_key exists. 401 if r20_key != None: 402 try: 403 val_dic = val[r20_key] 404 except KeyError: 405 print("The key:%s does not exist"%r20_key) 406 continue 407 408 if is_float(val_dic): 409 # Printout. 410 print("The spin '%s' parameter %s '%s[%i]' is pre-set to %s, skipping it in the grid search." % (spin_ids[si], r20_key, param_name, param_index, val_dic)) 411 412 # Turn of the grid search for this parameter. 413 inc[param_index] = 1 414 lower[param_index] = val_dic 415 upper[param_index] = val_dic 416 417 # The full grid size. 418 grid_size = 1 419 for i in range(n): 420 grid_size *= inc[i] 421 422 # Diagonal scaling of minimisation options. 423 lower_new = [] 424 upper_new = [] 425 for i in range(n): 426 lower_new.append(lower[i] / scaling_matrix[i, i]) 427 upper_new.append(upper[i] / scaling_matrix[i, i]) 428 429 # Return the data structures. 430 return grid_size, inc, lower_new, upper_new

431 432

433 -def minimise_r2eff(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):

434 """Optimise the R2eff model by fitting the 2-parameter exponential curves. 435 436 This mimics the R1 and R2 relax_fit analysis. 437 438 439 @keyword min_algor: The minimisation algorithm to use. 440 @type min_algor: str 441 @keyword min_options: An array of options to be used by the minimisation algorithm. 442 @type min_options: array of str 443 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 444 @type func_tol: None or float 445 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 446 @type grad_tol: None or float 447 @keyword max_iterations: The maximum number of iterations for the algorithm. 448 @type max_iterations: int 449 @keyword constraints: If True, constraints are used during optimisation. 450 @type constraints: bool 451 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 452 @type scaling: bool 453 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 454 @type verbosity: int 455 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 456 @type sim_index: None or int 457 @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. 458 @type lower: array of numbers 459 @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. 460 @type upper: array of numbers 461 @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. 462 @type inc: array of int 463 """ 464 465 # Check that the C modules have been compiled. 466 if not C_module_exp_fn: 467 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 468 469 # Loop over the spins. 470 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 471 # Skip spins which have no data. 472 if not hasattr(spin, 'peak_intensity'): 473 continue 474 475 # Loop over each spectrometer frequency and dispersion point. 476 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 477 # The parameter key. 478 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 479 480 # The initial parameter vector. 481 param_vector = assemble_param_vector(spins=[spin], key=param_key, sim_index=sim_index) 482 483 # Diagonal scaling. 484 scaling_matrix = assemble_scaling_matrix(spins=[spin], key=param_key, scaling=scaling) 485 if len(scaling_matrix): 486 param_vector = dot(inv(scaling_matrix), param_vector) 487 488 # Get the grid search minimisation options. 489 lower_new, upper_new = None, None 490 if match('^[Gg]rid', min_algor): 491 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) 492 493 # Linear constraints. 494 A, b = None, None 495 if constraints: 496 A, b = linear_constraints(spins=[spin], scaling_matrix=scaling_matrix) 497 498 # Print out. 499 if verbosity >= 1: 500 # Individual spin section. 501 top = 2 502 if verbosity >= 2: 503 top += 2 504 text = "Fitting to spin %s, frequency %s and dispersion point %s" % (spin_id, frq, point) 505 subsection(file=sys.stdout, text=text, prespace=top) 506 507 # Grid search printout. 508 if match('^[Gg]rid', min_algor): 509 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % grid_size) 510 511 # The peak intensities, errors and times. 512 values = [] 513 errors = [] 514 times = [] 515 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 516 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, sim_index=sim_index)) 517 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 518 times.append(time) 519 520 # The scaling matrix in a diagonalised list form. 521 scaling_list = [] 522 for i in range(len(scaling_matrix)): 523 scaling_list.append(scaling_matrix[i, i]) 524 525 # Initialise the function to minimise. 526 setup(num_params=len(param_vector), num_times=len(times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_list) 527 528 # Grid search. 529 if search('^[Gg]rid', min_algor): 530 results = grid(func=func, args=(), num_incs=inc_new, lower=lower_new, upper=upper_new, A=A, b=b, verbosity=verbosity) 531 532 # Unpack the results. 533 param_vector, chi2, iter_count, warning = results 534 f_count = iter_count 535 g_count = 0.0 536 h_count = 0.0 537 538 # Minimisation. 539 else: 540 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) 541 542 # Unpack the results. 543 if results == None: 544 return 545 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 546 547 # Scaling. 548 if scaling: 549 param_vector = dot(scaling_matrix, param_vector) 550 551 # Disassemble the parameter vector. 552 disassemble_param_vector(param_vector=param_vector, spins=[spin], key=param_key, sim_index=sim_index) 553 554 # Monte Carlo minimisation statistics. 555 if sim_index != None: 556 # Chi-squared statistic. 557 spin.chi2_sim[sim_index] = chi2 558 559 # Iterations. 560 spin.iter_sim[sim_index] = iter_count 561 562 # Function evaluations. 563 spin.f_count_sim[sim_index] = f_count 564 565 # Gradient evaluations. 566 spin.g_count_sim[sim_index] = g_count 567 568 # Hessian evaluations. 569 spin.h_count_sim[sim_index] = h_count 570 571 # Warning. 572 spin.warning_sim[sim_index] = warning 573 574 # Normal statistics. 575 else: 576 # Chi-squared statistic. 577 spin.chi2 = chi2 578 579 # Iterations. 580 spin.iter = iter_count 581 582 # Function evaluations. 583 spin.f_count = f_count 584 585 # Gradient evaluations. 586 spin.g_count = g_count 587 588 # Hessian evaluations. 589 spin.h_count = h_count 590 591 # Warning. 592 spin.warning = warning

593 594 595 596

597 -class Disp_memo(Memo):

598 """The relaxation dispersion memo class.""" 599

600 - def __init__(self, spins=None, spin_ids=None, sim_index=None, scaling_matrix=None, verbosity=None):

601 """Initialise the relaxation dispersion memo class. 602 603 This is used for handling the optimisation results returned from a slave processor. It runs on the master processor and is used to store data which is passed to the slave processor and then passed back to the master via the results command. 604 605 606 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored. 607 @type spins: list of SpinContainer instances 608 @keyword spin_ids: The spin ID strings for the cluster. 609 @type spin_ids: list of str 610 @keyword sim_index: The optional MC simulation index. 611 @type sim_index: int 612 @keyword scaling_matrix: The diagonal, square scaling matrix. 613 @type scaling_matrix: numpy diagonal matrix 614 @keyword verbosity: The verbosity level. This is used by the result command returned to the master for printouts. 615 @type verbosity: int 616 """ 617 618 # Execute the base class __init__() method. 619 super(Disp_memo, self).__init__() 620 621 # Store the arguments. 622 self.spins = spins 623 self.spin_ids = spin_ids 624 self.sim_index = sim_index 625 self.scaling_matrix = scaling_matrix 626 self.verbosity = verbosity

627 628 629

630 -class Disp_minimise_command(Slave_command):

631 """Command class for relaxation dispersion optimisation on the slave processor.""" 632

633 - def __init__(self, spins=None, spin_ids=None, sim_index=None, scaling_matrix=None, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, verbosity=0, lower=None, upper=None, inc=None, fields=None, param_names=None):

634 """Initialise the base class, storing all the master data to be sent to the slave processor. 635 636 This method is run on the master processor whereas the run() method is run on the slave processor. 637 638 639 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored. 640 @type spins: list of SpinContainer instances 641 @keyword spin_ids: The list of spin ID strings corresponding to the spins argument. 642 @type spin_ids: list of str 643 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 644 @type sim_index: None or int 645 @keyword scaling_matrix: The diagonal, square scaling matrix. 646 @type scaling_matrix: numpy diagonal matrix 647 @keyword min_algor: The minimisation algorithm to use. 648 @type min_algor: str 649 @keyword min_options: An array of options to be used by the minimisation algorithm. 650 @type min_options: array of str 651 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 652 @type func_tol: None or float 653 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 654 @type grad_tol: None or float 655 @keyword max_iterations: The maximum number of iterations for the algorithm. 656 @type max_iterations: int 657 @keyword constraints: If True, constraints are used during optimisation. 658 @type constraints: bool 659 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 660 @type verbosity: int 661 @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. 662 @type lower: array of numbers 663 @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. 664 @type upper: array of numbers 665 @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. 666 @type inc: array of int 667 @keyword fields: The list of unique of spectrometer field strengths. 668 @type fields: int 669 @keyword param_names: The list of parameter names to use in printouts. 670 @type param_names: str 671 """ 672 673 # Execute the base class __init__() method. 674 super(Disp_minimise_command, self).__init__() 675 676 # Store the arguments needed by the run() method. 677 self.spins = spins 678 self.spin_ids = spin_ids 679 self.sim_index = sim_index 680 self.scaling_matrix = scaling_matrix 681 self.verbosity = verbosity 682 self.min_algor = min_algor 683 self.min_options = min_options 684 self.func_tol = func_tol 685 self.grad_tol = grad_tol 686 self.max_iterations = max_iterations 687 self.fields = fields 688 self.param_names = param_names 689 690 # Create the initial parameter vector. 691 self.param_vector = assemble_param_vector(spins=self.spins) 692 if len(scaling_matrix): 693 self.param_vector = dot(inv(scaling_matrix), self.param_vector) 694 695 # Get the grid search minimisation options. 696 self.lower_new, self.upper_new = None, None 697 if search('^[Gg]rid', min_algor): 698 self.grid_size, self.inc_new, self.lower_new, self.upper_new = grid_search_setup(spins=spins, spin_ids=spin_ids, param_vector=self.param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=self.scaling_matrix) 699 700 # Linear constraints. 701 self.A, self.b = None, None 702 if constraints: 703 self.A, self.b = linear_constraints(spins=spins, scaling_matrix=scaling_matrix) 704 705 # Test if the spectrometer frequencies have been set. 706 if spins[0].model in [MODEL_LM63, MODEL_CR72, MODEL_CR72_FULL, MODEL_M61, MODEL_TP02, MODEL_TAP03, MODEL_MP05] and not hasattr(cdp, 'spectrometer_frq'): 707 raise RelaxError("The spectrometer frequency information has not been specified.") 708 709 # The R2eff/R1rho data. 710 self.values, self.errors, self.missing, self.frqs, self.frqs_H, self.exp_types, self.relax_times = return_r2eff_arrays(spins=spins, spin_ids=spin_ids, fields=fields, field_count=len(fields), sim_index=sim_index) 711 712 # The offset and R1 data. 713 self.chemical_shifts, self.offsets, self.tilt_angles, self.Delta_omega, self.w_eff = return_offset_data(spins=spins, spin_ids=spin_ids, field_count=len(fields)) 714 self.r1 = return_r1_data(spins=spins, spin_ids=spin_ids, field_count=len(fields), sim_index=sim_index) 715 716 # Parameter number. 717 self.param_num = param_num(spins=spins) 718 719 # The dispersion data. 720 self.dispersion_points = cdp.dispersion_points 721 self.cpmg_frqs = return_cpmg_frqs(ref_flag=False) 722 self.spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False)

723 724

725 - def run(self, processor, completed):

726 """Set up and perform the optimisation.""" 727 728 # Print out. 729 if self.verbosity >= 1: 730 # Individual spin block section. 731 top = 2 732 if self.verbosity >= 2: 733 top += 2 734 subsection(file=sys.stdout, text="Fitting to the spin block %s"%self.spin_ids, prespace=top) 735 736 # Grid search printout. 737 if search('^[Gg]rid', self.min_algor): 738 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % self.grid_size) 739 740 # Initialise the function to minimise. 741 model = Dispersion(model=self.spins[0].model, num_params=self.param_num, num_spins=count_spins(self.spins), num_frq=len(self.fields), exp_types=self.exp_types, values=self.values, errors=self.errors, missing=self.missing, frqs=self.frqs, frqs_H=self.frqs_H, cpmg_frqs=self.cpmg_frqs, spin_lock_nu1=self.spin_lock_nu1, chemical_shifts=self.chemical_shifts, offset=self.offsets, tilt_angles=self.tilt_angles, r1=self.r1, relax_times=self.relax_times, scaling_matrix=self.scaling_matrix) 742 743 # Grid search. 744 if search('^[Gg]rid', self.min_algor): 745 results = grid(func=model.func, args=(), num_incs=self.inc_new, lower=self.lower_new, upper=self.upper_new, A=self.A, b=self.b, verbosity=self.verbosity) 746 747 # Unpack the results. 748 param_vector, chi2, iter_count, warning = results 749 f_count = iter_count 750 g_count = 0.0 751 h_count = 0.0 752 753 # Minimisation. 754 else: 755 results = generic_minimise(func=model.func, args=(), x0=self.param_vector, min_algor=self.min_algor, min_options=self.min_options, func_tol=self.func_tol, grad_tol=self.grad_tol, maxiter=self.max_iterations, A=self.A, b=self.b, full_output=True, print_flag=self.verbosity) 756 757 # Unpack the results. 758 if results == None: 759 return 760 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 761 762 # Optimisation printout. 763 if self.verbosity: 764 print("\nOptimised parameter values:") 765 for i in range(len(param_vector)): 766 print("%-20s %25.15f" % (self.param_names[i], param_vector[i]*self.scaling_matrix[i, i])) 767 768 # Printout. 769 if self.sim_index != None: 770 print("Simulation %s, cluster %s" % (self.sim_index+1, self.spin_ids)) 771 772 # Create the result command object to send back to the master. 773 processor.return_object(Disp_result_command(processor=processor, memo_id=self.memo_id, param_vector=param_vector, chi2=chi2, iter_count=iter_count, f_count=f_count, g_count=g_count, h_count=h_count, warning=warning, missing=self.missing, back_calc=model.back_calc, completed=False))

774 775 776

777 -class Disp_result_command(Result_command):

778 """Class for processing the dispersion optimisation results. 779 780 This object will be sent from the slave back to the master to have its run() method executed. 781 """ 782

783 - def __init__(self, processor=None, memo_id=None, param_vector=None, chi2=None, iter_count=None, f_count=None, g_count=None, h_count=None, warning=None, missing=None, back_calc=None, completed=True):

784 """Set up this class object on the slave, placing the minimisation results here. 785 786 @keyword processor: The processor object. 787 @type processor: multi.processor.Processor instance 788 @keyword memo_id: The memo identification string. 789 @type memo_id: str 790 @keyword param_vector: The optimised parameter vector. 791 @type param_vector: numpy rank-1 array 792 @keyword chi2: The final target function value. 793 @type chi2: float 794 @keyword iter_count: The number of optimisation iterations. 795 @type iter_count: int 796 @keyword f_count: The total function call count. 797 @type f_count: int 798 @keyword g_count: The total gradient call count. 799 @type g_count: int 800 @keyword h_count: The total Hessian call count. 801 @type h_count: int 802 @keyword warning: Any optimisation warnings. 803 @type warning: str or None 804 @keyword missing: The data structure indicating which R2eff/R1rho' base data is missing. 805 @type missing: numpy rank-3 array 806 @keyword back_calc: The back-calculated R2eff/R1rho' data structure from the target function class. This is will be transfered to the master to be stored in the r2eff_bc data structure. 807 @type back_calc: numpy rank-3 array 808 @keyword completed: A flag which if True signals that the optimisation successfully completed. 809 @type completed: bool 810 """ 811 812 # Execute the base class __init__() method. 813 super(Disp_result_command, self).__init__(processor=processor, completed=completed) 814 815 # Store the arguments (to be sent back to the master). 816 self.memo_id = memo_id 817 self.param_vector = param_vector 818 self.chi2 = chi2 819 self.iter_count = iter_count 820 self.f_count = f_count 821 self.g_count = g_count 822 self.h_count = h_count 823 self.warning = warning 824 self.missing = missing 825 self.back_calc = back_calc 826 self.completed = completed

827 828

829 - def run(self, processor=None, memo=None):

830 """Disassemble the model-free optimisation results (on the master). 831 832 @param processor: Unused! 833 @type processor: None 834 @param memo: The dispersion memo. This holds a lot of the data and objects needed for processing the results from the slave. 835 @type memo: memo 836 """ 837 838 # Scaling. 839 if memo.scaling_matrix != None: 840 param_vector = dot(memo.scaling_matrix, self.param_vector) 841 842 # Disassemble the parameter vector. 843 disassemble_param_vector(param_vector=param_vector, spins=memo.spins, sim_index=memo.sim_index) 844 param_conversion(spins=memo.spins, sim_index=memo.sim_index) 845 846 # Monte Carlo minimisation statistics. 847 if memo.sim_index != None: 848 for spin in memo.spins: 849 # Skip deselected spins. 850 if not spin.select: 851 continue 852 853 # Chi-squared statistic. 854 spin.chi2_sim[memo.sim_index] = self.chi2 855 856 # Iterations. 857 spin.iter_sim[memo.sim_index] = self.iter_count 858 859 # Function evaluations. 860 spin.f_count_sim[memo.sim_index] = self.f_count 861 862 # Gradient evaluations. 863 spin.g_count_sim[memo.sim_index] = self.g_count 864 865 # Hessian evaluations. 866 spin.h_count_sim[memo.sim_index] = self.h_count 867 868 # Warning. 869 spin.warning_sim[memo.sim_index] = self.warning 870 871 # Normal statistics. 872 else: 873 for spin in memo.spins: 874 # Skip deselected spins. 875 if not spin.select: 876 continue 877 878 # Chi-squared statistic. 879 spin.chi2 = self.chi2 880 881 # Iterations. 882 spin.iter = self.iter_count 883 884 # Function evaluations. 885 spin.f_count = self.f_count 886 887 # Gradient evaluations. 888 spin.g_count = self.g_count 889 890 # Hessian evaluations. 891 spin.h_count = self.h_count 892 893 # Warning. 894 spin.warning = self.warning 895 896 # Store the back-calculated values. 897 if memo.sim_index == None: 898 # 1H MMQ flag. 899 proton_mmq_flag = has_proton_mmq_cpmg() 900 901 # Loop over each spin, packing the data. 902 si = 0 903 for spin_index in range(len(memo.spins)): 904 # Skip deselected spins. 905 if not memo.spins[spin_index].select: 906 continue 907 908 # Pack the data. 909 pack_back_calc_r2eff(spin=memo.spins[spin_index], spin_id=memo.spin_ids[spin_index], si=si, back_calc=self.back_calc, proton_mmq_flag=proton_mmq_flag) 910 911 # Increment the spin index. 912 si += 1

913

Source Code for Module specific_analyses.relax_disp.optimisation