specific_analyses.relax_disp.optimisation

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

52 -def back_calc_peak_intensities(spin=None, spin_id=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 spin_id: The optional spin ID string for use in warning messages. 58 @type spin_id: str or None 59 @keyword exp_type: The experiment type. 60 @type exp_type: str 61 @keyword frq: The spectrometer frequency. 62 @type frq: float 63 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm. 64 @type offset: None or float 65 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 66 @type point: float 67 @return: The back-calculated peak intensities for the given exponential curve. 68 @rtype: numpy rank-1 float array 69 """ 70 71 # Check. 72 if not has_exponential_exp_type(): 73 raise RelaxError("Back-calculation is not allowed for the fixed time experiment types.") 74 75 # The key. 76 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 77 78 # Create the initial parameter vector. 79 param_vector = assemble_param_vector(spins=[spin], key=param_key) 80 81 # The peak intensities and times. 82 values = [] 83 errors = [] 84 times = [] 85 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 86 # Check the peak intensity keys. 87 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 88 for i in range(len(int_keys)): 89 if int_keys[i] not in spin.peak_intensity: 90 if spin_id: 91 warn(RelaxWarning("The spin %s peak intensity key '%s' is not present, skipping the back-calculation." % (spin_id, int_keys[i]))) 92 else: 93 warn(RelaxWarning("The peak intensity key '%s' is not present, skipping the back-calculation." % int_keys[i])) 94 return 95 96 # The data. 97 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) 98 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 99 times.append(time) 100 101 # The scaling matrix in a diagonalised list form. 102 scaling_list = [] 103 for i in range(len(param_vector)): 104 scaling_list.append(1.0) 105 106 # Initialise the relaxation fit functions. 107 model = Relax_fit_opt(model='exp', num_params=len(spin.params), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 108 109 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 110 model.func(param_vector) 111 112 # Get the data back. 113 results = model.back_calc_data() 114 115 # Return the correct peak height. 116 return results

117 118

119 -def back_calc_r2eff(spins=None, spin_ids=None, cpmg_frqs=None, spin_lock_offset=None, spin_lock_nu1=None, relax_times_new=None, store_chi2=False):

120 """Back-calculation of R2eff/R1rho values for the given spin. 121 122 @keyword spins: The list of specific spin data container for cluster. 123 @type spins: List of SpinContainer instances 124 @keyword spin_ids: The list of spin ID strings for the spin containers in cluster. 125 @type spin_ids: list of str 126 @keyword cpmg_frqs: The CPMG frequencies to use instead of the user loaded values - to enable interpolation. 127 @type cpmg_frqs: list of lists of numpy rank-1 float arrays 128 @keyword spin_lock_offset: The spin-lock offsets to use instead of the user loaded values - to enable interpolation. 129 @type spin_lock_offset: list of lists of numpy rank-1 float arrays 130 @keyword spin_lock_nu1: The spin-lock field strengths to use instead of the user loaded values - to enable interpolation. 131 @type spin_lock_nu1: list of lists of numpy rank-1 float arrays 132 @keyword relax_times_new: The interpolated experiment specific fixed time period for relaxation (in seconds). The dimensions are {Ei, Mi, Oi, Di, Ti}. 133 @type relax_times_new: rank-4 list of floats 134 @keyword store_chi2: A flag which if True will cause the spin specific chi-squared value to be stored in the spin container. 135 @type store_chi2: bool 136 @return: The back-calculated R2eff/R1rho value for the given spin. 137 @rtype: numpy rank-1 float array 138 """ 139 140 # Create the initial parameter vector. 141 param_vector = assemble_param_vector(spins=spins) 142 143 # Number of spectrometer fields. 144 fields = [None] 145 field_count = 1 146 if hasattr(cdp, 'spectrometer_frq_count'): 147 fields = cdp.spectrometer_frq_list 148 field_count = cdp.spectrometer_frq_count 149 150 # Initialise the data structures for the target function. 151 values, errors, missing, frqs, frqs_H, exp_types, relax_times = return_r2eff_arrays(spins=spins, spin_ids=spin_ids, fields=fields, field_count=field_count) 152 153 # The offset and R1 data. 154 r1_setup() 155 offsets, spin_lock_fields_inter, chemical_shifts, tilt_angles, Delta_omega, w_eff = return_offset_data(spins=spins, spin_ids=spin_ids, field_count=field_count, spin_lock_offset=spin_lock_offset, fields=spin_lock_nu1) 156 r1 = return_r1_data(spins=spins, spin_ids=spin_ids, field_count=field_count) 157 r1_fit = is_r1_optimised(spins[0].model) 158 159 # The relaxation times. 160 if relax_times_new != None: 161 relax_times = relax_times_new 162 163 # The dispersion data. 164 recalc_tau = True 165 if cpmg_frqs == None and spin_lock_nu1 == None and spin_lock_offset == None: 166 cpmg_frqs = return_cpmg_frqs(ref_flag=False) 167 spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False) 168 169 # Reset the cpmg_frqs if interpolating R1rho models. 170 elif cpmg_frqs == None and spin_lock_offset != None: 171 cpmg_frqs = None 172 spin_lock_nu1 = spin_lock_fields_inter 173 174 recalc_tau = False 175 values = [] 176 errors = [] 177 missing = [] 178 for exp_type, ei in loop_exp(return_indices=True): 179 values.append([]) 180 errors.append([]) 181 missing.append([]) 182 for si in range(len(spins)): 183 values[ei].append([]) 184 errors[ei].append([]) 185 missing[ei].append([]) 186 for frq, mi in loop_frq(return_indices=True): 187 values[ei][si].append([]) 188 errors[ei][si].append([]) 189 missing[ei][si].append([]) 190 for oi, offset in enumerate(offsets[ei][si][mi]): 191 num = len(spin_lock_nu1[ei][mi][oi]) 192 193 values[ei][si][mi].append(zeros(num, float64)) 194 errors[ei][si][mi].append(ones(num, float64)) 195 missing[ei][si][mi].append(zeros(num, int32)) 196 197 # Reconstruct the structures for interpolation. 198 else: 199 recalc_tau = False 200 values = [] 201 errors = [] 202 missing = [] 203 for exp_type, ei in loop_exp(return_indices=True): 204 values.append([]) 205 errors.append([]) 206 missing.append([]) 207 for si in range(len(spins)): 208 values[ei].append([]) 209 errors[ei].append([]) 210 missing[ei].append([]) 211 for frq, mi in loop_frq(return_indices=True): 212 values[ei][si].append([]) 213 errors[ei][si].append([]) 214 missing[ei][si].append([]) 215 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 216 if exp_type in EXP_TYPE_LIST_CPMG: 217 num = len(cpmg_frqs[ei][mi][oi]) 218 else: 219 num = len(spin_lock_nu1[ei][mi][oi]) 220 values[ei][si][mi].append(zeros(num, float64)) 221 errors[ei][si][mi].append(ones(num, float64)) 222 missing[ei][si][mi].append(zeros(num, int32)) 223 224 # Initialise the relaxation dispersion fit functions. 225 model = Dispersion(model=spins[0].model, num_params=param_num(spins=spins), num_spins=len(spins), 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, recalc_tau=recalc_tau, r1_fit=r1_fit) 226 227 # Make a single function call. This will cause back calculation and the data will be stored in the class instance. 228 chi2 = model.func(param_vector) 229 230 # Store the chi-squared value. 231 if store_chi2: 232 for spin in spins: 233 spin.chi2 = chi2 234 235 # Return the structure. 236 return model.get_back_calc()

237 238

239 -def calculate_r2eff():

240 """Calculate the R2eff values for fixed relaxation time period data.""" 241 242 # Data checks. 243 check_exp_type() 244 check_disp_points() 245 check_exp_type_fixed_time() 246 247 # Printouts. 248 print("Calculating the R2eff/R1rho values for fixed relaxation time period data.") 249 250 # Loop over the spins. 251 for spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): 252 # Spin ID printout. 253 print("Spin '%s'." % spin_id) 254 255 # Skip spins which have no data. 256 if not hasattr(spin, 'peak_intensity'): 257 continue 258 259 # Initialise the data structures. 260 if not hasattr(spin, 'r2eff'): 261 spin.r2eff = {} 262 if not hasattr(spin, 'r2eff_err'): 263 spin.r2eff_err = {} 264 265 # Loop over all the data. 266 for exp_type, frq, offset, point, time in loop_exp_frq_offset_point_time(): 267 # The three keys. 268 ref_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=None, time=time) 269 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 270 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 271 272 # Check for missing data. 273 missing = False 274 for i in range(len(ref_keys)): 275 if ref_keys[i] not in spin.peak_intensity: 276 missing = True 277 for i in range(len(int_keys)): 278 if int_keys[i] not in spin.peak_intensity: 279 missing = True 280 if missing: 281 continue 282 283 # Average the reference intensity data and errors. 284 ref_intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time) 285 ref_intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=None, time=time, error=True) 286 287 # Average the intensity data and errors. 288 intensity = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 289 intensity_err = average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True) 290 291 # Check for math domain errors or log for values less than 0.0. 292 if ref_intensity == 0.0: 293 skip_data = True 294 elif float(intensity) / ref_intensity <= 0.0: 295 skip_data = True 296 else: 297 skip_data = False 298 299 if skip_data: 300 spin_string = generate_spin_string(spin=spin, mol_name=mol_name, res_num=resi, res_name=resn) 301 msg = "Math log(I / I_ref) domain error for spin '%s' in R2eff value calculation for fixed relaxation time period data. I=%3.3f, I_ref=%3.3f. The point is skipped." % (spin_string, intensity, ref_intensity) 302 warn(RelaxWarning("%s" % msg)) 303 point_info = "This happened for '%s' at %3.1f MHz, for offset=%3.1f ppm and dispersion point %3.1f Hz and time %1.2f s.\n" % (exp_type, frq/1E6, offset, point, time) 304 print(point_info) 305 else: 306 # Calculate the R2eff value. 307 spin.r2eff[param_key] = calc_two_point_r2eff(relax_time=time, I_ref=ref_intensity, I=intensity) 308 309 # Calculate the R2eff error. 310 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)

311 312

313 -def minimise_r2eff(spins=None, spin_ids=None, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling_matrix=None, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):

314 """Optimise the R2eff model by fitting the 2-parameter exponential curves. 315 316 This mimics the R1 and R2 relax_fit analysis. 317 318 319 @keyword spins: The list of spins for the cluster. 320 @type spins: list of SpinContainer instances 321 @keyword spin_ids: The list of spin IDs for the cluster. 322 @type spin_ids: list of str 323 @keyword min_algor: The minimisation algorithm to use. 324 @type min_algor: str 325 @keyword min_options: An array of options to be used by the minimisation algorithm. 326 @type min_options: array of str 327 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 328 @type func_tol: None or float 329 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 330 @type grad_tol: None or float 331 @keyword max_iterations: The maximum number of iterations for the algorithm. 332 @type max_iterations: int 333 @keyword constraints: If True, constraints are used during optimisation. 334 @type constraints: bool 335 @keyword scaling_matrix: The diagonal and square scaling matrix. 336 @type scaling_matrix: numpy rank-2, float64 array or None 337 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 338 @type verbosity: int 339 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 340 @type sim_index: None or int 341 @keyword lower: The model specific 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. 342 @type lower: list of numbers 343 @keyword upper: The model specific 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. 344 @type upper: list of numbers 345 @keyword inc: The model specific 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. 346 @type inc: list of int 347 """ 348 349 # Check that the C modules have been compiled. 350 if not C_module_exp_fn: 351 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 352 353 # Loop over the spins. 354 for si in range(len(spins)): 355 # Skip deselected spins. 356 if not spins[si].select: 357 continue 358 359 # Loop over each spectrometer frequency and dispersion point. 360 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 361 # The parameter key. 362 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 363 364 # The initial parameter vector. 365 param_vector = assemble_param_vector(spins=[spins[si]], key=param_key, sim_index=sim_index) 366 367 # Diagonal scaling. 368 if scaling_matrix is not None: 369 param_vector = dot(inv(scaling_matrix), param_vector) 370 371 # Linear constraints. 372 A, b = None, None 373 if constraints: 374 A, b = linear_constraints(spins=[spins[si]], scaling_matrix=scaling_matrix) 375 376 # Print out. 377 if verbosity >= 1: 378 # Individual spin section. 379 top = 2 380 if verbosity >= 2: 381 top += 2 382 text = "Fitting to spin %s, frequency %s and dispersion point %s" % (spin_ids[si], frq, point) 383 subsection(file=sys.stdout, text=text, prespace=top) 384 385 # Grid search printout. 386 if match('^[Gg]rid', min_algor): 387 result = 1 388 for x in inc: 389 result = mul(result, x) 390 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % result) 391 392 # The peak intensities, errors and times. 393 values = [] 394 errors = [] 395 times = [] 396 data_flag = True 397 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 398 # Check the peak intensity keys. 399 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 400 peak_intensities = spins[si].peak_intensity 401 if sim_index != None: 402 peak_intensities = spins[si].peak_intensity_sim 403 for i in range(len(int_keys)): 404 if int_keys[i] not in peak_intensities: 405 if verbosity: 406 warn(RelaxWarning("The spin %s peak intensity key '%s' is not present, skipping the optimisation." % (spin_ids[si], int_keys[i]))) 407 data_flag = False 408 break 409 410 if data_flag: 411 values.append(average_intensity(spin=spins[si], exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, sim_index=sim_index)) 412 errors.append(average_intensity(spin=spins[si], exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 413 times.append(time) 414 if not data_flag: 415 continue 416 417 # Raise errors if number of time points is less than 2. 418 if len(times) < 3: 419 subsection(file=sys.stdout, text="Exponential curve fitting error for point:", prespace=2) 420 point_info = "%s at %3.1f MHz, for offset=%3.3f ppm and dispersion point %-5.1f, with %i time points." % (exp_type, frq/1E6, offset, point, len(times)) 421 raise RelaxError("The data setup points to exponential curve fitting, but only %i time points was found, where 3 time points is minimum. If calculating R2eff values for fixed relaxation time period data, check that a reference intensity has been specified for each offset value."%(len(times))) 422 423 # The scaling matrix in a diagonalised list form. 424 scaling_list = [] 425 if scaling_matrix is None: 426 for i in range(len(param_vector)): 427 scaling_list.append(1.0) 428 else: 429 for i in range(len(scaling_matrix)): 430 scaling_list.append(scaling_matrix[i, i]) 431 432 # Initialise the function to minimise. 433 model = Relax_fit_opt(model='exp', num_params=len(param_vector), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 434 435 # Grid search. 436 if search('^[Gg]rid', min_algor): 437 results = grid(func=model.func, args=(), num_incs=inc, lower=lower, upper=upper, A=A, b=b, verbosity=verbosity) 438 439 # Unpack the results. 440 param_vector, chi2, iter_count, warning = results 441 f_count = iter_count 442 g_count = 0.0 443 h_count = 0.0 444 445 # Minimisation. 446 else: 447 results = generic_minimise(func=model.func, dfunc=model.dfunc, d2func=model.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) 448 449 # Unpack the results. 450 if results == None: 451 return 452 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 453 454 # Scaling. 455 if scaling_matrix is not None: 456 param_vector = dot(scaling_matrix, param_vector) 457 458 # Disassemble the parameter vector. 459 disassemble_param_vector(param_vector=param_vector, spins=[spins[si]], key=param_key, sim_index=sim_index) 460 461 # Monte Carlo minimisation statistics. 462 if sim_index != None: 463 # Chi-squared statistic. 464 spins[si].chi2_sim[sim_index] = chi2 465 466 # Iterations. 467 spins[si].iter_sim[sim_index] = iter_count 468 469 # Function evaluations. 470 spins[si].f_count_sim[sim_index] = f_count 471 472 # Gradient evaluations. 473 spins[si].g_count_sim[sim_index] = g_count 474 475 # Hessian evaluations. 476 spins[si].h_count_sim[sim_index] = h_count 477 478 # Warning. 479 spins[si].warning_sim[sim_index] = warning 480 481 # Normal statistics. 482 else: 483 # Chi-squared statistic. 484 spins[si].chi2 = chi2 485 486 # Iterations. 487 spins[si].iter = iter_count 488 489 # Function evaluations. 490 spins[si].f_count = f_count 491 492 # Gradient evaluations. 493 spins[si].g_count = g_count 494 495 # Hessian evaluations. 496 spins[si].h_count = h_count 497 498 # Warning. 499 spins[si].warning = warning

500 501 502

503 -class Disp_memo(Memo):

504 """The relaxation dispersion memo class.""" 505

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

507 """Initialise the relaxation dispersion memo class. 508 509 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. 510 511 512 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored. 513 @type spins: list of SpinContainer instances 514 @keyword spin_ids: The spin ID strings for the cluster. 515 @type spin_ids: list of str 516 @keyword sim_index: The optional MC simulation index. 517 @type sim_index: int 518 @keyword scaling_matrix: The diagonal, square scaling matrix. 519 @type scaling_matrix: numpy diagonal matrix 520 @keyword verbosity: The verbosity level. This is used by the result command returned to the master for printouts. 521 @type verbosity: int 522 """ 523 524 # Execute the base class __init__() method. 525 super(Disp_memo, self).__init__() 526 527 # Store the arguments. 528 self.spins = spins 529 self.spin_ids = spin_ids 530 self.sim_index = sim_index 531 self.scaling_matrix = scaling_matrix 532 self.verbosity = verbosity

533 534 535

536 -class Disp_minimise_command(Slave_command):

537 """Command class for relaxation dispersion optimisation on the slave processor.""" 538

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

540 """Initialise the base class, storing all the master data to be sent to the slave processor. 541 542 This method is run on the master processor whereas the run() method is run on the slave processor. 543 544 545 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored. 546 @type spins: list of SpinContainer instances 547 @keyword spin_ids: The list of spin ID strings corresponding to the spins argument. 548 @type spin_ids: list of str 549 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 550 @type sim_index: None or int 551 @keyword scaling_matrix: The diagonal, square scaling matrix. 552 @type scaling_matrix: numpy diagonal matrix 553 @keyword min_algor: The minimisation algorithm to use. 554 @type min_algor: str 555 @keyword min_options: An array of options to be used by the minimisation algorithm. 556 @type min_options: array of str 557 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 558 @type func_tol: None or float 559 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 560 @type grad_tol: None or float 561 @keyword max_iterations: The maximum number of iterations for the algorithm. 562 @type max_iterations: int 563 @keyword constraints: If True, constraints are used during optimisation. 564 @type constraints: bool 565 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 566 @type verbosity: int 567 @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. 568 @type lower: array of numbers 569 @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. 570 @type upper: array of numbers 571 @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. 572 @type inc: array of int 573 @keyword fields: The list of unique of spectrometer field strengths. 574 @type fields: int 575 @keyword param_names: The list of parameter names to use in printouts. 576 @type param_names: str 577 """ 578 579 # Execute the base class __init__() method. 580 super(Disp_minimise_command, self).__init__() 581 582 # Store the arguments needed by the run() method. 583 self.spins = spins 584 self.spin_ids = spin_ids 585 self.sim_index = sim_index 586 self.scaling_matrix = scaling_matrix 587 self.verbosity = verbosity 588 self.min_algor = min_algor 589 self.min_options = min_options 590 self.func_tol = func_tol 591 self.grad_tol = grad_tol 592 self.max_iterations = max_iterations 593 self.lower = lower 594 self.upper = upper 595 self.inc = inc 596 self.fields = fields 597 self.param_names = param_names 598 599 # Create the initial parameter vector. 600 self.param_vector = assemble_param_vector(spins=self.spins) 601 if len(scaling_matrix): 602 self.param_vector = dot(inv(scaling_matrix), self.param_vector) 603 604 # Linear constraints. 605 self.A, self.b = None, None 606 if constraints: 607 self.A, self.b = linear_constraints(spins=spins, scaling_matrix=scaling_matrix) 608 609 # Test if the spectrometer frequencies have been set. 610 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'): 611 raise RelaxError("The spectrometer frequency information has not been specified.") 612 613 # The R2eff/R1rho data. 614 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) 615 616 # The offset and R1 data. 617 r1_setup() 618 self.offsets, spin_lock_fields_inter, self.chemical_shifts, self.tilt_angles, self.Delta_omega, self.w_eff = return_offset_data(spins=spins, spin_ids=spin_ids, field_count=len(fields)) 619 self.r1 = return_r1_data(spins=spins, spin_ids=spin_ids, field_count=len(fields), sim_index=sim_index) 620 self.r1_fit = is_r1_optimised(spins[0].model) 621 622 # Parameter number. 623 self.param_num = param_num(spins=spins) 624 625 # The dispersion data. 626 self.dispersion_points = cdp.dispersion_points 627 self.cpmg_frqs = return_cpmg_frqs(ref_flag=False) 628 self.spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False)

629 630

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

632 """Set up and perform the optimisation.""" 633 634 # Print out. 635 if self.verbosity >= 1: 636 # Individual spin block section. 637 top = 2 638 if self.verbosity >= 2: 639 top += 2 640 subsection(file=sys.stdout, text="Fitting to the spin block %s"%self.spin_ids, prespace=top) 641 642 # Grid search printout. 643 if search('^[Gg]rid', self.min_algor): 644 result = 1 645 for x in self.inc: 646 result = mul(result, x) 647 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % result) 648 649 # Initialise the function to minimise. 650 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, r1_fit=self.r1_fit) 651 652 # Grid search. 653 if search('^[Gg]rid', self.min_algor): 654 results = grid(func=model.func, args=(), num_incs=self.inc, lower=self.lower, upper=self.upper, A=self.A, b=self.b, verbosity=self.verbosity) 655 656 # Unpack the results. 657 param_vector, chi2, iter_count, warning = results 658 f_count = iter_count 659 g_count = 0.0 660 h_count = 0.0 661 662 # Minimisation. 663 else: 664 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) 665 666 # Unpack the results. 667 if results == None: 668 return 669 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 670 671 # Optimisation printout. 672 if self.verbosity: 673 print("\nOptimised parameter values:") 674 for i in range(len(param_vector)): 675 print("%-20s %25.15f" % (self.param_names[i], param_vector[i]*self.scaling_matrix[i, i])) 676 677 # Create the result command object to send back to the master. 678 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.get_back_calc(), completed=False))

679 680 681

682 -class Disp_result_command(Result_command):

683 """Class for processing the dispersion optimisation results. 684 685 This object will be sent from the slave back to the master to have its run() method executed. 686 """ 687

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

689 """Set up this class object on the slave, placing the minimisation results here. 690 691 @keyword processor: The processor object. 692 @type processor: multi.processor.Processor instance 693 @keyword memo_id: The memo identification string. 694 @type memo_id: str 695 @keyword param_vector: The optimised parameter vector. 696 @type param_vector: numpy rank-1 array 697 @keyword chi2: The final target function value. 698 @type chi2: float 699 @keyword iter_count: The number of optimisation iterations. 700 @type iter_count: int 701 @keyword f_count: The total function call count. 702 @type f_count: int 703 @keyword g_count: The total gradient call count. 704 @type g_count: int 705 @keyword h_count: The total Hessian call count. 706 @type h_count: int 707 @keyword warning: Any optimisation warnings. 708 @type warning: str or None 709 @keyword missing: The data structure indicating which R2eff/R1rho' base data is missing. 710 @type missing: numpy rank-3 array 711 @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. 712 @type back_calc: numpy rank-3 array 713 @keyword completed: A flag which if True signals that the optimisation successfully completed. 714 @type completed: bool 715 """ 716 717 # Execute the base class __init__() method. 718 super(Disp_result_command, self).__init__(processor=processor, completed=completed) 719 720 # Store the arguments (to be sent back to the master). 721 self.memo_id = memo_id 722 self.param_vector = param_vector 723 self.chi2 = chi2 724 self.iter_count = iter_count 725 self.f_count = f_count 726 self.g_count = g_count 727 self.h_count = h_count 728 self.warning = warning 729 self.missing = missing 730 self.back_calc = back_calc 731 self.completed = completed

732 733

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

735 """Disassemble the model-free optimisation results (on the master). 736 737 @param processor: Unused! 738 @type processor: None 739 @param memo: The dispersion memo. This holds a lot of the data and objects needed for processing the results from the slave. 740 @type memo: memo 741 """ 742 743 # Printout. 744 if memo.sim_index != None: 745 print("Simulation %s, cluster %s" % (memo.sim_index+1, memo.spin_ids)) 746 747 # Scaling. 748 if memo.scaling_matrix is not None: 749 param_vector = dot(memo.scaling_matrix, self.param_vector) 750 751 # Disassemble the parameter vector. 752 disassemble_param_vector(param_vector=param_vector, spins=memo.spins, sim_index=memo.sim_index) 753 param_conversion(spins=memo.spins, sim_index=memo.sim_index) 754 755 # Monte Carlo minimisation statistics. 756 if memo.sim_index != None: 757 for spin in memo.spins: 758 # Skip deselected spins. 759 if not spin.select: 760 continue 761 762 # Chi-squared statistic. 763 spin.chi2_sim[memo.sim_index] = self.chi2 764 765 # Iterations. 766 spin.iter_sim[memo.sim_index] = self.iter_count 767 768 # Function evaluations. 769 spin.f_count_sim[memo.sim_index] = self.f_count 770 771 # Gradient evaluations. 772 spin.g_count_sim[memo.sim_index] = self.g_count 773 774 # Hessian evaluations. 775 spin.h_count_sim[memo.sim_index] = self.h_count 776 777 # Warning. 778 spin.warning_sim[memo.sim_index] = self.warning 779 780 # Normal statistics. 781 else: 782 for spin in memo.spins: 783 # Skip deselected spins. 784 if not spin.select: 785 continue 786 787 # Chi-squared statistic. 788 spin.chi2 = self.chi2 789 790 # Iterations. 791 spin.iter = self.iter_count 792 793 # Function evaluations. 794 spin.f_count = self.f_count 795 796 # Gradient evaluations. 797 spin.g_count = self.g_count 798 799 # Hessian evaluations. 800 spin.h_count = self.h_count 801 802 # Warning. 803 spin.warning = self.warning 804 805 # Store the back-calculated values. 806 if memo.sim_index == None: 807 # 1H MMQ flag. 808 proton_mmq_flag = has_proton_mmq_cpmg() 809 810 # Loop over each spin, packing the data. 811 si = 0 812 for spin_index in range(len(memo.spins)): 813 # Skip deselected spins. 814 if not memo.spins[spin_index].select: 815 continue 816 817 # Pack the data. 818 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) 819 820 # Increment the spin index. 821 si += 1

822

Source Code for Module specific_analyses.relax_disp.optimisation