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