Package specific_fns :: Module relax_fit
[hide private]
[frames] | no frames]

Source Code for Module specific_fns.relax_fit

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2004-2013 Edward d'Auvergne                                   # 
  4  #                                                                             # 
  5  # This file is part of the program relax (http://www.nmr-relax.com).          # 
  6  #                                                                             # 
  7  # This program is free software: you can redistribute it and/or modify        # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation, either version 3 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # This program is distributed in the hope that it will be useful,             # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
 19  #                                                                             # 
 20  ############################################################################### 
 21   
 22  # Module docstring. 
 23  """Analysis specific code for the R1 and R2 exponential relaxation curve fitting.""" 
 24   
 25  # Python module imports. 
 26  from minfx.generic import generic_minimise 
 27  from minfx.grid import grid 
 28  from numpy import array, average, dot, float64, identity, zeros 
 29  from numpy.linalg import inv 
 30  from re import match, search 
 31  from warnings import warn 
 32   
 33  # relax module imports. 
 34  from dep_check import C_module_exp_fn 
 35  from generic_fns import pipes 
 36  from generic_fns.mol_res_spin import exists_mol_res_spin_data, generate_spin_id_unique, return_spin, spin_loop 
 37  from relax_errors import RelaxError, RelaxFuncSetupError, RelaxLenError, RelaxNoModelError, RelaxNoSequenceError 
 38  from relax_warnings import RelaxDeselectWarning 
 39  from specific_fns.api_base import API_base 
 40  from specific_fns.api_common import API_common 
 41  from user_functions.data import Uf_tables; uf_tables = Uf_tables() 
 42  from user_functions.objects import Desc_container 
 43   
 44  # C modules. 
 45  if C_module_exp_fn: 
 46      from maths_fns.relax_fit import setup, func, dfunc, d2func, back_calc_I 
 47   
 48   
49 -class Relax_fit(API_base, API_common):
50 """Class containing functions for relaxation curve fitting.""" 51
52 - def __init__(self):
53 """Initialise the class by placing API_common methods into the API.""" 54 55 # Execute the base class __init__ method. 56 super(Relax_fit, self).__init__() 57 58 # Place methods into the API. 59 self.base_data_loop = self._base_data_loop_spin 60 self.model_loop = self._model_loop_spin 61 self.return_conversion_factor = self._return_no_conversion_factor 62 self.return_value = self._return_value_general 63 self.set_error = self._set_error_spin 64 self.set_param_values = self._set_param_values_spin 65 self.set_selected_sim = self._set_selected_sim_spin 66 self.sim_init_values = self._sim_init_values_spin 67 self.sim_return_param = self._sim_return_param_spin 68 self.sim_return_selected = self._sim_return_selected_spin 69 70 # Set up the spin parameters. 71 self.PARAMS.add('intensities', scope='spin', py_type=list, grace_string='\\qPeak intensities\\Q') 72 self.PARAMS.add('relax_times', scope='spin', py_type=list, grace_string='\\qRelaxation time period (s)\\Q') 73 self.PARAMS.add('rx', scope='spin', default=8.0, desc='Either the R1 or R2 relaxation rate', set='params', py_type=float, grace_string='\\qR\\sx\\Q', err=True, sim=True) 74 self.PARAMS.add('i0', scope='spin', default=10000.0, desc='The initial intensity', py_type=float, set='params', grace_string='\\qI\\s0\\Q', err=True, sim=True) 75 self.PARAMS.add('iinf', scope='spin', default=0.0, desc='The intensity at infinity', py_type=float, set='params', grace_string='\\qI\\sinf\\Q', err=True, sim=True) 76 self.PARAMS.add('params', scope='spin', desc='The model parameters', py_type=list) 77 78 # Add the minimisation data. 79 self.PARAMS.add_min_data(min_stats_global=False, min_stats_spin=True)
80 81
82 - def _assemble_param_vector(self, spin=None, sim_index=None):
83 """Assemble the exponential curve parameter vector (as a numpy array). 84 85 @keyword spin: The spin data container. 86 @type spin: SpinContainer instance 87 @keyword sim_index: The optional MC simulation index. 88 @type sim_index: int 89 @return: An array of the parameter values of the exponential model. 90 @rtype: numpy array 91 """ 92 93 # Initialise. 94 param_vector = [] 95 96 # Loop over the model parameters. 97 for i in range(len(spin.params)): 98 # Relaxation rate. 99 if spin.params[i] == 'Rx': 100 if sim_index != None: 101 param_vector.append(spin.rx_sim[sim_index]) 102 elif spin.rx == None: 103 param_vector.append(0.0) 104 else: 105 param_vector.append(spin.rx) 106 107 # Initial intensity. 108 elif spin.params[i] == 'I0': 109 if sim_index != None: 110 param_vector.append(spin.i0_sim[sim_index]) 111 elif spin.i0 == None: 112 param_vector.append(0.0) 113 else: 114 param_vector.append(spin.i0) 115 116 # Intensity at infinity. 117 elif spin.params[i] == 'Iinf': 118 if sim_index != None: 119 param_vector.append(spin.iinf_sim[sim_index]) 120 elif spin.iinf == None: 121 param_vector.append(0.0) 122 else: 123 param_vector.append(spin.iinf) 124 125 # Return a numpy array. 126 return array(param_vector, float64)
127 128
129 - def _assemble_scaling_matrix(self, spin=None, scaling=True):
130 """Create and return the scaling matrix. 131 132 @keyword spin: The spin data container. 133 @type spin: SpinContainer instance 134 @keyword scaling: A flag which if false will cause the identity matrix to be returned. 135 @type scaling: bool 136 @return: The diagonal and square scaling matrix. 137 @rtype: numpy diagonal matrix 138 """ 139 140 # Initialise. 141 scaling_matrix = identity(len(spin.params), float64) 142 i = 0 143 144 # No diagonal scaling. 145 if not scaling: 146 return scaling_matrix 147 148 # Loop over the parameters. 149 for i in range(len(spin.params)): 150 # Relaxation rate. 151 if spin.params[i] == 'Rx': 152 pass 153 154 # Intensity scaling. 155 elif search('^i', spin.params[i]): 156 # Find the position of the first time point. 157 pos = cdp.relax_times.index(min(cdp.relax_times)) 158 159 # Scaling. 160 scaling_matrix[i, i] = 1.0 / average(spin.intensities[pos]) 161 162 # Increment i. 163 i = i + 1 164 165 # Return the scaling matrix. 166 return scaling_matrix
167 168
169 - def _back_calc(self, spin=None, relax_time_id=None):
170 """Back-calculation of peak intensity for the given relaxation time. 171 172 @keyword spin: The spin container. 173 @type spin: SpinContainer instance 174 @keyword relax_time_id: The ID string for the desired relaxation time. 175 @type relax_time_id: str 176 @return: The peak intensity for the desired relaxation time. 177 @rtype: float 178 """ 179 180 # Create the initial parameter vector. 181 param_vector = self._assemble_param_vector(spin=spin) 182 183 # Create a scaling matrix. 184 scaling_matrix = self._assemble_scaling_matrix(spin=spin, scaling=False) 185 186 # The keys. 187 keys = list(spin.intensities.keys()) 188 189 # The peak intensities and times. 190 values = [] 191 errors = [] 192 times = [] 193 for key in keys: 194 values.append(spin.intensities[key]) 195 errors.append(spin.intensity_err[key]) 196 times.append(cdp.relax_times[key]) 197 198 # Initialise the relaxation fit functions. 199 setup(num_params=len(spin.params), num_times=len(cdp.relax_times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_matrix) 200 201 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 202 self._func(param_vector) 203 204 # Get the data back. 205 results = back_calc_I() 206 207 # Return the correct peak height. 208 return results[keys.index(relax_time_id)]
209 210
211 - def _disassemble_param_vector(self, param_vector=None, spin=None, sim_index=None):
212 """Disassemble the parameter vector. 213 214 @keyword param_vector: The parameter vector. 215 @type param_vector: numpy array 216 @keyword spin: The spin data container. 217 @type spin: SpinContainer instance 218 @keyword sim_index: The optional MC simulation index. 219 @type sim_index: int 220 """ 221 222 # Monte Carlo simulations. 223 if sim_index != None: 224 # The relaxation rate. 225 spin.rx_sim[sim_index] = param_vector[0] 226 227 # Initial intensity. 228 spin.i0_sim[sim_index] = param_vector[1] 229 230 # Intensity at infinity. 231 if cdp.curve_type == 'inv': 232 spin.iinf_sim[sim_index] = param_vector[2] 233 234 # Parameter values. 235 else: 236 # The relaxation rate. 237 spin.rx = param_vector[0] 238 239 # Initial intensity. 240 spin.i0 = param_vector[1] 241 242 # Intensity at infinity. 243 if cdp.curve_type == 'inv': 244 spin.iinf = param_vector[2]
245 246
247 - def _func(self, params):
248 """Wrapper function for the C module, for converting numpy arrays. 249 250 @param params: The parameter array from the minimisation code. 251 @type params: numpy array 252 @return: The function value generated by the C module. 253 @rtype: float 254 """ 255 256 # Call the C code. 257 chi2 = func(params.tolist()) 258 259 # Return the chi2 value. 260 return chi2
261 262
263 - def _dfunc(self, params):
264 """Wrapper function for the C module, for converting numpy arrays. 265 266 The currently does nothing. 267 """
268 269
270 - def _d2func(self, params):
271 """Wrapper function for the C module, for converting numpy arrays. 272 273 The currently does nothing. 274 """
275 276
277 - def _grid_search_setup(self, spin=None, param_vector=None, lower=None, upper=None, inc=None, scaling_matrix=None):
278 """The grid search setup function. 279 280 @keyword spin: The spin data container. 281 @type spin: SpinContainer instance 282 @keyword param_vector: The parameter vector. 283 @type param_vector: numpy array 284 @keyword lower: The lower bounds of the grid search which must be equal to the 285 number of parameters in the model. This optional argument is 286 only used when doing a grid search. 287 @type lower: array of numbers 288 @keyword upper: The upper bounds of the grid search which must be equal to the 289 number of parameters in the model. This optional argument is 290 only used when doing a grid search. 291 @type upper: array of numbers 292 @keyword inc: The increments for each dimension of the space for the grid 293 search. The number of elements in the array must equal to the 294 number of parameters in the model. This argument is only used 295 when doing a grid search. 296 @type inc: array of int 297 @keyword scaling_matrix: The scaling matrix. 298 @type scaling_matrix: numpy diagonal matrix 299 @return: A tuple of the grid size and the minimisation options. For the 300 minimisation options, the first dimension corresponds to the 301 model parameter. The second dimension is a list of the number 302 of increments, the lower bound, and upper bound. 303 @rtype: (int, list of lists [int, float, float]) 304 """ 305 306 # The length of the parameter array. 307 n = len(param_vector) 308 309 # Make sure that the length of the parameter array is > 0. 310 if n == 0: 311 raise RelaxError("Cannot run a grid search on a model with zero parameters.") 312 313 # Lower bounds. 314 if lower != None and len(lower) != n: 315 raise RelaxLenError('lower bounds', n) 316 317 # Upper bounds. 318 if upper != None and len(upper) != n: 319 raise RelaxLenError('upper bounds', n) 320 321 # Increments. 322 if isinstance(inc, list) and len(inc) != n: 323 raise RelaxLenError('increment', n) 324 elif isinstance(inc, int): 325 inc = [inc]*n 326 327 # Set up the default bounds. 328 if not lower: 329 # Init. 330 lower = [] 331 upper = [] 332 333 # Loop over the parameters. 334 for i in range(n): 335 # Relaxation rate (from 0 to 20 s^-1). 336 if spin.params[i] == 'Rx': 337 lower.append(0.0) 338 upper.append(20.0) 339 340 # Intensity 341 elif search('^I', spin.params[i]): 342 # Find the ID of the first time point. 343 min_time = min(cdp.relax_times.values()) 344 for key in list(cdp.relax_times.keys()): 345 if cdp.relax_times[key] == min_time: 346 id = key 347 break 348 349 # Defaults. 350 lower.append(0.0) 351 upper.append(average(spin.intensities[id])) 352 353 # Parameter scaling. 354 for i in range(n): 355 lower[i] = lower[i] / scaling_matrix[i, i] 356 upper[i] = upper[i] / scaling_matrix[i, i] 357 358 return inc, lower, upper
359 360
361 - def _linear_constraints(self, spin=None, scaling_matrix=None):
362 """Set up the relaxation curve fitting linear constraint matrices A and b. 363 364 Standard notation 365 ================= 366 367 The relaxation rate constraints are:: 368 369 Rx >= 0 370 371 The intensity constraints are:: 372 373 I0 >= 0 374 Iinf >= 0 375 376 377 Matrix notation 378 =============== 379 380 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter 381 values, and b is a vector of scalars, these inequality constraints are:: 382 383 | 1 0 0 | | Rx | | 0 | 384 | | | | | | 385 | 1 0 0 | . | I0 | >= | 0 | 386 | | | | | | 387 | 1 0 0 | | Iinf | | 0 | 388 389 390 @keyword spin: The spin data container. 391 @type spin: SpinContainer instance 392 @keyword scaling_matrix: The diagonal, square scaling matrix. 393 @type scaling_matrix: numpy diagonal matrix 394 """ 395 396 # Initialisation (0..j..m). 397 A = [] 398 b = [] 399 n = len(spin.params) 400 zero_array = zeros(n, float64) 401 i = 0 402 j = 0 403 404 # Loop over the parameters. 405 for k in range(len(spin.params)): 406 # Relaxation rate. 407 if spin.params[k] == 'Rx': 408 # Rx >= 0. 409 A.append(zero_array * 0.0) 410 A[j][i] = 1.0 411 b.append(0.0) 412 j = j + 1 413 414 # Intensity parameter. 415 elif search('^I', spin.params[k]): 416 # I0, Iinf >= 0. 417 A.append(zero_array * 0.0) 418 A[j][i] = 1.0 419 b.append(0.0) 420 j = j + 1 421 422 # Increment i. 423 i = i + 1 424 425 # Convert to numpy data structures. 426 A = array(A, float64) 427 b = array(b, float64) 428 429 return A, b
430 431
432 - def _model_setup(self, model, params):
433 """Update various model specific data structures. 434 435 @param model: The exponential curve type. 436 @type model: str 437 @param params: A list consisting of the model parameters. 438 @type params: list of str 439 """ 440 441 # Set the model. 442 cdp.curve_type = model 443 444 # Loop over the sequence. 445 for spin in spin_loop(): 446 # Skip deselected spins. 447 if not spin.select: 448 continue 449 450 # Initialise the data structures (if needed). 451 self.data_init(spin) 452 453 # The model and parameter names. 454 spin.model = model 455 spin.params = params
456 457
458 - def _relax_time(self, time=0.0, spectrum_id=None):
459 """Set the relaxation time period associated with a given spectrum. 460 461 @keyword time: The time, in seconds, of the relaxation period. 462 @type time: float 463 @keyword spectrum_id: The spectrum identification string. 464 @type spectrum_id: str 465 """ 466 467 # Test if the spectrum id exists. 468 if spectrum_id not in cdp.spectrum_ids: 469 raise RelaxError("The peak heights corresponding to spectrum id '%s' have not been loaded." % spectrum_id) 470 471 # Initialise the global relaxation time data structure if needed. 472 if not hasattr(cdp, 'relax_times'): 473 cdp.relax_times = {} 474 475 # Add the time at the correct position. 476 cdp.relax_times[spectrum_id] = time
477 478
479 - def _select_model(self, model='exp'):
480 """Function for selecting the model of the exponential curve. 481 482 @keyword model: The exponential curve type. Can be one of 'exp' or 'inv'. 483 @type model: str 484 """ 485 486 # Test if the current pipe exists. 487 pipes.test() 488 489 # Test if the pipe type is set to 'relax_fit'. 490 function_type = cdp.pipe_type 491 if function_type != 'relax_fit': 492 raise RelaxFuncSetupError(specific_setup.get_string(function_type)) 493 494 # Test if sequence data is loaded. 495 if not exists_mol_res_spin_data(): 496 raise RelaxNoSequenceError 497 498 # Two parameter exponential fit. 499 if model == 'exp': 500 print("Two parameter exponential fit.") 501 params = ['Rx', 'I0'] 502 503 # Three parameter inversion recovery fit. 504 elif model == 'inv': 505 print("Three parameter inversion recovery fit.") 506 params = ['Rx', 'I0', 'Iinf'] 507 508 # Invalid model. 509 else: 510 raise RelaxError("The model '" + model + "' is invalid.") 511 512 # Set up the model. 513 self._model_setup(model, params)
514 515
516 - def create_mc_data(self, data_id=None):
517 """Create the Monte Carlo peak intensity data. 518 519 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 520 @type data_id: str 521 @return: The Monte Carlo simulation data. 522 @rtype: list of floats 523 """ 524 525 # Initialise the MC data data structure. 526 mc_data = {} 527 528 # Get the spin container. 529 spin = return_spin(data_id) 530 531 # Skip deselected spins. 532 if not spin.select: 533 return 534 535 # Skip spins which have no data. 536 if not hasattr(spin, 'intensities'): 537 return 538 539 # Test if the model is set. 540 if not hasattr(spin, 'model') or not spin.model: 541 raise RelaxNoModelError 542 543 # Loop over the spectral time points. 544 for id in list(cdp.relax_times.keys()): 545 # Back calculate the value. 546 value = self._back_calc(spin=spin, relax_time_id=id) 547 548 # Append the value. 549 mc_data[id] = value 550 551 # Return the MC data. 552 return mc_data
553 554
555 - def data_init(self, data_cont, sim=False):
556 """Initialise the spin specific data structures. 557 558 @param data_cont: The spin container. 559 @type data_cont: SpinContainer instance 560 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure. 561 @type sim: bool 562 """ 563 564 # Loop over the data structure names. 565 for name in self.data_names(set='params'): 566 # Data structures which are initially empty arrays. 567 list_data = [ 'params' ] 568 if name in list_data: 569 init_data = [] 570 571 # Otherwise initialise the data structure to None. 572 else: 573 init_data = None 574 575 # If the name is not in 'data_cont', add it. 576 if not hasattr(data_cont, name): 577 setattr(data_cont, name, init_data)
578 579 580 default_value_doc = Desc_container("Relaxation curve fitting default values") 581 default_value_doc.add_paragraph("These values are completely arbitrary as peak heights (or volumes) are extremely variable and the Rx value is a compensation for both the R1 and R2 values.") 582 _table = uf_tables.add_table(label="table: curve-fit default values", caption="Relaxation curve fitting default values.") 583 _table.add_headings(["Data type", "Object name", "Value"]) 584 _table.add_row(["Relaxation rate", "'rx'", "8.0"]) 585 _table.add_row(["Initial intensity", "'i0'", "10000.0"]) 586 _table.add_row(["Intensity at infinity", "'iinf'", "0.0"]) 587 default_value_doc.add_table(_table.label) 588 589
590 - def grid_search(self, lower=None, upper=None, inc=None, constraints=True, verbosity=1, sim_index=None):
591 """The exponential curve fitting grid search method. 592 593 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 594 @type lower: array of numbers 595 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 596 @type upper: array of numbers 597 @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. 598 @type inc: array of int 599 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 600 @type constraints: bool 601 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 602 @type verbosity: int 603 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 604 @type sim_index: int 605 """ 606 607 # Minimisation. 608 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
609 610
611 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling=True, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
612 """Relaxation curve fitting minimisation method. 613 614 @keyword min_algor: The minimisation algorithm to use. 615 @type min_algor: str 616 @keyword min_options: An array of options to be used by the minimisation algorithm. 617 @type min_options: array of str 618 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 619 @type func_tol: None or float 620 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 621 @type grad_tol: None or float 622 @keyword max_iterations: The maximum number of iterations for the algorithm. 623 @type max_iterations: int 624 @keyword constraints: If True, constraints are used during optimisation. 625 @type constraints: bool 626 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 627 @type scaling: bool 628 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 629 @type verbosity: int 630 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 631 @type sim_index: None or int 632 @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. 633 @type lower: array of numbers 634 @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. 635 @type upper: array of numbers 636 @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. 637 @type inc: array of int 638 """ 639 640 # Test if sequence data is loaded. 641 if not exists_mol_res_spin_data(): 642 raise RelaxNoSequenceError 643 644 # Loop over the sequence. 645 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 646 # Skip deselected spins. 647 if not spin.select: 648 continue 649 650 # Skip spins which have no data. 651 if not hasattr(spin, 'intensities'): 652 continue 653 654 # Create the initial parameter vector. 655 param_vector = self._assemble_param_vector(spin=spin) 656 657 # Diagonal scaling. 658 scaling_matrix = self._assemble_scaling_matrix(spin=spin, scaling=scaling) 659 if len(scaling_matrix): 660 param_vector = dot(inv(scaling_matrix), param_vector) 661 662 # Get the grid search minimisation options. 663 if match('^[Gg]rid', min_algor): 664 inc, lower, upper = self._grid_search_setup(spin=spin, param_vector=param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix) 665 666 # Linear constraints. 667 if constraints: 668 A, b = self._linear_constraints(spin=spin, scaling_matrix=scaling_matrix) 669 else: 670 A, b = None, None 671 672 # Print out. 673 if verbosity >= 1: 674 # Get the spin id string. 675 spin_id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=spin.num, spin_name=spin.name) 676 677 # Individual spin printout. 678 if verbosity >= 2: 679 print("\n\n") 680 681 string = "Fitting to spin " + repr(spin_id) 682 print("\n\n" + string) 683 print(len(string) * '~') 684 685 686 # Initialise the function to minimise. 687 ###################################### 688 689 # The keys. 690 keys = list(spin.intensities.keys()) 691 692 # The peak intensities and times. 693 values = [] 694 errors = [] 695 times = [] 696 for key in keys: 697 # The values. 698 if sim_index == None: 699 values.append(spin.intensities[key]) 700 else: 701 values.append(spin.sim_intensities[sim_index][key]) 702 703 # The errors. 704 errors.append(spin.intensity_err[key]) 705 706 # The relaxation times. 707 times.append(cdp.relax_times[key]) 708 709 setup(num_params=len(spin.params), num_times=len(values), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_matrix.tolist()) 710 711 712 # Setup the minimisation algorithm when constraints are present. 713 ################################################################ 714 715 if constraints and not match('^[Gg]rid', min_algor): 716 algor = min_options[0] 717 else: 718 algor = min_algor 719 720 721 # Levenberg-Marquardt minimisation. 722 ################################### 723 724 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 725 # Reconstruct the error data structure. 726 lm_error = zeros(len(spin.relax_times), float64) 727 index = 0 728 for k in range(len(spin.relax_times)): 729 lm_error[index:index+len(relax_error[k])] = relax_error[k] 730 index = index + len(relax_error[k]) 731 732 min_options = min_options + (self.relax_fit.lm_dri, lm_error) 733 734 735 # Minimisation. 736 ############### 737 738 # Grid search. 739 if search('^[Gg]rid', min_algor): 740 results = grid(func=self._func, args=(), num_incs=inc, lower=lower, upper=upper, A=A, b=b, verbosity=verbosity) 741 742 # Unpack the results. 743 param_vector, chi2, iter_count, warning = results 744 f_count = iter_count 745 g_count = 0.0 746 h_count = 0.0 747 748 # Minimisation. 749 else: 750 results = generic_minimise(func=self._func, dfunc=self._dfunc, d2func=self._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) 751 752 # Unpack the results. 753 if results == None: 754 return 755 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 756 757 # Scaling. 758 if scaling: 759 param_vector = dot(scaling_matrix, param_vector) 760 761 # Disassemble the parameter vector. 762 self._disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index) 763 764 # Monte Carlo minimisation statistics. 765 if sim_index != None: 766 # Chi-squared statistic. 767 spin.chi2_sim[sim_index] = chi2 768 769 # Iterations. 770 spin.iter_sim[sim_index] = iter_count 771 772 # Function evaluations. 773 spin.f_count_sim[sim_index] = f_count 774 775 # Gradient evaluations. 776 spin.g_count_sim[sim_index] = g_count 777 778 # Hessian evaluations. 779 spin.h_count_sim[sim_index] = h_count 780 781 # Warning. 782 spin.warning_sim[sim_index] = warning 783 784 785 # Normal statistics. 786 else: 787 # Chi-squared statistic. 788 spin.chi2 = chi2 789 790 # Iterations. 791 spin.iter = iter_count 792 793 # Function evaluations. 794 spin.f_count = f_count 795 796 # Gradient evaluations. 797 spin.g_count = g_count 798 799 # Hessian evaluations. 800 spin.h_count = h_count 801 802 # Warning. 803 spin.warning = warning
804 805
806 - def overfit_deselect(self, data_check=True, verbose=True):
807 """Deselect spins which have insufficient data to support minimisation. 808 809 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 810 @type data_check: bool 811 @keyword verbose: A flag which if True will allow printouts. 812 @type verbose: bool 813 """ 814 815 # Print out. 816 if verbose: 817 print("\nOver-fit spin deselection:") 818 819 # Test the sequence data exists. 820 if not exists_mol_res_spin_data(): 821 raise RelaxNoSequenceError 822 823 # Loop over spin data. 824 deselect_flag = False 825 for spin, spin_id in spin_loop(return_id=True): 826 # Skip deselected spins. 827 if not spin.select: 828 continue 829 830 # Check if data exists. 831 if not hasattr(spin, 'intensities'): 832 warn(RelaxDeselectWarning(spin_id, 'missing intensity data')) 833 spin.select = False 834 deselect_flag = True 835 continue 836 837 # Require 3 or more data points. 838 elif len(spin.intensities) < 3: 839 warn(RelaxDeselectWarning(spin_id, 'insufficient data, 3 or more data points are required')) 840 spin.select = False 841 deselect_flag = True 842 continue 843 844 # Check that the number of relaxation times is complete. 845 if len(spin.intensities) != len(cdp.relax_times): 846 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.intensities), spin_id, len(cdp.relax_times), list(spin.intensities.keys()), list(cdp.relax_times.keys()))) 847 848 # Final printout. 849 if verbose and not deselect_flag: 850 print("No spins have been deselected.")
851 852
853 - def return_data(self, data_id=None):
854 """Function for returning the peak intensity data structure. 855 856 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 857 @type data_id: str 858 @return: The peak intensity data structure. 859 @rtype: list of float 860 """ 861 862 # Get the spin container. 863 spin = return_spin(data_id) 864 865 # Return the peak intensities. 866 return spin.intensities
867 868 869 return_data_name_doc = Desc_container("Relaxation curve fitting data type string matching patterns") 870 _table = uf_tables.add_table(label="table: curve-fit data type patterns", caption="Relaxation curve fitting data type string matching patterns.") 871 _table.add_headings(["Data type", "Object name"]) 872 _table.add_row(["Relaxation rate", "'rx'"]) 873 _table.add_row(["Peak intensities (series)", "'intensities'"]) 874 _table.add_row(["Initial intensity", "'i0'"]) 875 _table.add_row(["Intensity at infinity", "'iinf'"]) 876 _table.add_row(["Relaxation period times (series)", "'relax_times'"]) 877 return_data_name_doc.add_table(_table.label) 878 879
880 - def return_error(self, data_id):
881 """Return the standard deviation data structure. 882 883 @param data_id: The spin identification string, as yielded by the base_data_loop() generator 884 method. 885 @type data_id: str 886 @return: The standard deviation data structure. 887 @rtype: list of float 888 """ 889 890 # Get the spin container. 891 spin = return_spin(data_id) 892 893 # Return the error list. 894 return spin.intensity_err
895 896
897 - def return_units(self, param):
898 """Dummy method which returns None as the stats have no units. 899 900 @param param: The name of the parameter to return the units string for. 901 @type param: str 902 @return: Nothing. 903 @rtype: None 904 """ 905 906 # Unitless. 907 return None
908 909 910 set_doc = Desc_container("Relaxation curve fitting set details") 911 set_doc.add_paragraph("Only three parameters can be set, the relaxation rate (Rx), the initial intensity (I0), and the intensity at infinity (Iinf). Setting the parameter Iinf has no effect if the chosen model is that of the exponential curve which decays to zero.") 912 913
914 - def sim_pack_data(self, data_id, sim_data):
915 """Pack the Monte Carlo simulation data. 916 917 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 918 @type data_id: str 919 @param sim_data: The Monte Carlo simulation data. 920 @type sim_data: list of float 921 """ 922 923 # Get the spin container. 924 spin = return_spin(data_id) 925 926 # Create the data structure. 927 spin.sim_intensities = sim_data
928