specific_analyses.relax

51 """Class containing functions for relaxation curve fitting.""" 52 53 # Class variable for storing the class instance (for the singleton design pattern). 54 instance = None 55

56 - def __init__(self):

57 """Initialise the class by placing API_common methods into the API.""" 58 59 # Place methods into the API. 60 self.base_data_loop = self._base_data_loop_spin 61 self.model_loop = self._model_loop_spin 62 self.print_model_title = self._print_model_title_spin 63 self.return_conversion_factor = self._return_no_conversion_factor 64 self.return_value = self._return_value_general 65 self.set_error = self._set_error_spin 66 self.set_param_values = self._set_param_values_spin 67 self.set_selected_sim = self._set_selected_sim_spin 68 self.sim_init_values = self._sim_init_values_spin 69 self.sim_return_param = self._sim_return_param_spin 70 self.sim_return_selected = self._sim_return_selected_spin 71 72 # Place a copy of the parameter list object in the instance namespace. 73 self._PARAMS = Relax_fit_params()

74 75

76 - def covariance_matrix(self, model_info=None, verbosity=1):

77 """Return the Jacobian and weights required for parameter errors via the covariance matrix. 78 79 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method. 80 @type model_info: SpinContainer instance, str 81 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 82 @type verbosity: int 83 @return: The Jacobian and weight matrices for the given model. 84 @rtype: numpy rank-2 array, numpy rank-2 array 85 """ 86 87 # Unpack the data. 88 spin, spin_id = model_info 89 90 # Check that the C modules have been compiled. 91 if not C_module_exp_fn: 92 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 93 94 # Raise Error, if not optimised. 95 if not (hasattr(spin, 'rx') and hasattr(spin, 'i0')): 96 raise RelaxError("Spin '%s' does not contain optimised 'rx' and 'i0' values. Try execute: minimise.execute(min_algor='Newton', constraints=False)"%(spin_id)) 97 98 # Raise warning, if gradient count is 0. This could point to a lack of minimisation first. 99 if hasattr(spin, 'g_count'): 100 if getattr(spin, 'g_count') == 0.0: 101 text = "Spin %s contains a gradient count of 0.0. Is the rx parameter optimised? Try execute: minimise.execute(min_algor='Newton', constraints=False)" %(spin_id) 102 warn(RelaxWarning("%s." % text)) 103 104 # Print information. 105 if verbosity >= 1: 106 # Individual spin block section. 107 top = 2 108 if verbosity >= 2: 109 top += 2 110 subsection(file=sys.stdout, text="Estimating rx error for spin: %s"%spin_id, prespace=top) 111 112 # The peak intensities and times. 113 values = [] 114 errors = [] 115 times = [] 116 for key in spin.peak_intensity: 117 values.append(spin.peak_intensity[key]) 118 errors.append(spin.peak_intensity_err[key]) 119 times.append(cdp.relax_times[key]) 120 121 # Convert to numpy array. 122 values = asarray(values) 123 errors = asarray(errors) 124 times = asarray(times) 125 126 # Create the parameter vector and scaling matrix (as a diagonalised list). 127 param_vector = assemble_param_vector(spin=spin) 128 scaling_list = [] 129 for i in range(len(spin.params)): 130 scaling_list.append(1.0) 131 132 # Initialise data in C code. 133 model = Relax_fit_opt(model=spin.model, num_params=len(param_vector), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 134 135 # Use the direct Jacobian from function. 136 jacobian_matrix_exp = transpose(asarray( model.jacobian(param_vector) ) ) 137 weights = 1. / errors**2 138 139 # Return the matrices. 140 return jacobian_matrix_exp, weights

141 142

143 - def create_mc_data(self, data_id=None):

144 """Create the Monte Carlo peak intensity data. 145 146 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 147 @type data_id: str 148 @return: The Monte Carlo simulation data. 149 @rtype: list of floats 150 """ 151 152 # Initialise the MC data data structure. 153 mc_data = {} 154 155 # Get the spin container. 156 spin = return_spin(data_id) 157 158 # Skip deselected spins. 159 if not spin.select: 160 return 161 162 # Skip spins which have no data. 163 if not hasattr(spin, 'peak_intensity'): 164 return 165 166 # Test if the model is set. 167 if not hasattr(spin, 'model') or not spin.model: 168 raise RelaxNoModelError 169 170 # Loop over the spectral time points. 171 for id in cdp.relax_times: 172 # Back calculate the value. 173 value = back_calc(spin=spin, relax_time_id=id) 174 175 # Append the value. 176 mc_data[id] = value 177 178 # Return the MC data. 179 return mc_data

180 181

182 - def data_init(self, data, sim=False):

183 """Initialise the spin specific data structures. 184 185 @param data: The spin ID string from the _base_data_loop_spin() method. 186 @type data: str 187 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure. 188 @type sim: bool 189 """ 190 191 # Get the spin container. 192 spin = return_spin(data) 193 194 # Loop over the data structure names. 195 for name in self.data_names(set='params'): 196 # Data structures which are initially empty arrays. 197 list_data = [ 'params' ] 198 if name in list_data: 199 init_data = [] 200 201 # Otherwise initialise the data structure to None. 202 else: 203 init_data = None 204 205 # If the name is not in the spin container, add it. 206 if not hasattr(spin, name): 207 setattr(spin, name, init_data)

208 209

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

211 """Return a vector of parameter names. 212 213 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method. 214 @type model_info: SpinContainer instance, str 215 @return: The vector of parameter names. 216 @rtype: list of str 217 """ 218 219 # Check that the model is setup. 220 check_model_setup() 221 222 # Unpack the data. 223 spin, spin_id = model_info 224 225 # Return the parameter names. 226 return spin.params

227 228

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

230 """Return a vector of parameter values. 231 232 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method. 233 @type model_info: SpinContainer instance, str 234 @keyword sim_index: The optional Monte Carlo simulation index. 235 @type sim_index: int 236 @return: The vector of parameter values. 237 @rtype: list of str 238 """ 239 240 # Unpack the data. 241 spin, spin_id = model_info 242 243 # Return the vector. 244 return assemble_param_vector(spin=spin, sim_index=sim_index)

245 246

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

248 """The exponential curve fitting grid search method. 249 250 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model. 251 @type lower: list of lists of numbers 252 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model. 253 @type upper: list of lists of numbers 254 @keyword inc: The per-model 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. 255 @type inc: list of lists of int 256 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 257 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 258 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 259 @type constraints: bool 260 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 261 @type verbosity: int 262 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 263 @type sim_index: int 264 """ 265 266 # Minimisation. 267 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

268 269

270 - def minimise(self, 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):

271 """Relaxation curve fitting minimisation method. 272 273 @keyword min_algor: The minimisation algorithm to use. 274 @type min_algor: str 275 @keyword min_options: An array of options to be used by the minimisation algorithm. 276 @type min_options: array of str 277 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 278 @type func_tol: None or float 279 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 280 @type grad_tol: None or float 281 @keyword max_iterations: The maximum number of iterations for the algorithm. 282 @type max_iterations: int 283 @keyword constraints: If True, constraints are used during optimisation. 284 @type constraints: bool 285 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 286 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 287 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 288 @type verbosity: int 289 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 290 @type sim_index: None or int 291 @keyword lower: The per-model 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. 292 @type lower: list of lists of numbers 293 @keyword upper: The per-model 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. 294 @type upper: list of lists of numbers 295 @keyword inc: The per-model 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. 296 @type inc: list of lists of int 297 """ 298 299 # Checks. 300 check_mol_res_spin_data() 301 302 # Loop over the sequence. 303 model_index = 0 304 for spin, spin_id in self.model_loop(): 305 # Skip deselected spins. 306 if not spin.select: 307 continue 308 309 # Skip spins which have no data. 310 if not hasattr(spin, 'peak_intensity'): 311 continue 312 313 # Create the initial parameter vector. 314 param_vector = assemble_param_vector(spin=spin) 315 316 # Diagonal scaling. 317 if scaling_matrix[model_index] is not None: 318 param_vector = dot(inv(scaling_matrix[model_index]), param_vector) 319 320 # Linear constraints. 321 if constraints: 322 A, b = linear_constraints(spin=spin, scaling_matrix=scaling_matrix[model_index]) 323 else: 324 A, b = None, None 325 326 # Print out. 327 if verbosity >= 1: 328 # Individual spin printout. 329 if verbosity >= 2: 330 print("\n\n") 331 332 string = "Fitting to spin " + repr(spin_id) 333 print("\n\n" + string) 334 print(len(string) * '~') 335 336 337 # Initialise the function to minimise. 338 ###################################### 339 340 # The peak intensities and times. 341 values = [] 342 errors = [] 343 times = [] 344 for key in spin.peak_intensity: 345 # The values. 346 if sim_index == None: 347 values.append(spin.peak_intensity[key]) 348 else: 349 values.append(spin.peak_intensity_sim[sim_index][key]) 350 351 # The errors. 352 errors.append(spin.peak_intensity_err[key]) 353 354 # The relaxation times. 355 times.append(cdp.relax_times[key]) 356 357 # The scaling matrix in a diagonalised list form. 358 scaling_list = [] 359 if scaling_matrix[model_index] is None: 360 for i in range(len(param_vector)): 361 scaling_list.append(1.0) 362 else: 363 for i in range(len(scaling_matrix[model_index])): 364 scaling_list.append(scaling_matrix[model_index][i, i]) 365 366 # Set up the target function. 367 model = Relax_fit_opt(model=spin.model, num_params=len(spin.params), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 368 369 370 # Setup the minimisation algorithm when constraints are present. 371 ################################################################ 372 373 if constraints and not match('^[Gg]rid', min_algor): 374 algor = min_options[0] 375 else: 376 algor = min_algor 377 378 379 # Levenberg-Marquardt minimisation. 380 ################################### 381 382 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 383 # Reconstruct the error data structure. 384 lm_error = zeros(len(spin.relax_times), float64) 385 index = 0 386 for k in range(len(spin.relax_times)): 387 lm_error[index:index+len(relax_error[k])] = relax_error[k] 388 index = index + len(relax_error[k]) 389 390 min_options = min_options + (self.relax_fit.lm_dri, lm_error) 391 392 393 # Minimisation. 394 ############### 395 396 # Grid search. 397 if search('^[Gg]rid', min_algor): 398 results = grid(func=model.func, args=(), num_incs=inc[model_index], lower=lower[model_index], upper=upper[model_index], A=A, b=b, verbosity=verbosity) 399 400 # Unpack the results. 401 param_vector, chi2, iter_count, warning = results 402 f_count = iter_count 403 g_count = 0.0 404 h_count = 0.0 405 406 # Minimisation. 407 else: 408 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) 409 410 # Unpack the results. 411 if results == None: 412 return 413 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 414 415 # Scaling. 416 if scaling_matrix[model_index] is not None: 417 param_vector = dot(scaling_matrix[model_index], param_vector) 418 419 # Disassemble the parameter vector. 420 disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index) 421 422 # Monte Carlo minimisation statistics. 423 if sim_index != None: 424 # Chi-squared statistic. 425 spin.chi2_sim[sim_index] = chi2 426 427 # Iterations. 428 spin.iter_sim[sim_index] = iter_count 429 430 # Function evaluations. 431 spin.f_count_sim[sim_index] = f_count 432 433 # Gradient evaluations. 434 spin.g_count_sim[sim_index] = g_count 435 436 # Hessian evaluations. 437 spin.h_count_sim[sim_index] = h_count 438 439 # Warning. 440 spin.warning_sim[sim_index] = warning 441 442 443 # Normal statistics. 444 else: 445 # Chi-squared statistic. 446 spin.chi2 = chi2 447 448 # Iterations. 449 spin.iter = iter_count 450 451 # Function evaluations. 452 spin.f_count = f_count 453 454 # Gradient evaluations. 455 spin.g_count = g_count 456 457 # Hessian evaluations. 458 spin.h_count = h_count 459 460 # Warning. 461 spin.warning = warning 462 463 # Increment the model index. 464 model_index += 1

465 466

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

468 """Deselect spins which have insufficient data to support minimisation. 469 470 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 471 @type data_check: bool 472 @keyword verbose: A flag which if True will allow printouts. 473 @type verbose: bool 474 """ 475 476 # Print out. 477 if verbose: 478 print("\nOver-fit spin deselection:") 479 480 # Checks. 481 check_mol_res_spin_data() 482 check_model_setup() 483 484 # Loop over spin data. 485 deselect_flag = False 486 for spin, spin_id in spin_loop(return_id=True): 487 # Skip deselected spins. 488 if not spin.select: 489 continue 490 491 # Check if data exists. 492 if not hasattr(spin, 'peak_intensity'): 493 warn(RelaxDeselectWarning(spin_id, 'missing intensity data')) 494 spin.select = False 495 deselect_flag = True 496 continue 497 498 # Require 3 or more data points. 499 elif len(spin.peak_intensity) < 3: 500 warn(RelaxDeselectWarning(spin_id, 'insufficient data, 3 or more data points are required')) 501 spin.select = False 502 deselect_flag = True 503 continue 504 505 # Check that the number of relaxation times is complete. 506 if len(spin.peak_intensity) != len(cdp.relax_times): 507 raise RelaxError("The %s peak intensity points of the spin '%s' does not match the expected number of %s (the IDs %s do not match %s)." % (len(spin.peak_intensity), spin_id, len(cdp.relax_times), sorted(spin.peak_intensity.keys()), sorted(cdp.relax_times.keys()))) 508 509 # Final printout. 510 if verbose and not deselect_flag: 511 print("No spins have been deselected.")

512 513

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

515 """Function for returning the peak intensity data structure. 516 517 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 518 @type data_id: str 519 @return: The peak intensity data structure. 520 @rtype: list of float 521 """ 522 523 # Get the spin container. 524 spin = return_spin(data_id) 525 526 # Return the peak intensities. 527 return spin.peak_intensity

528 529

530 - def return_error(self, data_id):

531 """Return the standard deviation data structure. 532 533 @param data_id: The spin identification string, as yielded by the base_data_loop() generator 534 method. 535 @type data_id: str 536 @return: The standard deviation data structure. 537 @rtype: list of float 538 """ 539 540 # Get the spin container. 541 spin = return_spin(data_id) 542 543 # Return the error list. 544 return spin.peak_intensity_err

545 546

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

548 """Pack the Monte Carlo simulation data. 549 550 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 551 @type data_id: str 552 @param sim_data: The Monte Carlo simulation data. 553 @type sim_data: list of float 554 """ 555 556 # Get the spin container. 557 spin = return_spin(data_id) 558 559 # Create the data structure. 560 spin.peak_intensity_sim = sim_data

Source Code for Module specific_analyses.relax_fit.api