specific_fns.relax

49 """Class containing functions for relaxation curve fitting.""" 50

51 - def __init__(self):

52 """Initialise the class by placing API_common methods into the API.""" 53 54 # Place methods into the API. 55 self.base_data_loop = self._base_data_loop_spin 56 self.default_value = self._default_value_spin 57 self.model_loop = self._model_loop_spin 58 self.return_conversion_factor = self._return_no_conversion_factor 59 self.return_data_name = self._return_data_name_spin 60 self.return_grace_string = self._return_grace_string_spin 61 self.return_value = self._return_value_general 62 self.set_error = self._set_error_spin 63 self.set_param_values = self._set_param_values_spin 64 self.set_selected_sim = self._set_selected_sim_spin 65 self.sim_init_values = self._sim_init_values_spin 66 self.sim_return_param = self._sim_return_param_spin 67 self.sim_return_selected = self._sim_return_selected_spin 68 69 # Set up the spin parameters. 70 self.SPIN_PARAMS.add('rx', default=8.0, grace_string='\\qR\\sx\\Q') 71 self.SPIN_PARAMS.add('intensities', grace_string='\\qPeak intensities\\Q') 72 self.SPIN_PARAMS.add('i0', default=10000.0, grace_string='\\qI\\s0\\Q') 73 self.SPIN_PARAMS.add('iinf', default=0.0, grace_string='\\qI\\sinf\\Q') 74 self.SPIN_PARAMS.add('relax_times', grace_string='\\qRelaxation time period (s)\\Q')

75 76

77 - def _assemble_param_vector(self, spin=None, sim_index=None):

78 """Assemble the exponential curve parameter vector (as a numpy array). 79 80 @keyword spin: The spin data container. 81 @type spin: SpinContainer instance 82 @keyword sim_index: The optional MC simulation index. 83 @type sim_index: int 84 @return: An array of the parameter values of the exponential model. 85 @rtype: numpy array 86 """ 87 88 # Initialise. 89 param_vector = [] 90 91 # Loop over the model parameters. 92 for i in xrange(len(spin.params)): 93 # Relaxation rate. 94 if spin.params[i] == 'Rx': 95 if sim_index != None: 96 param_vector.append(spin.rx_sim[sim_index]) 97 elif spin.rx == None: 98 param_vector.append(0.0) 99 else: 100 param_vector.append(spin.rx) 101 102 # Initial intensity. 103 elif spin.params[i] == 'I0': 104 if sim_index != None: 105 param_vector.append(spin.i0_sim[sim_index]) 106 elif spin.i0 == None: 107 param_vector.append(0.0) 108 else: 109 param_vector.append(spin.i0) 110 111 # Intensity at infinity. 112 elif spin.params[i] == 'Iinf': 113 if sim_index != None: 114 param_vector.append(spin.iinf_sim[sim_index]) 115 elif spin.iinf == None: 116 param_vector.append(0.0) 117 else: 118 param_vector.append(spin.iinf) 119 120 # Return a numpy array. 121 return array(param_vector, float64)

122 123

124 - def _assemble_scaling_matrix(self, spin=None, scaling=True):

125 """Create and return the scaling matrix. 126 127 @keyword spin: The spin data container. 128 @type spin: SpinContainer instance 129 @keyword scaling: A flag which if false will cause the identity matrix to be returned. 130 @type scaling: bool 131 @return: The diagonal and square scaling matrix. 132 @rtype: numpy diagonal matrix 133 """ 134 135 # Initialise. 136 scaling_matrix = identity(len(spin.params), float64) 137 i = 0 138 139 # No diagonal scaling. 140 if not scaling: 141 return scaling_matrix 142 143 # Loop over the parameters. 144 for i in xrange(len(spin.params)): 145 # Relaxation rate. 146 if spin.params[i] == 'Rx': 147 pass 148 149 # Intensity scaling. 150 elif search('^i', spin.params[i]): 151 # Find the position of the first time point. 152 pos = cdp.relax_times.index(min(cdp.relax_times)) 153 154 # Scaling. 155 scaling_matrix[i, i] = 1.0 / average(spin.intensities[pos]) 156 157 # Increment i. 158 i = i + 1 159 160 # Return the scaling matrix. 161 return scaling_matrix

162 163

164 - def _back_calc(self, spin=None, relax_time_id=None):

165 """Back-calculation of peak intensity for the given relaxation time. 166 167 @keyword spin: The spin container. 168 @type spin: SpinContainer instance 169 @keyword relax_time_id: The ID string for the desired relaxation time. 170 @type relax_time_id: str 171 @return: The peak intensity for the desired relaxation time. 172 @rtype: float 173 """ 174 175 # Create the initial parameter vector. 176 param_vector = self._assemble_param_vector(spin=spin) 177 178 # Create a scaling matrix. 179 scaling_matrix = self._assemble_scaling_matrix(spin=spin, scaling=False) 180 181 # The keys. 182 keys = spin.intensities.keys() 183 184 # The peak intensities and times. 185 values = [] 186 errors = [] 187 times = [] 188 for key in keys: 189 values.append(spin.intensities[key]) 190 errors.append(spin.intensity_err[key]) 191 times.append(cdp.relax_times[key]) 192 193 # Initialise the relaxation fit functions. 194 setup(num_params=len(spin.params), num_times=len(cdp.relax_times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_matrix) 195 196 # Make a single function call. This will cause back calculation and the data will be stored in the C module. 197 self._func(param_vector) 198 199 # Get the data back. 200 results = back_calc_I() 201 202 # Return the correct peak height. 203 return results[keys.index(relax_time_id)]

204 205

206 - def _disassemble_param_vector(self, param_vector=None, spin=None, sim_index=None):

207 """Disassemble the parameter vector. 208 209 @keyword param_vector: The parameter vector. 210 @type param_vector: numpy array 211 @keyword spin: The spin data container. 212 @type spin: SpinContainer instance 213 @keyword sim_index: The optional MC simulation index. 214 @type sim_index: int 215 """ 216 217 # Monte Carlo simulations. 218 if sim_index != None: 219 # The relaxation rate. 220 spin.rx_sim[sim_index] = param_vector[0] 221 222 # Initial intensity. 223 spin.i0_sim[sim_index] = param_vector[1] 224 225 # Intensity at infinity. 226 if cdp.curve_type == 'inv': 227 spin.iinf_sim[sim_index] = param_vector[2] 228 229 # Parameter values. 230 else: 231 # The relaxation rate. 232 spin.rx = param_vector[0] 233 234 # Initial intensity. 235 spin.i0 = param_vector[1] 236 237 # Intensity at infinity. 238 if cdp.curve_type == 'inv': 239 spin.iinf = param_vector[2]

240 241

242 - def _func(self, params):

243 """Wrapper function for the C module, for converting numpy arrays. 244 245 @param params: The parameter array from the minimisation code. 246 @type params: numpy array 247 @return: The function value generated by the C module. 248 @rtype: float 249 """ 250 251 # Call the C code. 252 chi2 = func(params.tolist()) 253 254 # Return the chi2 value. 255 return chi2

256 257

258 - def _dfunc(self, params):

259 """Wrapper function for the C module, for converting numpy arrays. 260 261 The currently does nothing. 262 """

263 264

265 - def _d2func(self, params):

266 """Wrapper function for the C module, for converting numpy arrays. 267 268 The currently does nothing. 269 """

270 271

272 - def _grid_search_setup(self, spin=None, param_vector=None, lower=None, upper=None, inc=None, scaling_matrix=None):

273 """The grid search setup function. 274 275 @keyword spin: The spin data container. 276 @type spin: SpinContainer instance 277 @keyword param_vector: The parameter vector. 278 @type param_vector: numpy array 279 @keyword lower: The lower bounds of the grid search which must be equal to the 280 number of parameters in the model. This optional argument is 281 only used when doing a grid search. 282 @type lower: array of numbers 283 @keyword upper: The upper bounds of the grid search which must be equal to the 284 number of parameters in the model. This optional argument is 285 only used when doing a grid search. 286 @type upper: array of numbers 287 @keyword inc: The increments for each dimension of the space for the grid 288 search. The number of elements in the array must equal to the 289 number of parameters in the model. This argument is only used 290 when doing a grid search. 291 @type inc: array of int 292 @keyword scaling_matrix: The scaling matrix. 293 @type scaling_matrix: numpy diagonal matrix 294 @return: A tuple of the grid size and the minimisation options. For the 295 minimisation options, the first dimension corresponds to the 296 model parameter. The second dimension is a list of the number 297 of increments, the lower bound, and upper bound. 298 @rtype: (int, list of lists [int, float, float]) 299 """ 300 301 # The length of the parameter array. 302 n = len(param_vector) 303 304 # Make sure that the length of the parameter array is > 0. 305 if n == 0: 306 raise RelaxError("Cannot run a grid search on a model with zero parameters.") 307 308 # Lower bounds. 309 if lower != None and len(lower) != n: 310 raise RelaxLenError('lower bounds', n) 311 312 # Upper bounds. 313 if upper != None and len(upper) != n: 314 raise RelaxLenError('upper bounds', n) 315 316 # Increments. 317 if isinstance(inc, list) and len(inc) != n: 318 raise RelaxLenError('increment', n) 319 elif isinstance(inc, int): 320 inc = [inc]*n 321 322 # Set up the default bounds. 323 if not lower: 324 # Init. 325 lower = [] 326 upper = [] 327 328 # Loop over the parameters. 329 for i in range(n): 330 # Relaxation rate (from 0 to 20 s^-1). 331 if spin.params[i] == 'Rx': 332 lower.append(0.0) 333 upper.append(20.0) 334 335 # Intensity 336 elif search('^I', spin.params[i]): 337 # Find the ID of the first time point. 338 min_time = min(cdp.relax_times.values()) 339 for key in cdp.relax_times.keys(): 340 if cdp.relax_times[key] == min_time: 341 id = key 342 break 343 344 # Defaults. 345 lower.append(0.0) 346 upper.append(average(spin.intensities[id])) 347 348 # Parameter scaling. 349 for i in range(n): 350 lower[i] = lower[i] / scaling_matrix[i, i] 351 upper[i] = upper[i] / scaling_matrix[i, i] 352 353 return inc, lower, upper

354 355

356 - def _linear_constraints(self, spin=None, scaling_matrix=None):

357 """Set up the relaxation curve fitting linear constraint matrices A and b. 358 359 Standard notation 360 ================= 361 362 The relaxation rate constraints are:: 363 364 Rx >= 0 365 366 The intensity constraints are:: 367 368 I0 >= 0 369 Iinf >= 0 370 371 372 Matrix notation 373 =============== 374 375 In the notation A.x >= b, where A is an matrix of coefficients, x is an array of parameter 376 values, and b is a vector of scalars, these inequality constraints are:: 377 378 | 1 0 0 | | Rx | | 0 | 379 | | | | | | 380 | 1 0 0 | . | I0 | >= | 0 | 381 | | | | | | 382 | 1 0 0 | | Iinf | | 0 | 383 384 385 @keyword spin: The spin data container. 386 @type spin: SpinContainer instance 387 @keyword scaling_matrix: The diagonal, square scaling matrix. 388 @type scaling_matrix: numpy diagonal matrix 389 """ 390 391 # Initialisation (0..j..m). 392 A = [] 393 b = [] 394 n = len(spin.params) 395 zero_array = zeros(n, float64) 396 i = 0 397 j = 0 398 399 # Loop over the parameters. 400 for k in xrange(len(spin.params)): 401 # Relaxation rate. 402 if spin.params[k] == 'Rx': 403 # Rx >= 0. 404 A.append(zero_array * 0.0) 405 A[j][i] = 1.0 406 b.append(0.0) 407 j = j + 1 408 409 # Intensity parameter. 410 elif search('^I', spin.params[k]): 411 # I0, Iinf >= 0. 412 A.append(zero_array * 0.0) 413 A[j][i] = 1.0 414 b.append(0.0) 415 j = j + 1 416 417 # Increment i. 418 i = i + 1 419 420 # Convert to numpy data structures. 421 A = array(A, float64) 422 b = array(b, float64) 423 424 return A, b

425 426

427 - def _model_setup(self, model, params):

428 """Update various model specific data structures. 429 430 @param model: The exponential curve type. 431 @type model: str 432 @param params: A list consisting of the model parameters. 433 @type params: list of str 434 """ 435 436 # Set the model. 437 cdp.curve_type = model 438 439 # Loop over the sequence. 440 for spin in spin_loop(): 441 # Skip deselected spins. 442 if not spin.select: 443 continue 444 445 # Initialise the data structures (if needed). 446 self.data_init(spin) 447 448 # The model and parameter names. 449 spin.model = model 450 spin.params = params

451 452

453 - def _relax_time(self, time=0.0, spectrum_id=None):

454 """Set the relaxation time period associated with a given spectrum. 455 456 @keyword time: The time, in seconds, of the relaxation period. 457 @type time: float 458 @keyword spectrum_id: The spectrum identification string. 459 @type spectrum_id: str 460 """ 461 462 # Test if the spectrum id exists. 463 if spectrum_id not in cdp.spectrum_ids: 464 raise RelaxError("The peak heights corresponding to spectrum id '%s' have not been loaded." % spectrum_id) 465 466 # Initialise the global relaxation time data structure if needed. 467 if not hasattr(cdp, 'relax_times'): 468 cdp.relax_times = {} 469 470 # Add the time at the correct position. 471 cdp.relax_times[spectrum_id] = time

472 473

474 - def _select_model(self, model='exp'):

475 """Function for selecting the model of the exponential curve. 476 477 @keyword model: The exponential curve type. Can be one of 'exp' or 'inv'. 478 @type model: str 479 """ 480 481 # Test if the current pipe exists. 482 pipes.test() 483 484 # Test if the pipe type is set to 'relax_fit'. 485 function_type = cdp.pipe_type 486 if function_type != 'relax_fit': 487 raise RelaxFuncSetupError(specific_setup.get_string(function_type)) 488 489 # Test if sequence data is loaded. 490 if not exists_mol_res_spin_data(): 491 raise RelaxNoSequenceError 492 493 # Two parameter exponential fit. 494 if model == 'exp': 495 print("Two parameter exponential fit.") 496 params = ['Rx', 'I0'] 497 498 # Three parameter inversion recovery fit. 499 elif model == 'inv': 500 print("Three parameter inversion recovery fit.") 501 params = ['Rx', 'I0', 'Iinf'] 502 503 # Invalid model. 504 else: 505 raise RelaxError("The model '" + model + "' is invalid.") 506 507 # Set up the model. 508 self._model_setup(model, params)

509 510

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

512 """Create the Monte Carlo peak intensity data. 513 514 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 515 @type data_id: str 516 @return: The Monte Carlo simulation data. 517 @rtype: list of floats 518 """ 519 520 # Initialise the MC data data structure. 521 mc_data = {} 522 523 # Get the spin container. 524 spin = return_spin(data_id) 525 526 # Skip deselected spins. 527 if not spin.select: 528 return 529 530 # Skip spins which have no data. 531 if not hasattr(spin, 'intensities'): 532 return 533 534 # Test if the model is set. 535 if not hasattr(spin, 'model') or not spin.model: 536 raise RelaxNoModelError 537 538 # Loop over the spectral time points. 539 for id in cdp.relax_times.keys(): 540 # Back calculate the value. 541 value = self._back_calc(spin=spin, relax_time_id=id) 542 543 # Append the value. 544 mc_data[id] = value 545 546 # Return the MC data. 547 return mc_data

548 549

550 - def data_init(self, data_cont, sim=False):

551 """Initialise the spin specific data structures. 552 553 @param data_cont: The spin container. 554 @type data_cont: SpinContainer instance 555 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure. 556 @type sim: bool 557 """ 558 559 # Loop over the data structure names. 560 for name in self.data_names(): 561 # Data structures which are initially empty arrays. 562 list_data = [ 'params' ] 563 if name in list_data: 564 init_data = [] 565 566 # Otherwise initialise the data structure to None. 567 else: 568 init_data = None 569 570 # If the name is not in 'data_cont', add it. 571 if not hasattr(data_cont, name): 572 setattr(data_cont, name, init_data)

573 574

575 - def data_names(self, set='all', error_names=False, sim_names=False):

576 """Return a list of names of data structures. 577 578 Description 579 =========== 580 581 The names are as follows: 582 583 - 'params', an array of the parameter names associated with the model. 584 - 'rx', either the R1 or R2 relaxation rate. 585 - 'i0', the initial intensity. 586 - 'iinf', the intensity at infinity. 587 - 'chi2', chi-squared value. 588 - 'iter', iterations. 589 - 'f_count', function count. 590 - 'g_count', gradient count. 591 - 'h_count', hessian count. 592 - 'warning', minimisation warning. 593 594 595 @keyword set: The set of object names to return. This can be set to 'all' for all names, to 'generic' for generic object names, 'params' for analysis specific parameter names, or to 'min' for minimisation specific object names. 596 @type set: str 597 @keyword error_names: A flag which if True will add the error object names as well. 598 @type error_names: bool 599 @keyword sim_names: A flag which if True will add the Monte Carlo simulation object names as well. 600 @type sim_names: bool 601 @return: The list of object names. 602 @rtype: list of str 603 """ 604 605 # Initialise. 606 names = [] 607 608 # Generic. 609 if set == 'all' or set == 'generic': 610 names.append('params') 611 612 # Parameters. 613 if set == 'all' or set == 'params': 614 names.append('rx') 615 names.append('i0') 616 names.append('iinf') 617 618 # Minimisation statistics. 619 if set == 'all' or set == 'min': 620 names.append('chi2') 621 names.append('iter') 622 names.append('f_count') 623 names.append('g_count') 624 names.append('h_count') 625 names.append('warning') 626 627 # Parameter errors. 628 if error_names and (set == 'all' or set == 'params'): 629 names.append('rx_err') 630 names.append('i0_err') 631 names.append('iinf_err') 632 633 # Parameter simulation values. 634 if sim_names and (set == 'all' or set == 'params'): 635 names.append('rx_sim') 636 names.append('i0_sim') 637 names.append('iinf_sim') 638 639 # Return the names. 640 return names

641 642 643 default_value_doc = ["Relaxation curve fitting default values", """ 644 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. 645 ___________________________________________________________________ 646 | | | | 647 | Data type | Object name | Value | 648 |________________________|_______________|________________________| 649 | | | | 650 | Relaxation rate | 'rx' | 8.0 | 651 | | | | 652 | Initial intensity | 'i0' | 10000.0 | 653 | | | | 654 | Intensity at infinity | 'iinf' | 0.0 | 655 | | | | 656 |________________________|_______________|________________________| 657 658 """] 659 660

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

662 """The exponential curve fitting grid search method. 663 664 @keyword lower: The lower bounds of the grid search which must be equal to the number of parameters in the model. 665 @type lower: array of numbers 666 @keyword upper: The upper bounds of the grid search which must be equal to the number of parameters in the model. 667 @type upper: array of numbers 668 @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. 669 @type inc: array of int 670 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 671 @type constraints: bool 672 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 673 @type verbosity: int 674 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 675 @type sim_index: int 676 """ 677 678 # Minimisation. 679 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, constraints=constraints, verbosity=verbosity, sim_index=sim_index)

680 681

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

683 """Relaxation curve fitting minimisation method. 684 685 @keyword min_algor: The minimisation algorithm to use. 686 @type min_algor: str 687 @keyword min_options: An array of options to be used by the minimisation algorithm. 688 @type min_options: array of str 689 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 690 @type func_tol: None or float 691 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 692 @type grad_tol: None or float 693 @keyword max_iterations: The maximum number of iterations for the algorithm. 694 @type max_iterations: int 695 @keyword constraints: If True, constraints are used during optimisation. 696 @type constraints: bool 697 @keyword scaling: If True, diagonal scaling is enabled during optimisation to allow the problem to be better conditioned. 698 @type scaling: bool 699 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 700 @type verbosity: int 701 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 702 @type sim_index: None or int 703 @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. 704 @type lower: array of numbers 705 @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. 706 @type upper: array of numbers 707 @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. 708 @type inc: array of int 709 """ 710 711 # Test if sequence data is loaded. 712 if not exists_mol_res_spin_data(): 713 raise RelaxNoSequenceError 714 715 # Loop over the sequence. 716 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 717 # Skip deselected spins. 718 if not spin.select: 719 continue 720 721 # Skip spins which have no data. 722 if not hasattr(spin, 'intensities'): 723 continue 724 725 # Create the initial parameter vector. 726 param_vector = self._assemble_param_vector(spin=spin) 727 728 # Diagonal scaling. 729 scaling_matrix = self._assemble_scaling_matrix(spin=spin, scaling=scaling) 730 if len(scaling_matrix): 731 param_vector = dot(inv(scaling_matrix), param_vector) 732 733 # Get the grid search minimisation options. 734 if match('^[Gg]rid', min_algor): 735 inc, lower, upper = self._grid_search_setup(spin=spin, param_vector=param_vector, lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix) 736 737 # Linear constraints. 738 if constraints: 739 A, b = self._linear_constraints(spin=spin, scaling_matrix=scaling_matrix) 740 else: 741 A, b = None, None 742 743 # Print out. 744 if verbosity >= 1: 745 # Get the spin id string. 746 spin_id = generate_spin_id(mol_name, res_num, res_name, spin.num, spin.name) 747 748 # Individual spin print out. 749 if verbosity >= 2: 750 print("\n\n") 751 752 string = "Fitting to spin " + repr(spin_id) 753 print(("\n\n" + string)) 754 print((len(string) * '~')) 755 756 757 # Initialise the function to minimise. 758 ###################################### 759 760 # The keys. 761 keys = spin.intensities.keys() 762 763 # The peak intensities and times. 764 values = [] 765 errors = [] 766 times = [] 767 for key in keys: 768 # The values. 769 if sim_index == None: 770 values.append(spin.intensities[key]) 771 else: 772 values.append(spin.sim_intensities[sim_index][key]) 773 774 # The errors. 775 errors.append(spin.intensity_err[key]) 776 777 # The relaxation times. 778 times.append(cdp.relax_times[key]) 779 780 setup(num_params=len(spin.params), num_times=len(cdp.relax_times), values=values, sd=errors, relax_times=times, scaling_matrix=scaling_matrix.tolist()) 781 782 783 # Setup the minimisation algorithm when constraints are present. 784 ################################################################ 785 786 if constraints and not match('^[Gg]rid', min_algor): 787 algor = min_options[0] 788 else: 789 algor = min_algor 790 791 792 # Levenberg-Marquardt minimisation. 793 ################################### 794 795 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 796 # Reconstruct the error data structure. 797 lm_error = zeros(len(spin.relax_times), float64) 798 index = 0 799 for k in xrange(len(spin.relax_times)): 800 lm_error[index:index+len(relax_error[k])] = relax_error[k] 801 index = index + len(relax_error[k]) 802 803 min_options = min_options + (self.relax_fit.lm_dri, lm_error) 804 805 806 # Minimisation. 807 ############### 808 809 # Grid search. 810 if search('^[Gg]rid', min_algor): 811 results = grid(func=self._func, args=(), num_incs=inc, lower=lower, upper=upper, A=A, b=b, verbosity=verbosity) 812 813 # Unpack the results. 814 param_vector, chi2, iter_count, warning = results 815 f_count = iter_count 816 g_count = 0.0 817 h_count = 0.0 818 819 # Minimisation. 820 else: 821 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) 822 823 # Unpack the results. 824 if results == None: 825 return 826 param_vector, chi2, iter_count, f_count, g_count, h_count, warning = results 827 828 # Scaling. 829 if scaling: 830 param_vector = dot(scaling_matrix, param_vector) 831 832 # Disassemble the parameter vector. 833 self._disassemble_param_vector(param_vector=param_vector, spin=spin, sim_index=sim_index) 834 835 # Monte Carlo minimisation statistics. 836 if sim_index != None: 837 # Chi-squared statistic. 838 spin.chi2_sim[sim_index] = chi2 839 840 # Iterations. 841 spin.iter_sim[sim_index] = iter_count 842 843 # Function evaluations. 844 spin.f_count_sim[sim_index] = f_count 845 846 # Gradient evaluations. 847 spin.g_count_sim[sim_index] = g_count 848 849 # Hessian evaluations. 850 spin.h_count_sim[sim_index] = h_count 851 852 # Warning. 853 spin.warning_sim[sim_index] = warning 854 855 856 # Normal statistics. 857 else: 858 # Chi-squared statistic. 859 spin.chi2 = chi2 860 861 # Iterations. 862 spin.iter = iter_count 863 864 # Function evaluations. 865 spin.f_count = f_count 866 867 # Gradient evaluations. 868 spin.g_count = g_count 869 870 # Hessian evaluations. 871 spin.h_count = h_count 872 873 # Warning. 874 spin.warning = warning

875 876

877 - def overfit_deselect(self):

878 """Deselect spins which have insufficient data to support minimisation.""" 879 880 # Print out. 881 print("\n\nOver-fit spin deselection.\n") 882 883 # Test the sequence data exists. 884 if not exists_mol_res_spin_data(): 885 raise RelaxNoSequenceError 886 887 # Loop over spin data. 888 for spin, spin_id in spin_loop(return_id=True): 889 # Skip deselected spins. 890 if not spin.select: 891 continue 892 893 # Check if data exists. 894 if not hasattr(spin, 'intensities'): 895 warn(RelaxDeselectWarning(spin_id, 'missing intensity data')) 896 spin.select = False 897 898 # Require 3 or more data points. 899 elif len(spin.intensities) < 3: 900 warn(RelaxDeselectWarning(spin_id, 'insufficient data, 3 or more data points are required')) 901 spin.select = False

902 903

904 - def return_data(self, spin):

905 """Function for returning the peak intensity data structure. 906 907 @param spin: The spin container. 908 @type spin: SpinContainer instance 909 @return: The peak intensity data structure. 910 @rtype: list of float 911 """ 912 913 return spin.intensities

914 915 916 return_data_name_doc = ["Relaxation curve fitting data type string matching patterns", """ 917 ____________________________________________________________ 918 | | | 919 | Data type | Object name | 920 |___________________________________|______________________| 921 | | | 922 | Relaxation rate | 'rx' | 923 | | | 924 | Peak intensities (series) | 'intensities' | 925 | | | 926 | Initial intensity | 'i0' | 927 | | | 928 | Intensity at infinity | 'iinf' | 929 | | | 930 | Relaxation period times (series) | 'relax_times' | 931 |___________________________________|______________________| 932 933 """] 934 935

936 - def return_error(self, data_id):

937 """Return the standard deviation data structure. 938 939 @param data_id: The spin identification string, as yielded by the base_data_loop() generator 940 method. 941 @type data_id: str 942 @return: The standard deviation data structure. 943 @rtype: list of float 944 """ 945 946 # Get the spin container. 947 spin = return_spin(data_id) 948 949 # Return the error list. 950 return spin.intensity_err

951 952

953 - def return_units(self, param):

954 """Dummy method which returns None as the stats have no units. 955 956 @param param: The name of the parameter to return the units string for. 957 @type param: str 958 @return: Nothing. 959 @rtype: None 960 """ 961 962 # Unitless. 963 return None

964 965 966 set_doc = ["Relaxation curve fitting set details", """ 967 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. 968 """] 969 970

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

972 """Pack the Monte Carlo simulation data. 973 974 @param data_id: The spin identification string, as yielded by the base_data_loop() generator method. 975 @type data_id: str 976 @param sim_data: The Monte Carlo simulation data. 977 @type sim_data: list of float 978 """ 979 980 # Get the spin container. 981 spin = return_spin(data_id) 982 983 # Create the data structure. 984 spin.sim_intensities = sim_data

Source Code for Module specific_fns.relax_fit