Package specific_analyses :: Package relax_fit :: Module api
[hide private]
[frames] | no frames]

Source Code for Module specific_analyses.relax_fit.api

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003-2014 Edward d'Auvergne                                   # 
  4  # Copyright (C) 2006 Chris MacRaild                                           # 
  5  # Copyright (C) 2008-2009 Sebastien Morin                                     # 
  6  #                                                                             # 
  7  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  8  #                                                                             # 
  9  # This program is free software: you can redistribute it and/or modify        # 
 10  # it under the terms of the GNU General Public License as published by        # 
 11  # the Free Software Foundation, either version 3 of the License, or           # 
 12  # (at your option) any later version.                                         # 
 13  #                                                                             # 
 14  # This program is distributed in the hope that it will be useful,             # 
 15  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 16  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 17  # GNU General Public License for more details.                                # 
 18  #                                                                             # 
 19  # You should have received a copy of the GNU General Public License           # 
 20  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 21  #                                                                             # 
 22  ############################################################################### 
 23   
 24  # Module docstring. 
 25  """The R1 and R2 exponential relaxation curve fitting API object.""" 
 26   
 27  # Python module imports. 
 28  from minfx.generic import generic_minimise 
 29  from minfx.grid import grid 
 30  from numpy import asarray, dot, float64, transpose, zeros 
 31  from numpy.linalg import inv 
 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.errors import RelaxError, RelaxNoModelError 
 39  from lib.text.sectioning import subsection 
 40  from lib.warnings import RelaxDeselectWarning 
 41  from pipe_control.mol_res_spin import check_mol_res_spin_data, return_spin, spin_loop 
 42  from specific_analyses.api_base import API_base 
 43  from specific_analyses.api_common import API_common 
 44  from specific_analyses.relax_fit.checks import check_model_setup 
 45  from specific_analyses.relax_fit.optimisation import back_calc 
 46  from specific_analyses.relax_fit.parameter_object import Relax_fit_params 
 47  from specific_analyses.relax_fit.parameters import assemble_param_vector, disassemble_param_vector, linear_constraints 
 48  from target_functions.relax_fit_wrapper import Relax_fit_opt 
 49   
 50   
51 -class Relax_fit(API_base, API_common):
52 """Class containing functions for relaxation curve fitting.""" 53 54 # Class variable for storing the class instance (for the singleton design pattern). 55 instance = None 56
57 - def __init__(self):
58 """Initialise the class by placing API_common methods into the API.""" 59 60 # Place methods into the API. 61 self.base_data_loop = self._base_data_loop_spin 62 self.model_loop = self._model_loop_spin 63 self.print_model_title = self._print_model_title_spin 64 self.return_conversion_factor = self._return_no_conversion_factor 65 self.return_value = self._return_value_general 66 self.set_error = self._set_error_spin 67 self.set_param_values = self._set_param_values_spin 68 self.set_selected_sim = self._set_selected_sim_spin 69 self.sim_init_values = self._sim_init_values_spin 70 self.sim_return_param = self._sim_return_param_spin 71 self.sim_return_selected = self._sim_return_selected_spin 72 73 # Place a copy of the parameter list object in the instance namespace. 74 self._PARAMS = Relax_fit_params()
75 76
77 - def covariance_matrix(self, model_info=None, verbosity=1):
78 """Return the Jacobian and weights required for parameter errors via the covariance matrix. 79 80 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method. 81 @type model_info: SpinContainer instance, str 82 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 83 @type verbosity: int 84 @return: The Jacobian and weight matrices for the given model. 85 @rtype: numpy rank-2 array, numpy rank-2 array 86 """ 87 88 # Unpack the data. 89 spin, spin_id = model_info 90 91 # Check that the C modules have been compiled. 92 if not C_module_exp_fn: 93 raise RelaxError("Relaxation curve fitting is not available. Try compiling the C modules on your platform.") 94 95 # Raise Error, if not optimised. 96 if not (hasattr(spin, 'rx') and hasattr(spin, 'i0')): 97 raise RelaxError("Spin '%s' does not contain optimised 'rx' and 'i0' values. Try execute: minimise.execute(min_algor='Newton', constraints=False)"%(spin_id)) 98 99 # Raise warning, if gradient count is 0. This could point to a lack of minimisation first. 100 if hasattr(spin, 'g_count'): 101 if getattr(spin, 'g_count') == 0.0: 102 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) 103 warn(RelaxWarning("%s." % text)) 104 105 # Print information. 106 if verbosity >= 1: 107 # Individual spin block section. 108 top = 2 109 if verbosity >= 2: 110 top += 2 111 subsection(file=sys.stdout, text="Estimating rx error for spin: %s"%spin_id, prespace=top) 112 113 # The peak intensities and times. 114 values = [] 115 errors = [] 116 times = [] 117 for key in spin.peak_intensity: 118 values.append(spin.peak_intensity[key]) 119 errors.append(spin.peak_intensity_err[key]) 120 times.append(cdp.relax_times[key]) 121 122 # Convert to numpy array. 123 values = asarray(values) 124 errors = asarray(errors) 125 times = asarray(times) 126 127 # Create the parameter vector and scaling matrix (as a diagonalised list). 128 param_vector = assemble_param_vector(spin=spin) 129 scaling_list = [] 130 for i in range(len(spin.params)): 131 scaling_list.append(1.0) 132 133 # Initialise data in C code. 134 model = Relax_fit_opt(model=spin.model, num_params=len(param_vector), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 135 136 # Use the direct Jacobian from function. 137 jacobian_matrix_exp = transpose(asarray( model.jacobian(param_vector) ) ) 138 weights = 1. / errors**2 139 140 # Return the matrices. 141 return jacobian_matrix_exp, weights
142 143
144 - def create_mc_data(self, data_id=None):
145 """Create the Monte Carlo peak intensity data. 146 147 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 148 @type data_id: str 149 @return: The Monte Carlo simulation data. 150 @rtype: list of floats 151 """ 152 153 # Initialise the MC data data structure. 154 mc_data = {} 155 156 # Get the spin container. 157 spin = return_spin(spin_id=data_id) 158 159 # Skip deselected spins. 160 if not spin.select: 161 return 162 163 # Skip spins which have no data. 164 if not hasattr(spin, 'peak_intensity'): 165 return 166 167 # Test if the model is set. 168 if not hasattr(spin, 'model') or not spin.model: 169 raise RelaxNoModelError 170 171 # Loop over the spectral time points. 172 for id in cdp.relax_times: 173 # Back calculate the value. 174 value = back_calc(spin=spin, relax_time_id=id) 175 176 # Append the value. 177 mc_data[id] = value 178 179 # Return the MC data. 180 return mc_data
181 182
183 - def data_init(self, data, sim=False):
184 """Initialise the spin specific data structures. 185 186 @param data: The spin ID string from the _base_data_loop_spin() method. 187 @type data: str 188 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure. 189 @type sim: bool 190 """ 191 192 # Get the spin container. 193 spin = return_spin(spin_id=data) 194 195 # Loop over the data structure names. 196 for name in self.data_names(set='params'): 197 # Data structures which are initially empty arrays. 198 list_data = [ 'params' ] 199 if name in list_data: 200 init_data = [] 201 202 # Otherwise initialise the data structure to None. 203 else: 204 init_data = None 205 206 # If the name is not in the spin container, add it. 207 if not hasattr(spin, name): 208 setattr(spin, name, init_data)
209 210
211 - def get_param_names(self, model_info=None):
212 """Return a vector of parameter names. 213 214 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method. 215 @type model_info: SpinContainer instance, str 216 @return: The vector of parameter names. 217 @rtype: list of str 218 """ 219 220 # Check that the model is setup. 221 check_model_setup() 222 223 # Unpack the data. 224 spin, spin_id = model_info 225 226 # Return the parameter names. 227 return spin.params
228 229
230 - def get_param_values(self, model_info=None, sim_index=None):
231 """Return a vector of parameter values. 232 233 @keyword model_info: The spin container and the spin ID string from the _model_loop_spin() method. 234 @type model_info: SpinContainer instance, str 235 @keyword sim_index: The optional Monte Carlo simulation index. 236 @type sim_index: int 237 @return: The vector of parameter values. 238 @rtype: list of str 239 """ 240 241 # Unpack the data. 242 spin, spin_id = model_info 243 244 # Return the vector. 245 return assemble_param_vector(spin=spin, sim_index=sim_index)
246 247
248 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
249 """The exponential curve fitting grid search method. 250 251 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model. 252 @type lower: list of lists of numbers 253 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model. 254 @type upper: list of lists of numbers 255 @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. 256 @type inc: list of lists of int 257 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 258 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 259 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 260 @type constraints: bool 261 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 262 @type verbosity: int 263 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 264 @type sim_index: int 265 """ 266 267 # Minimisation. 268 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
269 270
271 - 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):
272 """Relaxation curve fitting minimisation method. 273 274 @keyword min_algor: The minimisation algorithm to use. 275 @type min_algor: str 276 @keyword min_options: An array of options to be used by the minimisation algorithm. 277 @type min_options: array of str 278 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 279 @type func_tol: None or float 280 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 281 @type grad_tol: None or float 282 @keyword max_iterations: The maximum number of iterations for the algorithm. 283 @type max_iterations: int 284 @keyword constraints: If True, constraints are used during optimisation. 285 @type constraints: bool 286 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 287 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 288 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 289 @type verbosity: int 290 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 291 @type sim_index: None or int 292 @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. 293 @type lower: list of lists of numbers 294 @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. 295 @type upper: list of lists of numbers 296 @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. 297 @type inc: list of lists of int 298 """ 299 300 # Checks. 301 check_mol_res_spin_data() 302 303 # Loop over the sequence. 304 model_index = 0 305 for spin, spin_id in self.model_loop(): 306 # Skip deselected spins. 307 if not spin.select: 308 continue 309 310 # Skip spins which have no data. 311 if not hasattr(spin, 'peak_intensity'): 312 continue 313 314 # Create the initial parameter vector. 315 param_vector = assemble_param_vector(spin=spin) 316 317 # Diagonal scaling. 318 if scaling_matrix[model_index] is not None: 319 param_vector = dot(inv(scaling_matrix[model_index]), param_vector) 320 321 # Linear constraints. 322 if constraints: 323 A, b = linear_constraints(spin=spin, scaling_matrix=scaling_matrix[model_index]) 324 else: 325 A, b = None, None 326 327 # Print out. 328 if verbosity >= 1: 329 # Individual spin printout. 330 if verbosity >= 2: 331 print("\n\n") 332 333 string = "Fitting to spin " + repr(spin_id) 334 print("\n\n" + string) 335 print(len(string) * '~') 336 337 338 # Initialise the function to minimise. 339 ###################################### 340 341 # The peak intensities and times. 342 values = [] 343 errors = [] 344 times = [] 345 for key in spin.peak_intensity: 346 # The values. 347 if sim_index == None: 348 values.append(spin.peak_intensity[key]) 349 else: 350 values.append(spin.peak_intensity_sim[sim_index][key]) 351 352 # The errors. 353 errors.append(spin.peak_intensity_err[key]) 354 355 # The relaxation times. 356 times.append(cdp.relax_times[key]) 357 358 # The scaling matrix in a diagonalised list form. 359 scaling_list = [] 360 if scaling_matrix[model_index] is None: 361 for i in range(len(param_vector)): 362 scaling_list.append(1.0) 363 else: 364 for i in range(len(scaling_matrix[model_index])): 365 scaling_list.append(scaling_matrix[model_index][i, i]) 366 367 # Set up the target function. 368 model = Relax_fit_opt(model=spin.model, num_params=len(spin.params), values=values, errors=errors, relax_times=times, scaling_matrix=scaling_list) 369 370 371 # Setup the minimisation algorithm when constraints are present. 372 ################################################################ 373 374 if constraints and not match('^[Gg]rid', min_algor): 375 algor = min_options[0] 376 else: 377 algor = min_algor 378 379 380 # Levenberg-Marquardt minimisation. 381 ################################### 382 383 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 384 # Reconstruct the error data structure. 385 lm_error = zeros(len(spin.relax_times), float64) 386 index = 0 387 for k in range(len(spin.relax_times)): 388 lm_error[index:index+len(relax_error[k])] = relax_error[k] 389 index = index + len(relax_error[k]) 390 391 min_options = min_options + (self.relax_fit.lm_dri, lm_error) 392 393 394 # Minimisation. 395 ############### 396 397 # Grid search. 398 if search('^[Gg]rid', min_algor): 399 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) 400 401 # Unpack the results. 402 param_vector, chi2, iter_count, warning = results 403 f_count = iter_count 404 g_count = 0.0 405 h_count = 0.0 406 407 # Minimisation. 408 else: 409 results = generic_minimise(func=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) 410 411 # Unpack the results. 412 if results == None: 413 return 414 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 415 416 # Scaling. 417 if scaling_matrix[model_index] is not None: 418 param_vector = dot(scaling_matrix[model_index], param_vector) 419 420 # Disassemble the parameter vector. 421 disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index) 422 423 # Monte Carlo minimisation statistics. 424 if sim_index != None: 425 # Chi-squared statistic. 426 spin.chi2_sim[sim_index] = chi2 427 428 # Iterations. 429 spin.iter_sim[sim_index] = iter_count 430 431 # Function evaluations. 432 spin.f_count_sim[sim_index] = f_count 433 434 # Gradient evaluations. 435 spin.g_count_sim[sim_index] = g_count 436 437 # Hessian evaluations. 438 spin.h_count_sim[sim_index] = h_count 439 440 # Warning. 441 spin.warning_sim[sim_index] = warning 442 443 444 # Normal statistics. 445 else: 446 # Chi-squared statistic. 447 spin.chi2 = chi2 448 449 # Iterations. 450 spin.iter = iter_count 451 452 # Function evaluations. 453 spin.f_count = f_count 454 455 # Gradient evaluations. 456 spin.g_count = g_count 457 458 # Hessian evaluations. 459 spin.h_count = h_count 460 461 # Warning. 462 spin.warning = warning 463 464 # Increment the model index. 465 model_index += 1
466 467
468 - def overfit_deselect(self, data_check=True, verbose=True):
469 """Deselect spins which have insufficient data to support minimisation. 470 471 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 472 @type data_check: bool 473 @keyword verbose: A flag which if True will allow printouts. 474 @type verbose: bool 475 """ 476 477 # Print out. 478 if verbose: 479 print("\nOver-fit spin deselection:") 480 481 # Checks. 482 check_mol_res_spin_data() 483 check_model_setup() 484 485 # Loop over spin data. 486 deselect_flag = False 487 for spin, spin_id in spin_loop(return_id=True): 488 # Skip deselected spins. 489 if not spin.select: 490 continue 491 492 # Check if data exists. 493 if not hasattr(spin, 'peak_intensity'): 494 warn(RelaxDeselectWarning(spin_id, 'missing intensity data')) 495 spin.select = False 496 deselect_flag = True 497 continue 498 499 # Require 3 or more data points. 500 elif len(spin.peak_intensity) < 3: 501 warn(RelaxDeselectWarning(spin_id, 'insufficient data, 3 or more data points are required')) 502 spin.select = False 503 deselect_flag = True 504 continue 505 506 # Check that the number of relaxation times is complete. 507 if len(spin.peak_intensity) != len(cdp.relax_times): 508 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()))) 509 510 # Final printout. 511 if verbose and not deselect_flag: 512 print("No spins have been deselected.")
513 514
515 - def return_data(self, data_id=None):
516 """Function for returning the peak intensity data structure. 517 518 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 519 @type data_id: str 520 @return: The peak intensity data structure. 521 @rtype: list of float 522 """ 523 524 # Get the spin container. 525 spin = return_spin(spin_id=data_id) 526 527 # Return the peak intensities. 528 return spin.peak_intensity
529 530
531 - def return_error(self, data_id):
532 """Return the standard deviation data structure. 533 534 @param data_id: The spin identification string, as yielded by the base_data_loop() generator 535 method. 536 @type data_id: str 537 @return: The standard deviation data structure. 538 @rtype: list of float 539 """ 540 541 # Get the spin container. 542 spin = return_spin(spin_id=data_id) 543 544 # Return the error list. 545 return spin.peak_intensity_err
546 547
548 - def sim_pack_data(self, data_id, sim_data):
549 """Pack the Monte Carlo simulation data. 550 551 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 552 @type data_id: str 553 @param sim_data: The Monte Carlo simulation data. 554 @type sim_data: list of float 555 """ 556 557 # Get the spin container. 558 spin = return_spin(spin_id=data_id) 559 560 # Create the data structure. 561 spin.peak_intensity_sim = sim_data
562