specific_analyses.relax_disp.optimisation

1 ############################################################################### 2 # # 3 # Copyright (C) 2013-2014 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_LIST_MMQ, 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, 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 # The peak intensities and times. 80 values = [] 81 errors = [] 82 times = [] 83 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 84 # The data. 85 values.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) 86 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 87 times.append(time) 88 89 # The scaling matrix in a diagonalised list form. 90 scaling_list = [] 91 for i in range(len(param_vector)): 92 scaling_list.append(1.0) 93 94 # Initialise the relaxation fit functions. 95 model = Relax_fit_opt(model='exp', num_params=len(spin.params), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 96 97 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 98 model.func(param_vector) 99 100 # Get the data back. 101 results = model.back_calc_data() 102 103 # Return the correct peak height. 104 return results

105 106

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

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

225 226

227 -def calculate_r2eff():

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

299 300

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

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

473 474 475

476 -class Disp_memo(Memo):

477 """The relaxation dispersion memo class.""" 478

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

480 """Initialise the relaxation dispersion memo class. 481 482 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. 483 484 485 @keyword spins: The list of spin data container for the cluster. If this argument is supplied, then the spin_id argument will be ignored. 486 @type spins: list of SpinContainer instances 487 @keyword spin_ids: The spin ID strings for the cluster. 488 @type spin_ids: list of str 489 @keyword sim_index: The optional MC simulation index. 490 @type sim_index: int 491 @keyword scaling_matrix: The diagonal, square scaling matrix. 492 @type scaling_matrix: numpy diagonal matrix 493 @keyword verbosity: The verbosity level. This is used by the result command returned to the master for printouts. 494 @type verbosity: int 495 """ 496 497 # Execute the base class __init__() method. 498 super(Disp_memo, self).__init__() 499 500 # Store the arguments. 501 self.spins = spins 502 self.spin_ids = spin_ids 503 self.sim_index = sim_index 504 self.scaling_matrix = scaling_matrix 505 self.verbosity = verbosity

506 507 508

509 -class Disp_minimise_command(Slave_command):

510 """Command class for relaxation dispersion optimisation on the slave processor.""" 511

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

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

602 603

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

605 """Set up and perform the optimisation.""" 606 607 # Print out. 608 if self.verbosity >= 1: 609 # Individual spin block section. 610 top = 2 611 if self.verbosity >= 2: 612 top += 2 613 subsection(file=sys.stdout, text="Fitting to the spin block %s"%self.spin_ids, prespace=top) 614 615 # Grid search printout. 616 if search('^[Gg]rid', self.min_algor): 617 result = 1 618 for x in self.inc: 619 result = mul(result, x) 620 print("Unconstrained grid search size: %s (constraints may decrease this size).\n" % result) 621 622 # Initialise the function to minimise. 623 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) 624 625 # Grid search. 626 if search('^[Gg]rid', self.min_algor): 627 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) 628 629 # Unpack the results. 630 param_vector, chi2, iter_count, warning = results 631 f_count = iter_count 632 g_count = 0.0 633 h_count = 0.0 634 635 # Minimisation. 636 else: 637 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) 638 639 # Unpack the results. 640 if results == None: 641 return 642 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 643 644 # Optimisation printout. 645 if self.verbosity: 646 print("\nOptimised parameter values:") 647 for i in range(len(param_vector)): 648 print("%-20s %25.15f" % (self.param_names[i], param_vector[i]*self.scaling_matrix[i, i])) 649 650 # Create the result command object to send back to the master. 651 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))

652 653 654

655 -class Disp_result_command(Result_command):

656 """Class for processing the dispersion optimisation results. 657 658 This object will be sent from the slave back to the master to have its run() method executed. 659 """ 660

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

662 """Set up this class object on the slave, placing the minimisation results here. 663 664 @keyword processor: The processor object. 665 @type processor: multi.processor.Processor instance 666 @keyword memo_id: The memo identification string. 667 @type memo_id: str 668 @keyword param_vector: The optimised parameter vector. 669 @type param_vector: numpy rank-1 array 670 @keyword chi2: The final target function value. 671 @type chi2: float 672 @keyword iter_count: The number of optimisation iterations. 673 @type iter_count: int 674 @keyword f_count: The total function call count. 675 @type f_count: int 676 @keyword g_count: The total gradient call count. 677 @type g_count: int 678 @keyword h_count: The total Hessian call count. 679 @type h_count: int 680 @keyword warning: Any optimisation warnings. 681 @type warning: str or None 682 @keyword missing: The data structure indicating which R2eff/R1rho' base data is missing. 683 @type missing: numpy rank-3 array 684 @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. 685 @type back_calc: numpy rank-3 array 686 @keyword completed: A flag which if True signals that the optimisation successfully completed. 687 @type completed: bool 688 """ 689 690 # Execute the base class __init__() method. 691 super(Disp_result_command, self).__init__(processor=processor, completed=completed) 692 693 # Store the arguments (to be sent back to the master). 694 self.memo_id = memo_id 695 self.param_vector = param_vector 696 self.chi2 = chi2 697 self.iter_count = iter_count 698 self.f_count = f_count 699 self.g_count = g_count 700 self.h_count = h_count 701 self.warning = warning 702 self.missing = missing 703 self.back_calc = back_calc 704 self.completed = completed

705 706

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

708 """Disassemble the model-free optimisation results (on the master). 709 710 @param processor: Unused! 711 @type processor: None 712 @param memo: The dispersion memo. This holds a lot of the data and objects needed for processing the results from the slave. 713 @type memo: memo 714 """ 715 716 # Printout. 717 if memo.sim_index != None: 718 print("Simulation %s, cluster %s" % (memo.sim_index+1, memo.spin_ids)) 719 720 # Scaling. 721 if memo.scaling_matrix != None: 722 param_vector = dot(memo.scaling_matrix, self.param_vector) 723 724 # Disassemble the parameter vector. 725 disassemble_param_vector(param_vector=param_vector, spins=memo.spins, sim_index=memo.sim_index) 726 param_conversion(spins=memo.spins, sim_index=memo.sim_index) 727 728 # Monte Carlo minimisation statistics. 729 if memo.sim_index != None: 730 for spin in memo.spins: 731 # Skip deselected spins. 732 if not spin.select: 733 continue 734 735 # Chi-squared statistic. 736 spin.chi2_sim[memo.sim_index] = self.chi2 737 738 # Iterations. 739 spin.iter_sim[memo.sim_index] = self.iter_count 740 741 # Function evaluations. 742 spin.f_count_sim[memo.sim_index] = self.f_count 743 744 # Gradient evaluations. 745 spin.g_count_sim[memo.sim_index] = self.g_count 746 747 # Hessian evaluations. 748 spin.h_count_sim[memo.sim_index] = self.h_count 749 750 # Warning. 751 spin.warning_sim[memo.sim_index] = self.warning 752 753 # Normal statistics. 754 else: 755 for spin in memo.spins: 756 # Skip deselected spins. 757 if not spin.select: 758 continue 759 760 # Chi-squared statistic. 761 spin.chi2 = self.chi2 762 763 # Iterations. 764 spin.iter = self.iter_count 765 766 # Function evaluations. 767 spin.f_count = self.f_count 768 769 # Gradient evaluations. 770 spin.g_count = self.g_count 771 772 # Hessian evaluations. 773 spin.h_count = self.h_count 774 775 # Warning. 776 spin.warning = self.warning 777 778 # Store the back-calculated values. 779 if memo.sim_index == None: 780 # 1H MMQ flag. 781 proton_mmq_flag = has_proton_mmq_cpmg() 782 783 # Loop over each spin, packing the data. 784 si = 0 785 for spin_index in range(len(memo.spins)): 786 # Skip deselected spins. 787 if not memo.spins[spin_index].select: 788 continue 789 790 # Pack the data. 791 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) 792 793 # Increment the spin index. 794 si += 1

795

Source Code for Module specific_analyses.relax_disp.optimisation