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

Source Code for Module specific_analyses.relax_disp.api

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2004-2016 Edward d'Auvergne                                   # 
   4  # Copyright (C) 2009 Sebastien Morin                                          # 
   5  # Copyright (C) 2013-2014 Troels E. Linnet                                    # 
   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 relaxation dispersion API object.""" 
  26   
  27  # Python module imports. 
  28  import bmrblib 
  29  from copy import deepcopy 
  30  from numpy import int32, sqrt, zeros 
  31  from re import match, search 
  32  import string 
  33  import sys 
  34  from types import MethodType 
  35   
  36  # relax module imports. 
  37  from lib.arg_check import is_list, is_str_list 
  38  from lib.dispersion.variables import EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ, MODEL_LIST_MMQ, MODEL_R2EFF, PARAMS_R20 
  39  from lib.errors import RelaxError, RelaxImplementError 
  40  from lib.text.sectioning import subsection 
  41  from multi import Processor_box 
  42  from pipe_control import pipes, sequence 
  43  from pipe_control.exp_info import bmrb_write_citations, bmrb_write_methods, bmrb_write_software 
  44  from pipe_control.mol_res_spin import bmrb_write_entity, check_mol_res_spin_data, get_molecule_names, return_spin, spin_loop 
  45  from pipe_control.pipes import check_pipe 
  46  from pipe_control.spectrometer import check_spectrometer_setup 
  47  from pipe_control.sequence import return_attached_protons 
  48  from specific_analyses.api_base import API_base 
  49  from specific_analyses.api_common import API_common 
  50  from specific_analyses.relax_disp.checks import check_model_type 
  51  from specific_analyses.relax_disp.data import average_intensity, calc_rotating_frame_params, find_intensity_keys, generate_r20_key, has_exponential_exp_type, has_proton_mmq_cpmg, loop_cluster, loop_exp_frq, loop_exp_frq_offset_point, loop_time, pack_back_calc_r2eff, return_param_key_from_data, spin_ids_to_containers 
  52  from specific_analyses.relax_disp.optimisation import Disp_memo, Disp_minimise_command, back_calc_peak_intensities, back_calc_r2eff, calculate_r2eff, minimise_r2eff 
  53  from specific_analyses.relax_disp.parameter_object import Relax_disp_params 
  54  from specific_analyses.relax_disp.parameters import get_param_names, get_value, loop_parameters, param_index_to_param_info, param_num, r1_setup 
  55   
  56   
57 -class Relax_disp(API_base, API_common):
58 """Class containing functions for relaxation dispersion curve fitting.""" 59 60 # Class variable for storing the class instance (for the singleton design pattern). 61 instance = None 62
63 - def __init__(self):
64 """Initialise the class by placing API_common methods into the API.""" 65 66 # Place methods into the API. 67 self.data_init = self._data_init_spin 68 self.model_type = self._model_type_local 69 self.return_conversion_factor = self._return_no_conversion_factor 70 self.return_value = self.return_value 71 72 # Place a copy of the parameter list object in the instance namespace. 73 self._PARAMS = Relax_disp_params()
74 75
76 - def base_data_loop(self):
77 """Custom generator method for looping over the base data. 78 79 For the R2eff model, the base data is the peak intensity data defining a single exponential curve. For all other models, the base data is the R2eff/R1rho values for individual spins. 80 81 82 @return: For the R2eff model, a tuple of the spin container and the exponential curve identifying key (the CPMG frequency or R1rho spin-lock field strength). For all other models, the spin container and ID from the spin loop. 83 @rtype: (tuple of SpinContainer instance and float) or (SpinContainer instance and str) 84 """ 85 86 # The R2eff model data (the base data is peak intensities). 87 if cdp.model_type == MODEL_R2EFF: 88 # Loop over the sequence. 89 for spin, spin_id in spin_loop(return_id=True): 90 # Skip deselected spins. 91 if not spin.select: 92 continue 93 94 # Skip spins with no peak intensity data. 95 if not hasattr(spin, 'peak_intensity'): 96 continue 97 98 # Loop over each spectrometer frequency and dispersion point. 99 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 100 yield spin, spin_id, exp_type, frq, offset, point 101 102 # All other models (the base data is the R2eff/R1rho values). 103 else: 104 # 1H MMQ flag. 105 proton_mmq_flag = has_proton_mmq_cpmg() 106 107 # Loop over the sequence. 108 for spin, spin_id in spin_loop(return_id=True): 109 # Skip deselected spins. 110 if not spin.select: 111 continue 112 113 # Skip protons for MMQ data. 114 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 115 continue 116 117 # Get the attached proton. 118 proton = None 119 if proton_mmq_flag: 120 proton = return_attached_protons(spin_id)[0] 121 122 # Skip spins with no R2eff/R1rho values. 123 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'): 124 continue 125 126 # Yield the spin container and ID. 127 yield spin, spin_id
128 129
130 - def bmrb_write(self, file_path, version=None):
131 """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file. 132 133 @param file_path: The full file path. 134 @type file_path: str 135 @keyword version: The BMRB NMR-STAR dictionary format to output to. 136 @type version: str 137 """ 138 139 # This function is not yet implemented, as it would require a re-write of the relax_data bmrb_write(star) function, and proper handling of cdp.ri_ids. 140 # It was also not readily possible to find examples of submitted CPMG data in the BMRB database. 141 142 # Not implemented. 143 raise RelaxImplementError('bmrb_write') 144 145 # Checks. 146 check_spectrometer_setup(escalate=2) 147 148 # Alias the current data pipe. 149 cdp = pipes.get_pipe() 150 151 # Initialise the NMR-STAR data object. 152 star = bmrblib.create_nmr_star('relax_relaxation_dispersion_results', file_path, version) 153 154 # Initialise the spin specific data lists. 155 mol_name_list = [] 156 res_num_list = [] 157 res_name_list = [] 158 atom_name_list = [] 159 160 isotope_list = [] 161 element_list = [] 162 163 chi2_list = [] 164 model_list = [] 165 166 # Store the spin specific data in lists for later use. 167 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True): 168 # Skip the protons. 169 if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'): 170 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id)) 171 continue 172 173 # Check the data for None (not allowed in BMRB!). 174 if res_num == None: 175 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id) 176 if res_name == None: 177 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id) 178 if spin.name == None: 179 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id) 180 if not hasattr(spin, 'isotope') or spin.isotope == None: 181 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id) 182 if not hasattr(spin, 'element') or spin.element == None: 183 raise RelaxError("For the BMRB, the spin element type of '%s' must be specified. Please use the spin user function for setting the element type." % spin_id) 184 185 # The molecule/residue/spin info. 186 mol_name_list.append(mol_name) 187 res_num_list.append(res_num) 188 res_name_list.append(res_name) 189 atom_name_list.append(spin.name) 190 191 # The nuclear isotope. 192 if hasattr(spin, 'isotope'): 193 isotope_list.append(int(spin.isotope.strip(string.ascii_letters))) 194 else: 195 isotope_list.append(None) 196 197 # The element. 198 if hasattr(spin, 'element'): 199 element_list.append(spin.element) 200 else: 201 element_list.append(None) 202 203 # Opt stats. 204 if hasattr(spin, 'chi2'): 205 chi2_list.append(spin.chi2) 206 else: 207 chi2_list.append(None) 208 209 # Model-free model. 210 model_list.append(spin.model) 211 212 # Convert the molecule names into the entity IDs. 213 entity_ids = zeros(len(mol_name_list), int32) 214 mol_names = get_molecule_names() 215 for i in range(len(mol_name_list)): 216 for j in range(len(mol_names)): 217 if mol_name_list[i] == mol_names[j]: 218 entity_ids[i] = j+1 219 220 221 # Create Supergroup 2 : The citations. 222 ###################################### 223 224 # Generate the citations saveframe. 225 bmrb_write_citations(star) 226 227 228 # Create Supergroup 3 : The molecular assembly saveframes. 229 ########################################################## 230 231 # Generate the entity saveframe. 232 bmrb_write_entity(star) 233 234 235 # Create Supergroup 4: The experimental descriptions saveframes. 236 ################################################################# 237 238 # Generate the method saveframes. 239 bmrb_write_methods(star) 240 241 # Generate the software saveframe. 242 software_ids, software_labels = bmrb_write_software(star)
243 244
245 - def calculate(self, spin_id=None, scaling_matrix=None, verbosity=1, sim_index=None):
246 """Calculate the model chi-squared value or the R2eff values for fixed time period data. 247 248 @keyword spin_id: The spin identification string. 249 @type spin_id: None or str 250 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 251 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 252 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 253 @type verbosity: int 254 @keyword sim_index: The MC simulation index (unused). 255 @type sim_index: None 256 """ 257 258 # Data checks. 259 check_pipe() 260 check_mol_res_spin_data() 261 check_model_type() 262 263 # Special exponential curve-fitting for the R2eff model. 264 if cdp.model_type == MODEL_R2EFF: 265 calculate_r2eff() 266 267 # Calculate the chi-squared value. 268 else: 269 # 1H MMQ flag. 270 proton_mmq_flag = has_proton_mmq_cpmg() 271 272 # Loop over the spin blocks. 273 model_index = -1 274 for spin_ids in self.model_loop(): 275 # Increment the model index. 276 model_index += 1 277 278 # The spin containers. 279 spins = spin_ids_to_containers(spin_ids) 280 281 # Skip deselected clusters. 282 skip = True 283 for spin in spins: 284 if spin.select: 285 skip = False 286 if skip: 287 continue 288 289 # The back calculated values. 290 back_calc = back_calc_r2eff(spins=spins, spin_ids=spin_ids, store_chi2=True) 291 292 # Pack the data. 293 for i, spin_id in enumerate(spin_ids): 294 spin = spins[i] 295 pack_back_calc_r2eff(spin=spin, spin_id=spin_id, si=i, back_calc=back_calc, proton_mmq_flag=proton_mmq_flag)
296 297
298 - def constraint_algorithm(self):
299 """Return the 'Log barrier' optimisation constraint algorithm. 300 301 @return: The 'Log barrier' constraint algorithm. 302 @rtype: str 303 """ 304 305 # The log barrier algorithm, as required by minfx. 306 return 'Log barrier'
307 308
309 - def create_mc_data(self, data_id):
310 """Create the Monte Carlo peak intensity data. 311 312 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 313 @type data_id: SpinContainer instance and float 314 @return: The Monte Carlo simulation data. 315 @rtype: list of floats 316 """ 317 318 # The R2eff model (with peak intensity base data). 319 if cdp.model_type == MODEL_R2EFF: 320 # Unpack the data. 321 spin, spin_id, exp_type, frq, offset, point = data_id 322 323 # Back calculate the peak intensities. 324 values = back_calc_peak_intensities(spin=spin, spin_id=spin_id, exp_type=exp_type, frq=frq, offset=offset, point=point) 325 326 # All other models (with R2eff/R1rho base data). 327 else: 328 # 1H MMQ flag. 329 proton_mmq_flag = has_proton_mmq_cpmg() 330 331 # Unpack the data. 332 spin, spin_id = data_id 333 334 # Back calculate the R2eff/R1rho data. 335 back_calc = back_calc_r2eff(spins=[spin], spin_ids=[spin_id]) 336 337 # Get the attached proton data. 338 if proton_mmq_flag: 339 proton = return_attached_protons(spin_id)[0] 340 proton_back_calc = back_calc_r2eff(spins=[proton], spin_ids=[spin_id]) 341 342 # Convert to a dictionary matching the R2eff data structure. 343 values = {} 344 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 345 # Alias the correct data. 346 current_bc = back_calc 347 current_spin = spin 348 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 349 current_spin = proton 350 current_bc = proton_back_calc 351 352 # The parameter key. 353 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 354 355 # Skip missing data. 356 if not hasattr(current_spin, 'r2eff') or param_key not in current_spin.r2eff: 357 continue 358 359 # Store the result. 360 values[param_key] = back_calc[ei][0][mi][oi][di] 361 362 # Return the MC data. 363 return values
364 365
366 - def deselect(self, sim_index=None, model_info=None):
367 """Deselect models or simulations. 368 369 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will. 370 @type sim_index: None or int 371 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 372 @type model_info: list of SpinContainer instances, list of str 373 """ 374 375 # Loop over all the spins and deselect them. 376 for spin_id in model_info: 377 # Get the spin. 378 spin = return_spin(spin_id) 379 380 # Spin deselection. 381 if sim_index == None: 382 spin.select = False 383 384 # Simulation deselection. 385 else: 386 spin.select_sim[sim_index] = False
387 388
389 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):
390 """Duplicate the data specific to a single model. 391 392 @keyword pipe_from: The data pipe to copy the data from. 393 @type pipe_from: str 394 @keyword pipe_to: The data pipe to copy the data to. 395 @type pipe_to: str 396 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 397 @type model_info: list of SpinContainer instances, list of str 398 @keyword global_stats: The global statistics flag. 399 @type global_stats: bool 400 @keyword verbose: A flag which if True will cause info to be printed out. 401 @type verbose: bool 402 """ 403 404 # First create the pipe_to data pipe, if it doesn't exist, but don't switch to it. 405 if not pipes.has_pipe(pipe_to): 406 pipes.create(pipe_to, pipe_type='relax_disp', switch=False) 407 408 # Get the data pipes. 409 dp_from = pipes.get_pipe(pipe_from) 410 dp_to = pipes.get_pipe(pipe_to) 411 412 # Duplicate all non-sequence specific data. 413 for data_name in dir(dp_from): 414 # Skip the container objects. 415 if data_name in ['mol', 'interatomic', 'structure', 'exp_info', 'result_files']: 416 continue 417 418 # Skip dispersion specific parameters. 419 if data_name in ['model']: 420 continue 421 422 # Skip special objects. 423 if search('^_', data_name) or data_name in dp_from.__class__.__dict__: 424 continue 425 426 # Get the original object. 427 data_from = getattr(dp_from, data_name) 428 429 # The data already exists. 430 if hasattr(dp_to, data_name): 431 # Get the object in the target pipe. 432 data_to = getattr(dp_to, data_name) 433 434 # The data must match! 435 if data_from != data_to: 436 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 437 438 # Skip the data. 439 continue 440 441 # Duplicate the data. 442 setattr(dp_to, data_name, deepcopy(data_from)) 443 444 # No sequence data, so skip the rest. 445 if dp_from.mol.is_empty(): 446 return 447 448 # Duplicate the sequence data if it doesn't exist. 449 if dp_to.mol.is_empty(): 450 sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose) 451 452 # Loop over the cluster. 453 for id in model_info: 454 # The original spin container. 455 spin = return_spin(id, pipe=pipe_from) 456 457 # Duplicate the spin specific data. 458 for name in dir(spin): 459 # Skip special objects. 460 if search('^__', name): 461 continue 462 463 # Get the object. 464 obj = getattr(spin, name) 465 466 # Skip methods. 467 if isinstance(obj, MethodType): 468 continue 469 470 # Duplicate the object. 471 new_obj = deepcopy(getattr(spin, name)) 472 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj)
473 474
475 - def eliminate(self, name, value, args, sim=None, model_info=None):
476 """Relaxation dispersion model elimination, parameter by parameter. 477 478 @param name: The parameter name. 479 @type name: str 480 @param value: The parameter value. 481 @type value: float 482 @param args: The c1 and c2 elimination constant overrides. 483 @type args: None or tuple of float 484 @keyword sim: The Monte Carlo simulation index. 485 @type sim: int 486 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 487 @type model_info: list of SpinContainer instances, list of str 488 @return: True if the model is to be eliminated, False otherwise. 489 @rtype: bool 490 """ 491 492 # Skip the R2eff model parameters. 493 if name in ['r2eff', 'i0']: 494 return False 495 496 # Default limits. 497 c1 = 0.501 498 c2 = 0.999 499 c3 = 1.0 500 501 # Depack the arguments. 502 if args != None: 503 c1, c2, c3 = args 504 505 # Elimination text. 506 elim_text = "Data pipe '%s': The %s parameter of %.5f is %s, eliminating " % (pipes.cdp_name(), name, value, "%s") 507 if sim == None: 508 elim_text += "the spin cluster %s." % model_info 509 else: 510 elim_text += "simulation %i of the spin cluster %s." % (sim, model_info) 511 512 # The pA parameter. 513 if name == 'pA': 514 if value < c1: 515 print(elim_text % ("less than %.5f" % c1)) 516 return True 517 if value > c2: 518 print(elim_text % ("greater than %.5f" % c2)) 519 return True 520 521 # The tex parameter. 522 if name == 'tex': 523 if value > c3: 524 print(elim_text % ("greater than %.5f" % c3)) 525 return True 526 527 # Accept model. 528 return False
529 530
531 - def get_param_names(self, model_info=None):
532 """Return a vector of parameter names. 533 534 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 535 @type model_info: list of SpinContainer instances, list of str 536 @return: The vector of parameter names. 537 @rtype: list of str 538 """ 539 540 # Get the spin containers. 541 spins = [] 542 for spin_id in model_info: 543 # Get the spin. 544 spin = return_spin(spin_id) 545 546 # Skip deselected spins. 547 if not spin.select: 548 continue 549 550 # Add the spin. 551 spins.append(spin) 552 553 # No spins. 554 if not len(spins): 555 return None 556 557 # Call the get_param_names() function. 558 return get_param_names(spins)
559 560
561 - def get_param_values(self, model_info=None, sim_index=None):
562 """Return a vector of parameter values. 563 564 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 565 @type model_info: list of SpinContainer instances, list of str 566 @keyword sim_index: The Monte Carlo simulation index. 567 @type sim_index: int 568 @return: The vector of parameter values. 569 @rtype: list of str 570 """ 571 572 # Get the spin containers. 573 spins = [] 574 for spin_id in model_info: 575 # Get the spin. 576 spin = return_spin(spin_id) 577 578 # Skip deselected spins. 579 if not spin.select: 580 continue 581 582 # Add the spin. 583 spins.append(spin) 584 585 # No spins. 586 if not len(spins): 587 return None 588 589 # Set up the R1 parameter, if needed. 590 r1_setup() 591 592 # Loop over the parameters of the cluster, fetching their values. 593 values = [] 594 for param_name, param_index, si, r20_key in loop_parameters(spins=spins): 595 values.append(get_value(spins=spins, sim_index=sim_index, param_name=param_name, spin_index=si, r20_key=r20_key)) 596 597 # Return all values. 598 return values
599 600
601 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
602 """The relaxation dispersion curve fitting grid search function. 603 604 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model. 605 @type lower: list of list of numbers 606 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model. 607 @type upper: list of list of numbers 608 @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. 609 @type inc: list of list of int 610 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 611 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 612 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 613 @type constraints: bool 614 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 615 @type verbosity: int 616 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 617 @type sim_index: int 618 """ 619 620 # Minimisation. 621 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
622 623
624 - def map_bounds(self, param, spin_id=None):
625 """Create bounds for the OpenDX mapping function. 626 627 @param param: The name of the parameter to return the lower and upper bounds of. 628 @type param: str 629 @param spin_id: The spin identification string (unused). 630 @type spin_id: None 631 @return: The upper and lower bounds of the parameter. 632 @rtype: list of float 633 """ 634 635 # Is the parameter is valid? 636 if not self._PARAMS.contains(param): 637 raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param) 638 639 # Return the spin. 640 spin = return_spin(spin_id) 641 642 # Loop over each spectrometer frequency and dispersion point to collect param_keys. 643 param_keys = [] 644 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 645 # The parameter key. 646 param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 647 648 # Collect the key. 649 param_keys.append(param_key) 650 651 # The initial parameter vector. 652 param_vector = [] 653 654 # Collect param_names. 655 param_names = [] 656 for param_name, param_index, si, r20_key in loop_parameters(spins=[spin]): 657 # Add to the param vector. 658 param_vector.append([0.0]) 659 660 # Collect parameter names. 661 param_names.append(param_name) 662 663 # Loop over the parameter names. 664 for i in range(len(param_names)): 665 # Test if the parameter is in the list: 666 667 if param_names[i] == param: 668 return [self._PARAMS.grid_lower(param, incs=0, model_info=[spin_id]), self._PARAMS.grid_upper(param, incs=0, model_info=[spin_id])]
669 670
671 - 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):
672 """Relaxation dispersion curve fitting function. 673 674 @keyword min_algor: The minimisation algorithm to use. 675 @type min_algor: str 676 @keyword min_options: An array of options to be used by the minimisation algorithm. 677 @type min_options: array of str 678 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 679 @type func_tol: None or float 680 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 681 @type grad_tol: None or float 682 @keyword max_iterations: The maximum number of iterations for the algorithm. 683 @type max_iterations: int 684 @keyword constraints: If True, constraints are used during optimisation. 685 @type constraints: bool 686 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 687 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 688 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 689 @type verbosity: int 690 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 691 @type sim_index: None or int 692 @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. 693 @type lower: list of list of numbers 694 @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. 695 @type upper: list of list of numbers 696 @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. 697 @type inc: list of list of int 698 """ 699 700 # Data checks. 701 check_mol_res_spin_data() 702 check_model_type() 703 704 # Check the optimisation algorithm. 705 algor = min_algor 706 if min_algor == 'Log barrier': 707 algor = min_options[0] 708 709 allow = False 710 # Check the model type. 711 if hasattr(cdp, 'model_type'): 712 # Set the model type: 713 model_type = cdp.model_type 714 715 if model_type == MODEL_R2EFF: 716 if match('^[Gg]rid$', algor): 717 allow = True 718 719 elif match('^[Ss]implex$', algor): 720 allow = True 721 722 # Quasi-Newton BFGS minimisation. 723 elif match('^[Bb][Ff][Gg][Ss]$', algor): 724 allow = True 725 726 # Newton minimisation. 727 elif match('^[Nn]ewton$', algor): 728 allow = True 729 730 # Newton minimisation. 731 elif match('^[Nn]ewton$', algor): 732 allow = True 733 734 # Constrained method, Method of Multipliers. 735 elif match('^[Mm][Oo][Mm]$', algor) or match('[Mm]ethod of [Mm]ultipliers$', algor): 736 allow = True 737 738 # Constrained method, Logarithmic barrier function. 739 elif match('^[Ll]og [Bb]arrier$', algor): 740 allow = True 741 742 # If the Jacobian and Hessian matrix have not been specified for fitting, 'simplex' should be used. 743 else: 744 if match('^[Gg]rid$', algor): 745 allow = True 746 747 elif match('^[Ss]implex$', algor): 748 allow = True 749 750 # Do not allow, if no model has been specified. 751 else: 752 model_type = 'None' 753 # Do not allow. 754 allow = False 755 756 if not allow: 757 raise RelaxError("Minimisation algorithm '%s' is not allowed, since function gradients for model '%s' is not implemented. Only the 'simplex' minimisation algorithm is supported for the relaxation dispersion analysis of this model."%(algor, model_type)) 758 759 # Initialise some empty data pipe structures so that the target function set up does not fail. 760 if not hasattr(cdp, 'cpmg_frqs_list'): 761 cdp.cpmg_frqs_list = [] 762 if not hasattr(cdp, 'spin_lock_nu1_list'): 763 cdp.spin_lock_nu1_list = [] 764 765 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 766 processor_box = Processor_box() 767 processor = processor_box.processor 768 769 # The number of time points for the exponential curves (if present). 770 num_time_pts = 1 771 if hasattr(cdp, 'num_time_pts'): 772 num_time_pts = cdp.num_time_pts 773 774 # Number of spectrometer fields. 775 fields = [None] 776 field_count = 1 777 if hasattr(cdp, 'spectrometer_frq'): 778 fields = cdp.spectrometer_frq_list 779 field_count = cdp.spectrometer_frq_count 780 781 # Loop over the spin blocks. 782 model_index = -1 783 for spin_ids in self.model_loop(): 784 # Increment the model index. 785 model_index += 1 786 787 # The spin containers. 788 spins = spin_ids_to_containers(spin_ids) 789 790 # Skip deselected clusters. 791 skip = True 792 for spin in spins: 793 if spin.select: 794 skip = False 795 if skip: 796 continue 797 798 # Alias the grid options. 799 lower_i, upper_i, inc_i = None, None, None 800 if min_algor == 'grid': 801 lower_i = lower[model_index] 802 upper_i = upper[model_index] 803 inc_i = inc[model_index] 804 805 # Special exponential curve-fitting for the R2eff model. 806 if cdp.model_type == MODEL_R2EFF: 807 # Sanity checks. 808 if not has_exponential_exp_type(): 809 raise RelaxError("The R2eff model with the fixed time period dispersion experiments cannot be optimised.") 810 811 # Optimisation. 812 minimise_r2eff(spins=spins, spin_ids=spin_ids, min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, scaling_matrix=scaling_matrix[model_index], verbosity=verbosity, sim_index=sim_index, lower=lower_i, upper=upper_i, inc=inc_i) 813 814 # Skip the rest. 815 continue 816 817 # Set up the slave command object. 818 command = Disp_minimise_command(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix[model_index], min_algor=min_algor, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, max_iterations=max_iterations, constraints=constraints, verbosity=verbosity, lower=lower_i, upper=upper_i, inc=inc_i, fields=fields, param_names=get_param_names(spins=spins, full=True)) 819 820 # Set up the memo. 821 memo = Disp_memo(spins=spins, spin_ids=spin_ids, sim_index=sim_index, scaling_matrix=scaling_matrix[model_index], verbosity=verbosity) 822 823 # Add the slave command and memo to the processor queue. 824 processor.add_to_queue(command, memo)
825 826
827 - def model_desc(self, model_info=None):
828 """Return a description of the model. 829 830 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 831 @type model_info: list of SpinContainer instances, list of str 832 @return: The model description. 833 @rtype: str 834 """ 835 836 # The model loop is over the spin clusters, so return a description of the cluster. 837 return "The spin cluster %s." % model_info
838 839
840 - def model_loop(self):
841 """Loop over the spin groupings for one model applied to multiple spins. 842 843 @return: The list of spins per block will be yielded, as well as the list of spin IDs. 844 @rtype: tuple of list of SpinContainer instances and list of str 845 """ 846 847 # Loop over individual spins for the R2eff model. 848 if not hasattr(cdp, 'model_type') or cdp.model_type == MODEL_R2EFF: 849 # The spin loop. 850 for spin, spin_id in spin_loop(return_id=True): 851 # Skip deselected spins 852 if not spin.select: 853 continue 854 855 # Yield the spin ID as a list. 856 yield [spin_id] 857 858 # The cluster loop. 859 else: 860 for spin_ids in loop_cluster(skip_desel=False): 861 yield spin_ids
862 863
864 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):
865 """Return the k, n, and chi2 model statistics. 866 867 k - number of parameters. 868 n - number of data points. 869 chi2 - the chi-squared value. 870 871 872 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 873 @type model_info: list of SpinContainer instances, list of str 874 @keyword spin_id: The spin ID string to override the model_info argument. This is ignored in the N-state model. 875 @type spin_id: None or str 876 @keyword global_stats: A parameter which determines if global or local statistics are returned. For the N-state model, this argument is ignored. 877 @type global_stats: None or bool 878 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2). 879 @rtype: tuple of (int, int, float) 880 """ 881 882 # Get the spin containers (the spin ID overrides the model info). 883 if spin_id != None: 884 spins = [return_spin(spin_id)] 885 else: 886 spins = spin_ids_to_containers(model_info) 887 888 # The number of parameters for the cluster. 889 k = param_num(spins=spins) 890 891 # The number of points from all spins. 892 n = 0 893 for spin in spins: 894 # Skip deselected spins. 895 if not spin.select: 896 continue 897 898 n += len(spin.r2eff) 899 900 # Take the chi-squared from the first spin of the cluster (which has a value). 901 chi2 = None 902 for spin in spins: 903 # Skip deselected spins. 904 if not spin.select: 905 continue 906 907 if hasattr(spin, 'chi2'): 908 chi2 = spin.chi2 909 break 910 911 # Return the values. 912 return k, n, chi2
913 914
915 - def overfit_deselect(self, data_check=True, verbose=True):
916 """Deselect spins which have insufficient data to support minimisation. 917 918 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 919 @type data_check: bool 920 @keyword verbose: A flag which if True will allow printouts. 921 @type verbose: bool 922 """ 923 924 # Test the sequence data exists and the model is setup. 925 check_mol_res_spin_data() 926 check_model_type() 927 928 # 1H MMQ flag. 929 proton_mmq_flag = has_proton_mmq_cpmg() 930 931 # Loop over spin data. 932 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 933 # Skip protons for MMQ data. 934 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 935 continue 936 937 # Get the attached proton. 938 proton = None 939 if proton_mmq_flag: 940 # Get all protons. 941 proton_spins = return_attached_protons(spin_id) 942 943 # Only one allowed. 944 if len(proton_spins) > 1: 945 print("Multiple protons attached to the spin '%s', but one one attached proton is supported for the MMQ-type models." % spin_id) 946 spin.select = False 947 continue 948 949 # Alias the single proton. 950 if len(proton_spins): 951 proton = proton_spins[0] 952 953 # Check if data exists. 954 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'): 955 print("No R2eff data could be found, deselecting the '%s' spin." % spin_id) 956 spin.select = False 957 continue
958 959
960 - def print_model_title(self, prefix=None, model_info=None):
961 """Print out the model title. 962 963 @keyword prefix: The starting text of the title. This should be printed out first, followed by the model information text. 964 @type prefix: str 965 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 966 @type model_info: list of SpinContainer instances, list of str 967 """ 968 969 # The printout. 970 subsection(file=sys.stdout, text=prefix+"the spin block %s"%model_info, prespace=2)
971 972
973 - def return_data(self, data_id=None):
974 """Return the peak intensity data structure. 975 976 @param data_id: The spin ID string, as yielded by the base_data_loop() generator method. 977 @type data_id: str 978 @return: The peak intensity data structure. 979 @rtype: list of float 980 """ 981 982 # The R2eff model. 983 if cdp.model_type == MODEL_R2EFF: 984 # Unpack the data. 985 spin, spin_id, exp_type, frq, offset, point = data_id 986 987 # Return the data. 988 return spin.peak_intensity 989 990 # All other models. 991 else: 992 raise RelaxImplementError
993 994
995 - def return_error(self, data_id=None):
996 """Return the standard deviation data structure. 997 998 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 999 @type data_id: SpinContainer instance and float 1000 @return: The standard deviation data structure. 1001 @rtype: list of float 1002 """ 1003 1004 # The R2eff model. 1005 if cdp.model_type == MODEL_R2EFF: 1006 # Unpack the data. 1007 spin, spin_id, exp_type, frq, offset, point = data_id 1008 1009 # Generate the data structure to return. 1010 errors = [] 1011 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 1012 errors.append(average_intensity(spin=spin, exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, error=True)) 1013 1014 # All other models. 1015 else: 1016 # Unpack the data. 1017 spin, spin_id = data_id 1018 1019 # 1H MMQ flag. 1020 proton_mmq_flag = has_proton_mmq_cpmg() 1021 1022 # Get the attached proton. 1023 proton = None 1024 if proton_mmq_flag: 1025 proton = return_attached_protons(spin_id)[0] 1026 1027 # The errors. 1028 r2eff_err = {} 1029 if hasattr(spin, 'r2eff_err'): 1030 r2eff_err.update(spin.r2eff_err) 1031 if hasattr(proton, 'r2eff_err'): 1032 r2eff_err.update(proton.r2eff_err) 1033 return r2eff_err 1034 1035 # Return the error list. 1036 return errors
1037 1038
1039 - def return_error_red_chi2(self, data_id=None):
1040 """Return the standard deviation data structure, where standard deviation is from the overall gauss distribution described by the STD_fit of the goodness of fit, where STD_fit = sqrt(chi2/(N-p)) 1041 1042 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 1043 @type data_id: SpinContainer instance and float 1044 @return: The standard deviation data structure. 1045 @rtype: list of float 1046 """ 1047 1048 # The R2eff model. 1049 if cdp.model_type == MODEL_R2EFF: 1050 raise RelaxError("Drawing errors from the gauss distribution described by the STD_fit of the goodness of fit, is not possible for the '%s' model."%MODEL_R2EFF) 1051 1052 # Get the errors structure as above. 1053 errors = self.return_error(data_id=data_id) 1054 1055 # Unpack the data. 1056 spin, spin_id = data_id 1057 1058 # Loop over the spin groupings for the model. 1059 for spin_ids in self.model_loop(): 1060 # If the spin of interest is in the returned spin cluster. 1061 if spin_id in spin_ids: 1062 # Get the statistics 1063 k, n, chi2 = self.model_statistics(model_info=spin_ids) 1064 1065 # Calculate degrees of freedom. 1066 dof = n - k 1067 1068 # Calculate reduced chi2, or named as the variance of the squared residuals. 1069 red_chi2 = chi2 / float(dof) 1070 1071 # Calculated the standard deviation. 1072 std_red_chi2 = sqrt(red_chi2) 1073 1074 # Replace values with the stored value. 1075 for id in errors: 1076 errors[id] = std_red_chi2 1077 1078 return errors
1079 1080
1081 - def return_value(self, spin, param, sim=None, bc=False):
1082 """Return the value and error corresponding to the parameter. 1083 1084 If sim is set to an integer, return the value of the simulation and None. 1085 1086 1087 @param spin: The SpinContainer object. 1088 @type spin: SpinContainer 1089 @param param: The name of the parameter to return values for. 1090 @type param: str 1091 @keyword sim: The Monte Carlo simulation index. 1092 @type sim: None or int 1093 @keyword bc: The back-calculated data flag. If True, then the back-calculated data will be returned rather than the actual data. 1094 @type bc: bool 1095 @return: The value and error corresponding to 1096 @rtype: tuple of length 2 of floats or None 1097 """ 1098 1099 # Define list of special parameters. 1100 special_parameters = ['theta', 'w_eff'] 1101 1102 # Use api_common function for paramets not defined in special_parameters. 1103 if param not in special_parameters: 1104 returnval = self._return_value_general(spin=spin, param=param, sim=sim, bc=bc) 1105 return returnval 1106 1107 # If parameter in special_parameters, the do the following. 1108 1109 # Initial values. 1110 value = None 1111 error = None 1112 1113 # Return the data 1114 theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=spin) 1115 1116 # Return for parameter theta 1117 if param == "theta": 1118 value = theta_spin_dic 1119 1120 # Return for parameter theta 1121 elif param == "w_eff": 1122 value = w_eff_spin_dic 1123 1124 # Return the data. 1125 return value, error
1126 1127
1128 - def set_error(self, index, error, model_info=None):
1129 """Set the parameter errors. 1130 1131 @param index: The index of the parameter to set the errors for. 1132 @type index: int 1133 @param error: The error value. 1134 @type error: float 1135 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 1136 @type model_info: list of SpinContainer instances, list of str 1137 """ 1138 1139 # Unpack the data. 1140 spin_ids = model_info 1141 spins = spin_ids_to_containers(spin_ids) 1142 1143 # The number of parameters. 1144 total_param_num = param_num(spins=spins) 1145 1146 # No more model parameters. 1147 model_param = True 1148 if index >= total_param_num: 1149 model_param = False 1150 1151 # The auxiliary cluster parameters. 1152 aux_params = [] 1153 if 'pA' in spins[0].params: 1154 aux_params.append('pB') 1155 if 'pB' in spins[0].params: 1156 aux_params.append('pA') 1157 if 'kex' in spins[0].params: 1158 aux_params.append('tex') 1159 if 'tex' in spins[0].params: 1160 aux_params.append('kex') 1161 1162 # Convert the parameter index. 1163 if model_param: 1164 param_name, si, mi = param_index_to_param_info(index=index, spins=spins) 1165 else: 1166 param_name = aux_params[index - total_param_num] 1167 si = None 1168 mi = None 1169 1170 # The parameter error name. 1171 err_name = param_name + "_err" 1172 1173 # The exponential curve parameters. 1174 if param_name in ['r2eff', 'i0']: 1175 # Initialise if needed. 1176 if not hasattr(spins[si], err_name): 1177 setattr(spins[si], err_name, {}) 1178 1179 # Set the value. 1180 setattr(spins[si], err_name, error) 1181 1182 # Model and auxiliary parameters. 1183 else: 1184 # If clustered paramater: 1185 if si == None: 1186 for spin in spins: 1187 setattr(spin, err_name, error) 1188 1189 # If independent value. 1190 else: 1191 # Set the value. 1192 setattr(spins[si], err_name, error)
1193 1194
1195 - def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
1196 """Set the spin specific parameter values. 1197 1198 @keyword param: The parameter name list. 1199 @type param: list of str 1200 @keyword value: The parameter value list. 1201 @type value: list 1202 @keyword index: The index for parameters which are of the list-type. 1203 @type index: None or int 1204 @keyword spin_id: The spin identification string, only used for spin specific parameters. 1205 @type spin_id: None or str 1206 @keyword error: A flag which if True will allow the parameter errors to be set instead of the values. 1207 @type error: bool 1208 @keyword force: A flag which if True will cause current values to be overwritten. If False, a RelaxError will raised if the parameter value is already set. 1209 @type force: bool 1210 """ 1211 1212 # Checks. 1213 is_str_list(param, 'parameter name') 1214 is_list(value, 'parameter value') 1215 1216 # Loop over the spin blocks. 1217 model_index = -1 1218 for spin_ids in self.model_loop(): 1219 # Increment the model index. 1220 model_index += 1 1221 1222 # The spin containers. 1223 spins = spin_ids_to_containers(spin_ids) 1224 1225 # Skip deselected clusters. 1226 skip = True 1227 for spin in spins: 1228 if spin.select: 1229 skip = False 1230 if skip: 1231 continue 1232 1233 # Loop over the parameters. 1234 for i in range(len(param)): 1235 param_i = param[i] 1236 value_i = value[i] 1237 1238 # Is the parameter is valid? 1239 if not self._PARAMS.contains(param_i): 1240 raise RelaxError("The parameter '%s' is not valid for this data pipe type." % param_i) 1241 1242 # If the parameter is a global parameter, then change for all spins part of the cluster. 1243 if param_i in ['pA', 'pB', 'pC', 'kex', 'k_AB', 'k_BC', 'k_AC', 'tex', 'kB', 'kC', 'kex_AB', 'kex_BC', 'kex_AC']: 1244 selection_list = spin_ids 1245 else: 1246 selection_list = [spin_id] 1247 1248 # Now loop over selections in the list. 1249 for selection in selection_list: 1250 # Spin loop. 1251 for spin in spin_loop(selection=selection): 1252 # The object name. 1253 obj_name = param_i 1254 if error: 1255 obj_name += '_err' 1256 1257 # Handle the R10 parameters. 1258 if param_i in ['r1']: 1259 # Loop over the current keys. 1260 for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True): 1261 # The parameter key. 1262 key = generate_r20_key(exp_type=exp_type, frq=frq) 1263 1264 # Initialise the structure if needed. 1265 if not hasattr(spin, obj_name): 1266 setattr(spin, obj_name, {}) 1267 1268 # Set the value. 1269 if index == None: 1270 obj = getattr(spin, obj_name) 1271 obj[key] = value_i 1272 1273 # If the index is specified, let it match the frequency index 1274 elif mi == index: 1275 obj = getattr(spin, obj_name) 1276 obj[key] = value_i 1277 1278 # Handle the R20 parameters. 1279 elif param_i in PARAMS_R20: 1280 # Loop over the current keys. 1281 for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True): 1282 # The parameter key. 1283 key = generate_r20_key(exp_type=exp_type, frq=frq) 1284 1285 # Initialise the structure if needed. 1286 if not hasattr(spin, obj_name): 1287 setattr(spin, obj_name, {}) 1288 1289 # Set the value. 1290 if index == None: 1291 obj = getattr(spin, obj_name) 1292 obj[key] = value_i 1293 1294 # If the index is specified, let it match the frequency index 1295 elif mi == index: 1296 obj = getattr(spin, obj_name) 1297 obj[key] = value_i 1298 1299 # Handle the R2eff and I0 parameters. 1300 elif param_i in ['r2eff', 'i0'] and not isinstance(value_i, dict): 1301 # Loop over all the data. 1302 for exp_type, frq, offset, point in loop_exp_frq_offset_point(): 1303 # The parameter key. 1304 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 1305 1306 # Initialise the structure if needed. 1307 if not hasattr(spin, obj_name): 1308 setattr(spin, obj_name, {}) 1309 1310 # Set the value. 1311 obj = getattr(spin, obj_name) 1312 obj[key] = value_i 1313 1314 # Set the other parameters. 1315 else: 1316 setattr(spin, obj_name, value_i)
1317 1318
1319 - def set_selected_sim(self, select_sim, model_info=None):
1320 """Set the simulation selection flag. 1321 1322 @param select_sim: The selection flag for the simulations. 1323 @type select_sim: bool 1324 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 1325 @type model_info: list of SpinContainer instances, list of str 1326 """ 1327 1328 # Unpack the data. 1329 spin_ids = model_info 1330 spins = spin_ids_to_containers(spin_ids) 1331 1332 # Loop over the spins, storing the structure for each spin. 1333 for spin in spins: 1334 spin.select_sim = deepcopy(select_sim)
1335 1336
1337 - def sim_init_values(self):
1338 """Initialise the Monte Carlo parameter values.""" 1339 1340 # Get the parameter object names. 1341 param_names = self.data_names(set='params') 1342 1343 # Add the names of kex-tex pair. 1344 pairs = {} 1345 pairs['kex'] = 'tex' 1346 pairs['tex'] = 'kex' 1347 1348 # Add the names of pA-pB pair. 1349 pairs['pA'] = 'pB' 1350 pairs['pB'] = 'pA' 1351 1352 # Add the names of kex-k_AB pair and kex-k_BA pair. 1353 pairs['k_AB'] = 'kex' 1354 pairs['k_BA'] = 'kex' 1355 1356 # Get the minimisation statistic object names. 1357 min_names = self.data_names(set='min') 1358 1359 1360 # Test if Monte Carlo parameter values have already been set. 1361 ############################################################# 1362 1363 # Loop over the spins. 1364 for spin in spin_loop(): 1365 # Skip deselected spins. 1366 if not spin.select: 1367 continue 1368 1369 # Loop over all the parameter names. 1370 for object_name in param_names: 1371 # Name for the simulation object. 1372 sim_object_name = object_name + '_sim' 1373 1374 1375 # Set the Monte Carlo parameter values. 1376 ####################################### 1377 1378 # Loop over the spins. 1379 for spin in spin_loop(): 1380 # Skip deselected spins. 1381 if not spin.select: 1382 continue 1383 1384 # Skip protons for MMQ data. 1385 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 1386 continue 1387 1388 # Loop over all the data names. 1389 for object_name in param_names: 1390 # Not a parameter of the model. 1391 if not (object_name in spin.params or (object_name in pairs and pairs[object_name] in spin.params)): 1392 continue 1393 1394 # Name for the simulation object. 1395 sim_object_name = object_name + '_sim' 1396 1397 # Create the simulation object. 1398 setattr(spin, sim_object_name, []) 1399 1400 # Get the simulation object. 1401 sim_object = getattr(spin, sim_object_name) 1402 1403 # Loop over the simulations. 1404 for j in range(cdp.sim_number): 1405 # The non-simulation value. 1406 if object_name not in spin.params: 1407 value = deepcopy(getattr(spin, pairs[object_name])) 1408 else: 1409 value = deepcopy(getattr(spin, object_name)) 1410 1411 # Copy and append the data. 1412 sim_object.append(value) 1413 1414 # Loop over all the minimisation object names. 1415 for object_name in min_names: 1416 # Name for the simulation object. 1417 sim_object_name = object_name + '_sim' 1418 1419 # Create the simulation object. 1420 setattr(spin, sim_object_name, []) 1421 1422 # Get the simulation object. 1423 sim_object = getattr(spin, sim_object_name) 1424 1425 # Loop over the simulations. 1426 for j in range(cdp.sim_number): 1427 # Copy and append the data. 1428 sim_object.append(deepcopy(getattr(spin, object_name)))
1429 1430
1431 - def sim_pack_data(self, data_id, sim_data):
1432 """Pack the Monte Carlo simulation data. 1433 1434 @param data_id: The tuple of the spin container and the exponential curve identifying key, as yielded by the base_data_loop() generator method. 1435 @type data_id: SpinContainer instance and float 1436 @param sim_data: The Monte Carlo simulation data. 1437 @type sim_data: list of float 1438 """ 1439 1440 # The R2eff model (with peak intensity base data). 1441 if cdp.model_type == MODEL_R2EFF: 1442 # Unpack the data. 1443 spin, spin_id, exp_type, frq, offset, point = data_id 1444 1445 # Initialise the data structure if needed. 1446 if not hasattr(spin, 'peak_intensity_sim'): 1447 spin.peak_intensity_sim = {} 1448 1449 # Loop over each time point. 1450 ti = 0 1451 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 1452 # Get the intensity keys. 1453 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 1454 1455 # Loop over the intensity keys. 1456 for int_key in int_keys: 1457 # Test if the simulation data point already exists. 1458 if int_key in spin.peak_intensity_sim: 1459 raise RelaxError("Monte Carlo simulation data for the key '%s' already exists." % int_key) 1460 1461 # Initialise the list. 1462 spin.peak_intensity_sim[int_key] = [] 1463 1464 # Loop over the simulations, appending the corresponding data. 1465 for i in range(cdp.sim_number): 1466 spin.peak_intensity_sim[int_key].append(sim_data[i][ti]) 1467 1468 # Increment the time index. 1469 ti += 1 1470 1471 # All other models (with R2eff/R1rho base data). 1472 else: 1473 # Unpack the data. 1474 spin, spin_id = data_id 1475 1476 # 1H MMQ flag. 1477 proton_mmq_flag = has_proton_mmq_cpmg() 1478 1479 # Get the attached proton. 1480 proton = None 1481 if proton_mmq_flag: 1482 proton = return_attached_protons(spin_id)[0] 1483 1484 # Pack the data. 1485 spin.r2eff_sim = sim_data 1486 if proton != None: 1487 proton.r2eff_sim = sim_data
1488 1489
1490 - def sim_return_param(self, index, model_info=None):
1491 """Return the array of simulation parameter values. 1492 1493 @param index: The index of the parameter to return the array of values for. 1494 @type index: int 1495 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 1496 @type model_info: list of SpinContainer instances, list of str 1497 @return: The array of simulation parameter values. 1498 @rtype: list of float 1499 """ 1500 1501 # Unpack the data. 1502 spin_ids = model_info 1503 spins = spin_ids_to_containers(spin_ids) 1504 1505 # The number of parameters. 1506 total_param_num = param_num(spins=spins) 1507 1508 # No more model parameters. 1509 model_param = True 1510 if index >= total_param_num: 1511 model_param = False 1512 1513 # The auxiliary cluster parameters. 1514 aux_params = [] 1515 for spin in spins: 1516 if not spin.select: 1517 continue 1518 if 'pA' in spin.params: 1519 aux_params.append('pB') 1520 if 'pB' in spin.params: 1521 aux_params.append('pA') 1522 if 'kex' in spin.params: 1523 aux_params.append('tex') 1524 if 'tex' in spin.params: 1525 aux_params.append('kex') 1526 break 1527 1528 # No more auxiliary parameters. 1529 total_aux_num = total_param_num + len(aux_params) 1530 if index >= total_aux_num: 1531 return 1532 1533 # Convert the parameter index. 1534 if model_param: 1535 param_name, si, mi = param_index_to_param_info(index=index, spins=spins) 1536 if not param_name in ['r2eff', 'i0']: 1537 if si == None: 1538 si = 0 1539 1540 else: 1541 param_name = aux_params[index - total_param_num] 1542 si = 0 1543 mi = 0 1544 1545 # The exponential curve parameters. 1546 sim_data = [] 1547 if param_name == 'r2eff': 1548 for i in range(cdp.sim_number): 1549 sim_data.append(spins[si].r2eff_sim[i]) 1550 elif param_name == 'i0': 1551 for i in range(cdp.sim_number): 1552 sim_data.append(spins[si].i0_sim[i]) 1553 1554 # Model and auxiliary parameters. 1555 else: 1556 sim_data = getattr(spins[si], param_name + "_sim") 1557 1558 # Set the sim data to None if empty. 1559 if sim_data == []: 1560 sim_data = None 1561 1562 # Return the object. 1563 return sim_data
1564 1565
1566 - def sim_return_selected(self, model_info=None):
1567 """Return the array of selected simulation flags. 1568 1569 @keyword model_info: The list of spins and spin IDs per cluster originating from model_loop(). 1570 @type model_info: list of SpinContainer instances, list of str 1571 @return: The array of selected simulation flags. 1572 @rtype: list of int 1573 """ 1574 1575 # Unpack the data. 1576 spin_ids = model_info 1577 spins = spin_ids_to_containers(spin_ids) 1578 1579 # Return the array from the first spin, as this array will be identical for all spins in the cluster. 1580 return spins[0].select_sim
1581