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

Source Code for Module specific_analyses.relax_disp.api

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