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