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

Source Code for Module specific_analyses.relax_disp.disp_data

   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  """Module for handling relaxation dispersion data within the relax data store. 
  26   
  27  Ordering of data 
  28  ================ 
  29   
  30  The dispersion data model is based on the following concepts, in order of importance: 
  31   
  32      - 'exp', the experiment type, 
  33      - 'spin', the spins of the cluster, 
  34      - 'frq', the spectrometer frequency (if multiple field data is present), 
  35      - 'offset', the spin-lock offsets, 
  36      - 'point', the dispersion point (nu_CPMG value or spin-lock nu1 field strength), 
  37      - 'time', the relaxation time point (if exponential curve data has been collected). 
  38   
  39   
  40  Indices 
  41  ======= 
  42   
  43  The data structures used in this module consist of many different index types which follow the data ordering above.  These are abbreviated as: 
  44   
  45      - Ei or ei:  The index for each experiment type. 
  46      - Si or si:  The index for each spin of the spin cluster. 
  47      - Mi or mi:  The index for each magnetic field strength. 
  48      - Oi or oi:  The index for each spin-lock offset.  In the case of CPMG-type data, this index is always zero. 
  49      - Di or di:  The index for each dispersion point (either the spin-lock field strength or the nu_CPMG frequency). 
  50      - Ti or ti:  The index for each dispersion point (either the spin-lock field strength or the nu_CPMG frequency). 
  51  """ 
  52   
  53  # Python module imports. 
  54  from math import atan, floor, pi, sqrt 
  55  from numpy import array, float64, int32, ones, zeros 
  56  from os.path import expanduser 
  57  from random import gauss 
  58  from re import search 
  59  import sys 
  60  from warnings import warn 
  61   
  62  # relax module imports. 
  63  from lib.errors import RelaxError, RelaxNoSpectraError, RelaxNoSpinError, RelaxSpinTypeError 
  64  from lib.float import isNaN 
  65  from lib.io import extract_data, get_file_path, open_write_file, read_spin_data, strip, write_data, write_spin_data 
  66  from lib.nmr import frequency_to_ppm, frequency_to_rad_per_s 
  67  from lib.physical_constants import g1H, return_gyromagnetic_ratio 
  68  from lib.software.grace import write_xy_data, write_xy_header, script_grace2images 
  69  from lib.warnings import RelaxWarning, RelaxNoSpinWarning 
  70  from pipe_control import pipes 
  71  from pipe_control.mol_res_spin import check_mol_res_spin_data, exists_mol_res_spin_data, find_index, generate_spin_id_unique, get_spin_ids, return_spin, spin_loop 
  72  from pipe_control.result_files import add_result_file 
  73  from pipe_control.selection import desel_spin 
  74  from pipe_control.sequence import return_attached_protons 
  75  from pipe_control.spectrum import add_spectrum_id, get_ids 
  76  from pipe_control.spectrometer import check_frequency, get_frequency, set_frequency 
  77  import specific_analyses 
  78  from specific_analyses.relax_disp.checks import check_exp_type, check_mixed_curve_types 
  79  from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_DQ, EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_ZQ, EXP_TYPE_DESC_CPMG_DQ, EXP_TYPE_DESC_CPMG_MQ, EXP_TYPE_DESC_CPMG_PROTON_MQ, EXP_TYPE_DESC_CPMG_PROTON_SQ, EXP_TYPE_DESC_CPMG_SQ, EXP_TYPE_DESC_CPMG_ZQ, EXP_TYPE_DESC_R1RHO, EXP_TYPE_LIST, EXP_TYPE_LIST_CPMG, EXP_TYPE_LIST_R1RHO, EXP_TYPE_R1RHO, MODEL_DPL94, MODEL_LIST_MMQ, MODEL_LIST_NUMERIC_CPMG, MODEL_LIST_R1RHO_FULL, MODEL_MP05, MODEL_NS_R1RHO_2SITE, MODEL_R2EFF, MODEL_TAP03, MODEL_TP02 
  80  from stat import S_IRWXU, S_IRGRP, S_IROTH 
  81  from os import chmod, sep 
  82   
  83   
  84  # Module variables. 
  85  R20_KEY_FORMAT = "%s - %.8f MHz" 
  86   
  87   
88 -def average_intensity(spin=None, exp_type=None, frq=None, offset=None, point=None, time=None, sim_index=None, error=False):
89 """Return the average peak intensity for the spectrometer frequency, dispersion point, and relaxation time. 90 91 This is for handling replicate peak intensity data. 92 93 94 @keyword spin: The spin container to average the peak intensities for. 95 @type spin: SpinContainer instance 96 @keyword exp_type: The experiment type. 97 @type exp_type: str 98 @keyword frq: The spectrometer frequency. 99 @type frq: float 100 @keyword offset: The spin-lock or hard pulse offset. 101 @type offset: float 102 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 103 @type point: float 104 @keyword time: The relaxation time period. 105 @type time: float 106 @keyword sim_index: The simulation index. This should be None for the measured intensity values. 107 @type sim_index: None or int 108 @keyword error: A flag which if True will average and return the peak intensity errors. 109 @type error: bool 110 @return: The average peak intensity value. 111 @rtype: float 112 """ 113 114 # The keys. 115 int_keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 116 117 # Initialise. 118 intensity = 0.0 119 120 # Loop over the replicates. 121 for i in range(len(int_keys)): 122 # Simulation intensity data. 123 if sim_index != None: 124 # Error checking. 125 if not int_keys[i] in spin.intensity_sim: 126 raise RelaxError("The peak intensity simulation data is missing the key '%s'." % int_keys[i]) 127 128 # Sum. 129 intensity += spin.intensity_sim[int_keys[i]][sim_index] 130 131 # Error intensity data. 132 if error: 133 # Error checking. 134 if not int_keys[i] in spin.intensity_err: 135 raise RelaxError("The peak intensity errors are missing the key '%s'." % int_keys[i]) 136 137 # Sum. 138 intensity += spin.intensity_err[int_keys[i]]**2 139 140 # Normal intensity data. 141 else: 142 # Error checking. 143 if not int_keys[i] in spin.intensities: 144 raise RelaxError("The peak intensity data is missing the key '%s'." % int_keys[i]) 145 146 # Sum. 147 intensity += spin.intensities[int_keys[i]] 148 149 # Average. 150 if error: 151 intensity = sqrt(intensity / len(int_keys)) 152 else: 153 intensity /= len(int_keys) 154 155 # Return the value. 156 return intensity
157 158
159 -def calc_rotating_frame_params(spin=None, spin_id=None, fields=None, verbosity=0):
160 """Calculates and rotating frame parameters, calculated from: 161 - The spectrometer frequency. 162 - The spin-lock or hard pulse offset. 163 - The dispersion point data (the spin-lock field strength in Hz). 164 165 The return will be for each spin, 166 - Rotating frame tilt angle ( theta = arctan(w_1 / Omega) ) [rad] 167 - The average resonance offset in the rotating frame ( Domega = w_{pop_ave} - w_rf ) [rad/s] 168 - Effective field in rotating frame ( w_eff = sqrt( Omega^2 + w_1^2 ) ) [rad/s] 169 170 Calculations are mentioned in the U{manual<http://www.nmr-relax.com/manual/Dispersion_model_summary.html>} 171 172 @keyword spin: The spin system specific data container 173 @type spin: SpinContainer instance 174 @keyword spin_id: The spin ID string. 175 @type spin_id: None or str 176 @keyword fields: The spin-lock field strengths to use instead of the user loaded values - to enable interpolation. The dimensions are {Ei, Mi}. 177 @type fields: rank-2 list of floats 178 @keyword verbosity: A flag specifying to print calculations. 179 @type verbosity: int 180 @return: List with dict() of theta, Domega, w_eff and list of dict() keys. 181 @rtype: List of dict() 182 """ 183 184 # If the spin is not selected, return None 185 if not spin.select: 186 return None, None, None, None 187 188 # If the spin does not have isotope, return None 189 if not hasattr(spin, 'isotope'): 190 return None, None, None, None 191 192 # Get the field count 193 field_count = count_frq() 194 195 # Check the experiment type 196 if not has_r1rho_exp_type(): 197 raise RelaxError("The experiment type is not of R1rho type.") 198 199 # Get the spin_lock_field points 200 if fields == None: 201 spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False) 202 else: 203 spin_lock_nu1 = fields 204 205 # The offset and R1 data. 206 chemical_shifts, offsets, tilt_angles, Delta_omega, w_eff = return_offset_data(spins=[spin], spin_ids=[spin_id], field_count=field_count, fields=spin_lock_nu1) 207 208 # Loop over the index of spins, then exp_type, frq, offset 209 if verbosity: 210 print("Printing the following") 211 print("exp_type spin_id frq offset{ppm} offsets[ei][si][mi][oi]{rad/s} ei mi oi si di cur_spin.chemical_shift{ppm} chemical_shifts[ei][si][mi]{rad/s} spin_lock_nu1{Hz} tilt_angles[ei][si][mi][oi]{rad} av_res_offset[ei][si][mi][oi]{rad/s}") 212 213 si = 0 214 theta_spin_dic = dict() 215 Domega_spin_dic = dict() 216 w_eff_spin_dic = dict() 217 dic_key_list = [] 218 219 for exp_type, frq, offset, ei, mi, oi in loop_exp_frq_offset(return_indices=True): 220 # Loop over the dispersion points. 221 spin_lock_fields = spin_lock_nu1[ei][mi][oi] 222 for di in range(len(spin_lock_fields)): 223 if verbosity: 224 print("%-8s %-10s %11.1f %8.4f %12.5f %i %i %i %i %i %7.3f %12.5f %12.5f %12.5f %12.5f"%(exp_type, spin_id, frq, offset, offsets[ei][si][mi][oi], ei, mi, oi, si, di, spin.chemical_shift, chemical_shifts[ei][si][mi], spin_lock_fields[di], tilt_angles[ei][si][mi][oi][di], Delta_omega[ei][si][mi][oi][di])) 225 dic_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=spin_lock_fields[di]) 226 dic_key_list.append(dic_key) 227 theta_spin_dic["%s"%(dic_key)] = tilt_angles[ei][si][mi][oi][di] 228 Domega_spin_dic["%s"%(dic_key)] = Delta_omega[ei][si][mi][oi][di] 229 w_eff_spin_dic["%s"%(dic_key)] = w_eff[ei][si][mi][oi][di] 230 231 # Return the dic and list of keys 232 return [theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list]
233 234
235 -def count_exp():
236 """Count the number of experiments present. 237 238 @return: The experiment count 239 @rtype: int 240 """ 241 242 # The normal count variable. 243 return len(cdp.exp_type_list)
244 245
246 -def count_frq():
247 """Count the number of spectrometer frequencies present. 248 249 @return: The spectrometer frequency count 250 @rtype: int 251 """ 252 253 # Handle missing frequency data. 254 if not hasattr(cdp, 'spectrometer_frq'): 255 return 1 256 257 # The normal count variable. 258 return cdp.spectrometer_frq_count
259 260
261 -def count_relax_times(exp_type=None, frq=None, offset=None, point=None, ei=None):
262 """Count the number of relaxation times present. 263 264 @keyword exp_type: The experiment type. 265 @type exp_type: str 266 @keyword frq: The spectrometer frequency in Hz. 267 @type frq: float 268 @keyword offset: The spin-lock or hard pulse offset value in ppm. 269 @type offset: None or float 270 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 271 @type point: float 272 @keyword ei: The experiment type index. 273 @type ei: str 274 @return: The relaxation time count for the given experiment. 275 @rtype: int 276 """ 277 278 # Loop over the times. 279 count = 0 280 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 281 # Find a matching experiment ID. 282 found = False 283 for id in cdp.exp_type.keys(): 284 # Skip non-matching experiments. 285 if cdp.exp_type[id] != cdp.exp_type_list[ei]: 286 continue 287 288 # Found. 289 found = True 290 break 291 292 # No data. 293 if not found: 294 continue 295 296 # A new time. 297 count += 1 298 299 # Return the count. 300 return count
301 302
303 -def count_spins(spins=None):
304 """Count the number of selected spins in the spin cluster.""" 305 306 # Count the selected spins. 307 spin_num = 0 308 for spin in spins: 309 if spin.select: 310 spin_num += 1 311 312 # Return the count. 313 return spin_num
314 315
316 -def cpmg_frq(spectrum_id=None, cpmg_frq=None):
317 """Set the CPMG frequency associated with a given spectrum. 318 319 @keyword spectrum_id: The spectrum identification string. 320 @type spectrum_id: str 321 @keyword cpmg_frq: The frequency, in Hz, of the CPMG pulse train. 322 @type cpmg_frq: float 323 """ 324 325 # Test if the spectrum id exists. 326 if spectrum_id not in cdp.spectrum_ids: 327 raise RelaxNoSpectraError(spectrum_id) 328 329 # Initialise the global CPMG frequency data structures if needed. 330 if not hasattr(cdp, 'cpmg_frqs'): 331 cdp.cpmg_frqs = {} 332 if not hasattr(cdp, 'cpmg_frqs_list'): 333 cdp.cpmg_frqs_list = [] 334 335 # Add the frequency at the correct position, converting to a float if needed. 336 if cpmg_frq == None: 337 cdp.cpmg_frqs[spectrum_id] = cpmg_frq 338 else: 339 cdp.cpmg_frqs[spectrum_id] = float(cpmg_frq) 340 341 # The unique curves for the R2eff fitting (CPMG). 342 if cdp.cpmg_frqs[spectrum_id] not in cdp.cpmg_frqs_list: 343 cdp.cpmg_frqs_list.append(cdp.cpmg_frqs[spectrum_id]) 344 345 # Sort the list (handling None for Python 3). 346 flag = False 347 if None in cdp.cpmg_frqs_list: 348 cdp.cpmg_frqs_list.pop(cdp.cpmg_frqs_list.index(None)) 349 flag = True 350 cdp.cpmg_frqs_list.sort() 351 if flag: 352 cdp.cpmg_frqs_list.insert(0, None) 353 354 # Update the exponential curve count (skipping the reference if present). 355 cdp.dispersion_points = len(cdp.cpmg_frqs_list) 356 if None in cdp.cpmg_frqs_list: 357 cdp.dispersion_points -= 1 358 359 # Printout. 360 print("The spectrum ID '%s' CPMG frequency is set to %s Hz." % (spectrum_id, cdp.cpmg_frqs[spectrum_id]))
361 362
363 -def decompose_r20_key(key=None):
364 """Decompose the unique R20 key into the experiment type and spectrometer frequency. 365 366 @keyword key: The unique R20 key. 367 @type key: str 368 @return: The experiment and the spectrometer frequency in Hz. 369 @rtype: str, float 370 """ 371 372 # Loop over the experiments and frequencies until the matching key is found. 373 for exp_type, frq in loop_exp_frq(): 374 if key == generate_r20_key(exp_type=exp_type, frq=frq): 375 return exp_type, frq
376 377
378 -def find_intensity_keys(exp_type=None, frq=None, offset=None, point=None, time=None, raise_error=True):
379 """Return the key corresponding to the spectrometer frequency, dispersion point, and relaxation time. 380 381 @keyword exp_type: The experiment type. 382 @type exp_type: str 383 @keyword frq: The spectrometer frequency. 384 @type frq: float 385 @keyword offset: The optional offset value for off-resonance R1rho-type data. 386 @type offset: None or float 387 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 388 @type point: float 389 @keyword time: The relaxation time period. 390 @type time: float 391 @keyword raise_error: A flag which if True will cause a RelaxError to be raised if no keys could be found. 392 @type raise_error: bool 393 @return: The keys corresponding to the spectrometer frequency, dispersion point, and relaxation time. 394 @rtype: list of str 395 """ 396 397 # Check. 398 if exp_type == None: 399 raise RelaxError("The experiment type has not been supplied.") 400 401 # Catch NaNs. 402 if isNaN(point): 403 point = None 404 405 # The dispersion data. 406 if exp_type in EXP_TYPE_LIST_CPMG: 407 disp_data = cdp.cpmg_frqs 408 else: 409 disp_data = cdp.spin_lock_nu1 410 411 # Loop over all spectrum IDs, returning the matching ID. 412 ids = [] 413 for id in cdp.exp_type.keys(): 414 # Skip non-matching experiments. 415 if cdp.exp_type[id] != exp_type: 416 continue 417 418 # Skip non-matching spectrometer frequencies. 419 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 420 continue 421 422 # Skip non-matching offsets. 423 if offset != None and hasattr(cdp, 'spin_lock_offset') and cdp.spin_lock_offset[id] != offset: 424 continue 425 426 # Skip non-matching dispersion points. 427 if disp_data[id] != point: 428 continue 429 430 # The reference point, so checking the time is pointless (and can fail as specifying the time should not be necessary). 431 if point == None or isNaN(point): 432 ids.append(id) 433 434 # Matching time. 435 elif time == None: 436 ids.append(id) 437 elif cdp.relax_times[id] == time: 438 ids.append(id) 439 440 # Check for missing IDs. 441 if raise_error and len(ids) == 0: 442 if point == None or isNaN(point): 443 raise RelaxError("No reference intensity data could be found corresponding to the spectrometer frequency of %s MHz and relaxation time of %s s." % (frq*1e-6, time)) 444 else: 445 raise RelaxError("No intensity data could be found corresponding to the spectrometer frequency of %s MHz, dispersion point of %s and relaxation time of %s s." % (frq*1e-6, point, time)) 446 447 # Return the IDs. 448 return ids
449 450
451 -def generate_r20_key(exp_type=None, frq=None):
452 """Generate the unique R20 key from the experiment type and spectrometer frequency. 453 454 @keyword exp_type: The experiment type. 455 @type exp_type: str 456 @keyword frq: The spectrometer frequency in Hz. 457 @type frq: float 458 @return: The unique R20 key. 459 @rtype: str 460 """ 461 462 # Generate and return the unique key. 463 return R20_KEY_FORMAT % (exp_type, frq/1e6)
464 465
466 -def get_curve_type(id=None):
467 """Return the unique curve type. 468 469 @keyword id: The spectrum ID. If not supplied, then all data will be assumed. 470 @type id: str 471 @return: The curve type - either 'fixed time' or 'exponential'. 472 @rtype: str 473 """ 474 475 # All data. 476 if id == None: 477 # Data checks. 478 check_mixed_curve_types() 479 480 # Determine the curve type. 481 curve_type = 'fixed time' 482 if has_exponential_exp_type(): 483 curve_type = 'exponential' 484 485 # A specific ID. 486 else: 487 # Determine the curve type. 488 curve_type = 'exponential' 489 exp_type = cdp.exp_type[id] 490 frq = cdp.spectrometer_frq[id] 491 if count_relax_times(exp_type = exp_type, frq = frq, ei = cdp.exp_type_list.index(cdp.exp_type[id])) == 1: 492 curve_type = 'fixed time' 493 494 # Return the type. 495 return curve_type
496 497
498 -def get_exp_type(id=None):
499 """Return the experiment type for the given ID. 500 501 @keyword id: The spectrum ID. 502 @type id: str 503 @return: The experiment type corresponding to the ID. 504 @rtype: str 505 """ 506 507 # Data check. 508 check_exp_type(id=id) 509 510 # Return the type. 511 return cdp.exp_type[id]
512 513
514 -def has_cpmg_exp_type():
515 """Determine if the current data pipe contains CPMG experiment types. 516 517 @return: True if CPMG experiment types exist, False otherwise. 518 @rtype: bool 519 """ 520 521 # No experiment types set. 522 if not hasattr(cdp, 'exp_type'): 523 return False 524 525 # Loop over all experiment types. 526 for exp_type in cdp.exp_type_list: 527 if exp_type in EXP_TYPE_LIST_CPMG: 528 return True 529 530 # No CPMG experiment types. 531 return False
532 533
534 -def has_disp_data(spins=None, spin_ids=None, exp_type=None, frq=None, offset=None, point=None):
535 """Determine if dispersion data exists for the given data combination. 536 537 @keyword spins: The list of spin containers in the cluster. 538 @type spins: list of SpinContainer instances 539 @keyword spin_ids: The list of spin IDs for the cluster. 540 @type spin_ids: list of str 541 @keyword exp_type: The experiment type. 542 @type exp_type: str 543 @keyword frq: The spectrometer frequency. 544 @type frq: float 545 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm. 546 @type offset: None or float 547 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 548 @type point: float 549 @return: True if dispersion data exists, False otherwise. 550 @rtype: bool 551 """ 552 553 # Skip reference spectra. 554 if point == None: 555 return False 556 557 # The key. 558 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 559 560 # Loop over the spins. 561 for si in range(len(spins)): 562 # Alias the correct spin. 563 current_spin = spins[si] 564 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 565 current_spin = return_attached_protons(spin_ids[si])[0] 566 567 # The data is present. 568 if key in current_spin.r2eff.keys(): 569 return True 570 571 # No data. 572 return False
573 574
575 -def has_exponential_exp_type():
576 """Determine if the current data pipe contains exponential curves. 577 578 @return: True if spectral data for exponential curves exist, False otherwise. 579 @rtype: bool 580 """ 581 582 # No experiment types set. 583 if not hasattr(cdp, 'exp_type'): 584 return False 585 586 # Loop over all spectra IDs. 587 for id in cdp.exp_type.keys(): 588 if get_curve_type(id) == 'exponential': 589 return True 590 591 # No exponential data. 592 return False
593 594
595 -def has_fixed_time_exp_type():
596 """Determine if the current data pipe contains fixed time data. 597 598 @return: True if spectral data for fixed time data exists, False otherwise. 599 @rtype: bool 600 """ 601 602 # No experiment types set. 603 if not hasattr(cdp, 'exp_type'): 604 return False 605 606 # Loop over all experiment types. 607 for id in cdp.exp_type.keys(): 608 if get_curve_type(id) == 'fixed time': 609 return True 610 611 # No exponential data. 612 return False
613 614
615 -def has_proton_mmq_cpmg():
616 """Determine if the current data pipe contains either proton SQ or MQ (MMQ) CPMG data. 617 618 This is only for the MMQ models. 619 620 621 @return: True if either proton SQ or MQ CPMG data exists, False otherwise. 622 @rtype: bool 623 """ 624 625 # 1H MMQ data exists. 626 if has_proton_sq_cpmg(): 627 return True 628 if has_proton_mq_cpmg(): 629 return True 630 631 # No 1H MMQ CPMG data. 632 return False
633 634
635 -def has_proton_mq_cpmg():
636 """Determine if the current data pipe contains proton MQ CPMG data. 637 638 This is only for the MMQ models. 639 640 641 @return: True if proton MQ CPMG data exists, False otherwise. 642 @rtype: bool 643 """ 644 645 # Proton MQ CPMG data is present. 646 if EXP_TYPE_CPMG_PROTON_MQ in cdp.exp_type_list: 647 return True 648 649 # No 1H MQ CPMG data. 650 return False
651 652
653 -def has_proton_sq_cpmg():
654 """Determine if the current data pipe contains proton SQ CPMG data. 655 656 This is only for the MMQ models. 657 658 659 @return: True if proton SQ CPMG data exists, False otherwise. 660 @rtype: bool 661 """ 662 663 # Proton SQ CPMG data is present. 664 if EXP_TYPE_CPMG_PROTON_SQ in cdp.exp_type_list: 665 return True 666 667 # No 1H SQ CPMG data. 668 return False
669 670
671 -def has_r1rho_exp_type():
672 """Determine if the current data pipe contains R1rho experiment types. 673 674 @return: True if R1rho experiment types exist, False otherwise. 675 @rtype: bool 676 """ 677 678 # No experiment types set. 679 if not hasattr(cdp, 'exp_type'): 680 return False 681 682 # Loop over all experiment types. 683 for exp_type in cdp.exp_type_list: 684 if exp_type in EXP_TYPE_LIST_R1RHO: 685 return True 686 687 # No CPMG experiment types. 688 return False
689 690
691 -def insignificance(level=0.0):
692 """Deselect all spins with insignificant dispersion profiles. 693 694 @keyword level: The R2eff/R1rho value in rad/s by which to judge insignificance. If the maximum difference between two points on all dispersion curves for a spin is less than this value, that spin will be deselected. 695 @type level: float 696 """ 697 698 # Nothing to do. 699 if level == 0.0: 700 return 701 702 # Number of spectrometer fields. 703 fields = [None] 704 field_count = 1 705 if hasattr(cdp, 'spectrometer_frq_count'): 706 fields = cdp.spectrometer_frq_list 707 field_count = cdp.spectrometer_frq_count 708 709 # Loop over all spins. 710 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 711 # Nothing to do (the R2eff model has no dispersion curves). 712 if spin.model == 'R2eff': 713 continue 714 715 # Get all the data. 716 try: 717 values, errors, missing, frqs, frqs_H, exp_types, relax_times = return_r2eff_arrays(spins=[spin], spin_ids=[spin_id], fields=fields, field_count=field_count) 718 719 # No R2eff data, so skip the rest. 720 except RelaxError: 721 continue 722 723 # The flag. 724 desel = True 725 726 # Loop over the experiments, magnetic fields, and offsets. 727 max_diff = 0.0 728 for exp_type, frq, offset, ei, mi, oi in loop_exp_frq_offset(return_indices=True): 729 # No data. 730 if not len(values[ei][0][mi][oi]): 731 continue 732 733 # The difference. 734 diff = values[ei][0][mi][oi].max() - values[ei][0][mi][oi].min() 735 736 # Significance detected. 737 if diff > level: 738 desel = False 739 740 # Store the maximum for the deselection printout. 741 if diff > max_diff: 742 max_diff = diff 743 744 # Deselect the spin. 745 if desel: 746 # Printout. 747 print("Deselecting spin '%s', the maximum dispersion curve difference for all curves is %s rad/s." % (spin_id, max_diff)) 748 749 # Deselection. 750 desel_spin(spin_id)
751 752
753 -def is_cpmg_exp_type(id=None):
754 """Determine if the given spectrum ID corresponds to a CPMG experiment type. 755 756 @keyword id: The spectrum ID string. 757 @type id: str 758 @return: True if the spectrum ID corresponds to a CPMG experiment type, False otherwise. 759 @rtype: bool 760 """ 761 762 # No experiment type set. 763 if not hasattr(cdp, 'exp_type') or id not in cdp.exp_type: 764 return False 765 766 # CPMG experiment type. 767 if cdp.exp_type[id] in EXP_TYPE_LIST_CPMG: 768 return True 769 770 # Not a CPMG experiment type. 771 return False
772 773
774 -def is_r1rho_exp_type(id=None):
775 """Determine if the given spectrum ID corresponds to a R1rho experiment type. 776 777 @keyword id: The spectrum ID string. 778 @type id: str 779 @return: True if the spectrum ID corresponds to a R1rho experiment type, False otherwise. 780 @rtype: bool 781 """ 782 783 # No experiment type set. 784 if not hasattr(cdp, 'exp_type') or id not in cdp.exp_type: 785 return False 786 787 # R1rho experiment type. 788 if cdp.exp_type[id] in EXP_TYPE_LIST_R1RHO: 789 return True 790 791 # Not a R1rho experiment type. 792 return False
793 794
795 -def loop_cluster(skip_desel=True):
796 """Loop over the spin groupings for one model applied to multiple spins. 797 798 @keyword skip_desel: A flag which if True will cause deselected spins or spin clusters to be skipped. 799 @type skip_desel: bool 800 @return: The list of spin IDs per block will be yielded. 801 @rtype: list of str 802 """ 803 804 # No clustering, so loop over the sequence. 805 if not hasattr(cdp, 'clustering'): 806 for spin, spin_id in spin_loop(return_id=True, skip_desel=skip_desel): 807 # Skip protons for MMQ data. 808 if hasattr(spin, 'model') and spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 809 continue 810 811 # Return the spin ID as a list. 812 yield [spin_id] 813 814 # Loop over the clustering. 815 else: 816 # The clusters. 817 for key in cdp.clustering.keys(): 818 # Skip the free spins. 819 if key == 'free spins': 820 continue 821 822 # Create the spin ID lists. 823 spin_id_list = [] 824 for spin_id in cdp.clustering[key]: 825 # Skip deselected spins. 826 spin = return_spin(spin_id) 827 if skip_desel and not spin.select: 828 continue 829 830 # Skip protons for MMQ data. 831 if hasattr(spin, 'model') and spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 832 continue 833 834 # Add the spin ID. 835 spin_id_list.append(spin_id) 836 837 # Yield the cluster. 838 yield spin_id_list 839 840 # The free spins. 841 for spin_id in cdp.clustering['free spins']: 842 # Skip deselected spins. 843 spin = return_spin(spin_id) 844 if skip_desel and not spin.select: 845 continue 846 847 # Skip protons for MMQ data. 848 if hasattr(spin, 'model') and spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 849 continue 850 851 # Yield each spin individually. 852 yield [spin_id]
853 854
855 -def loop_exp(return_indices=False):
856 """Generator method for looping over all experiment types. 857 858 @keyword return_indices: A flag which if True will cause the experiment type index to be returned as well. 859 @type return_indices: bool 860 @return: The experiment type, and the index if asked. 861 @rtype: str, (int) 862 """ 863 864 # Initialise the index. 865 ei = -1 866 867 # Loop over the experiment types. 868 for exp_type in cdp.exp_type_list: 869 # Increment the index. 870 ei += 1 871 872 # Yield each unique experiment type. 873 if return_indices: 874 yield exp_type, ei 875 else: 876 yield exp_type
877 878
879 -def loop_exp_frq(return_indices=False):
880 """Generator method for looping over the exp and frq data. 881 882 These are the experiment types and spectrometer frequencies. 883 884 885 @keyword return_indices: A flag which if True will cause the experiment type and spectrometer frequency indices to be returned as well. 886 @type return_indices: bool 887 @return: The experiment type and spectrometer frequency in Hz, and the indices if asked. 888 @rtype: str, float, (int, int) 889 """ 890 891 # First loop over the experiment types. 892 for exp_type, ei in loop_exp(return_indices=True): 893 # Then loop over the spectrometer frequencies. 894 for frq, mi in loop_frq(return_indices=True): 895 # Yield the data. 896 if return_indices: 897 yield exp_type, frq, ei, mi 898 else: 899 yield exp_type, frq
900 901
902 -def loop_exp_frq_offset(return_indices=False):
903 """Generator method for looping over the exp, frq, and offset data. 904 905 These are the experiment types, spectrometer frequencies and spin-lock offset data. 906 907 908 @keyword return_indices: A flag which if True will cause the experiment type, spectrometer frequency and spin-lock offset indices to be returned as well. 909 @type return_indices: bool 910 @return: The experiment type, spectrometer frequency in Hz and spin-lock offset data, and the indices if asked. 911 @rtype: str, float, float, (int, int, int) 912 """ 913 914 # First loop over the experiment types. 915 for exp_type, ei in loop_exp(return_indices=True): 916 # Then loop over the spectrometer frequencies. 917 for frq, mi in loop_frq(return_indices=True): 918 # And finally the offset data. 919 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 920 # Yield the data. 921 if return_indices: 922 yield exp_type, frq, offset, ei, mi, oi 923 else: 924 yield exp_type, frq, offset
925 926
927 -def loop_exp_frq_offset_point(return_indices=False):
928 """Generator method for looping over the exp, frq, offset, and point data. 929 930 These are the experiment types, spectrometer frequencies, spin-lock offset data, and dispersion points. 931 932 933 @keyword return_indices: A flag which if True will cause the experiment type, spectrometer frequency, spin-lock offset and dispersion point indices to be returned as well. 934 @type return_indices: bool 935 @return: The experiment type, spectrometer frequency in Hz, spin-lock offset data and dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), and the indices if asked. 936 @rtype: str, float, float, float, (int, int, int, int) 937 """ 938 939 # First loop over the experiment types. 940 for exp_type, ei in loop_exp(return_indices=True): 941 # Then loop over the spectrometer frequencies. 942 for frq, mi in loop_frq(return_indices=True): 943 # Then loop over the offset data. 944 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 945 # And finally the dispersion points. 946 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 947 # Yield the data. 948 if return_indices: 949 yield exp_type, frq, offset, point, ei, mi, oi, di 950 else: 951 yield exp_type, frq, offset, point
952 953
954 -def loop_exp_frq_offset_point_time(return_indices=False):
955 """Generator method for looping over the exp, frq, offset, and point data. 956 957 These are the experiment types, spectrometer frequencies, spin-lock offset data, and dispersion points. 958 959 960 @keyword return_indices: A flag which if True will cause the experiment type, spectrometer frequency, spin-lock offset and dispersion point indices to be returned as well. 961 @type return_indices: bool 962 @return: The experiment type, spectrometer frequency in Hz, spin-lock offset data and dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), and the indices if asked. 963 @rtype: str, float, float, float, (int, int, int, int) 964 """ 965 966 # First loop over the experiment types. 967 for exp_type, ei in loop_exp(return_indices=True): 968 # Then loop over the spectrometer frequencies. 969 for frq, mi in loop_frq(return_indices=True): 970 # Then loop over the offset data. 971 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 972 # Then the dispersion points. 973 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 974 # Finally the relaxation times. 975 for time, ti in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point, return_indices=True): 976 # Yield the data. 977 if return_indices: 978 yield exp_type, frq, offset, point, time, ei, mi, oi, di, ti 979 else: 980 yield exp_type, frq, offset, point, time
981 982
983 -def loop_exp_frq_point(return_indices=False):
984 """Generator method for looping over the exp, frq, and point data. 985 986 These are the experiment types, spectrometer frequencies and dispersion points. 987 988 989 @keyword return_indices: A flag which if True will cause the experiment type, spectrometer frequency and dispersion point indices to be returned as well. 990 @type return_indices: bool 991 @return: The experiment type, spectrometer frequency in Hz and dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), and the indices if asked. 992 @rtype: str, float, float, (int, int, int) 993 """ 994 995 # First loop over the experiment types. 996 for exp_type, ei in loop_exp(return_indices=True): 997 # Then loop over the spectrometer frequencies. 998 for frq, mi in loop_frq(return_indices=True): 999 # And finally the dispersion points. 1000 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=0.0, return_indices=True): 1001 # Yield the data. 1002 if return_indices: 1003 yield exp_type, frq, point, ei, mi, di 1004 else: 1005 yield exp_type, frq, point
1006 1007
1008 -def loop_exp_frq_point_time(return_indices=False):
1009 """Generator method for looping over the exp, frq, point, and time data. 1010 1011 These are the experiment types, spectrometer frequencies, dispersion points, and relaxation times. 1012 1013 1014 @keyword return_indices: A flag which if True will cause the experiment type, spectrometer frequency, dispersion point, and relaxation time indices to be returned as well. 1015 @type return_indices: bool 1016 @return: The experiment type, spectrometer frequency in Hz, dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), the relaxation time, and the indices if asked. 1017 @rtype: str, float, float, float, (int, int, int, int) 1018 """ 1019 1020 # First loop over the experiment types. 1021 for exp_type, ei in loop_exp(return_indices=True): 1022 # Then the spectrometer frequencies. 1023 for frq, mi in loop_frq(return_indices=True): 1024 # Then the dispersion points. 1025 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=0.0, return_indices=True): 1026 # Finally the relaxation times. 1027 for time, ti in loop_time(exp_type=exp_type, frq=frq, point=point, return_indices=True): 1028 # Yield all data. 1029 if return_indices: 1030 yield exp_type, frq, point, time, ei, mi, di, ti 1031 else: 1032 yield exp_type, frq, point, time
1033 1034
1035 -def loop_frq(return_indices=False):
1036 """Generator method for looping over all spectrometer frequencies. 1037 1038 @keyword return_indices: A flag which if True will cause the spectrometer frequency index to be returned as well. 1039 @type return_indices: bool 1040 @return: The spectrometer frequency in Hz, and the index if asked. 1041 @rtype: float, (int) 1042 """ 1043 1044 # Handle missing frequency data. 1045 frqs = [None] 1046 if hasattr(cdp, 'spectrometer_frq_list'): 1047 frqs = cdp.spectrometer_frq_list 1048 1049 # Initialise the index. 1050 mi = -1 1051 1052 # Loop over the spectrometer frequencies. 1053 for field in frqs: 1054 # Increment the index. 1055 mi += 1 1056 1057 # Yield each unique spectrometer field strength. 1058 if return_indices: 1059 yield field, mi 1060 else: 1061 yield field
1062 1063
1064 -def loop_frq_offset(exp_type=None, return_indices=False):
1065 """Generator method for looping over the spectrometer frequencies and dispersion points. 1066 1067 @keyword exp_type: The experiment type. 1068 @type exp_type: str 1069 @keyword return_indices: A flag which if True will cause the spectrometer frequency and dispersion point indices to be returned as well. 1070 @type return_indices: bool 1071 @return: The spectrometer frequency in Hz and dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 1072 @rtype: float, float, (int, int) 1073 """ 1074 1075 # Checks. 1076 if exp_type == None: 1077 raise RelaxError("The experiment type must be supplied.") 1078 1079 # First loop over the spectrometer frequencies. 1080 for frq, mi in loop_frq(return_indices=True): 1081 # Then the offset points. 1082 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 1083 # Yield the data. 1084 if return_indices: 1085 yield frq, offset, mi, oi 1086 else: 1087 yield frq, offset
1088 1089
1090 -def loop_frq_point(exp_type=None, return_indices=False):
1091 """Generator method for looping over the spectrometer frequencies and dispersion points. 1092 1093 @keyword exp_type: The experiment type. 1094 @type exp_type: str 1095 @keyword return_indices: A flag which if True will cause the spectrometer frequency and dispersion point indices to be returned as well. 1096 @type return_indices: bool 1097 @return: The spectrometer frequency in Hz and dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 1098 @rtype: float, float, (int, int) 1099 """ 1100 1101 # First loop over the spectrometer frequencies. 1102 for frq, mi in loop_frq(return_indices=True): 1103 # Then the dispersion points. 1104 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=0.0, return_indices=True): 1105 # Yield the data. 1106 if return_indices: 1107 yield frq, point, mi, di 1108 else: 1109 yield frq, point
1110 1111
1112 -def loop_frq_offset_point_key(exp_type=None):
1113 """Generator method for looping over the spectrometer frequencies, spin-lock offsets and dispersion points (returning the key). 1114 1115 @keyword exp_type: The experiment type. 1116 @type exp_type: str 1117 @return: The key corresponding to the spectrometer frequency, offset and dispersion point. 1118 @rtype: str 1119 """ 1120 1121 # First loop over the spectrometer frequencies, offsets and dispersion points. 1122 for frq, offset, point in loop_frq_offset_point(return_indices=True): 1123 # Generate and yield the key. 1124 yield return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point)
1125 1126
1127 -def loop_frq_point_time(exp_type=None, return_indices=False):
1128 """Generator method for looping over the spectrometer frequencies, dispersion points, and relaxation times. 1129 1130 @keyword exp_type: The experiment type. 1131 @type exp_type: str 1132 @keyword return_indices: A flag which if True will cause the spectrometer frequency, dispersion point, and relaxation time indices to be returned as well. 1133 @type return_indices: bool 1134 @return: The spectrometer frequency in Hz, dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), and the relaxation time. 1135 @rtype: float, float, float 1136 """ 1137 1138 # First loop over the spectrometer frequencies. 1139 for frq, mi in loop_frq(return_indices=True): 1140 # Then the dispersion points. 1141 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=0.0, return_indices=True): 1142 # Finally the relaxation times. 1143 for time, ti in loop_time(exp_type=exp_type, frq=frq, point=point, return_indices=True): 1144 # Yield all data. 1145 if return_indices: 1146 yield frq, point, time, mi, di, ti 1147 else: 1148 yield frq, point, time
1149 1150
1151 -def loop_offset(exp_type=None, frq=None, return_indices=False):
1152 """Generator method for looping over the spin-lock offset values. 1153 1154 @keyword exp_type: The experiment type. 1155 @type exp_type: str 1156 @keyword frq: The spectrometer frequency. 1157 @type frq: float 1158 @keyword return_indices: A flag which if True will cause the offset index to be returned as well. 1159 @type return_indices: bool 1160 @return: The spin-lock offset value and the index if asked. 1161 @rtype: float, (int) 1162 """ 1163 1164 # Checks. 1165 if exp_type == None: 1166 raise RelaxError("The experiment type must be supplied.") 1167 if frq == None: 1168 raise RelaxError("The spectrometer frequency must be supplied.") 1169 1170 # Initialise the index. 1171 oi = -1 1172 1173 # CPMG-type data. 1174 if exp_type in EXP_TYPE_LIST_CPMG: 1175 # Yield a single set of dummy values until hard pulse offset handling is implemented. 1176 yield 0.0, 0 1177 1178 # R1rho-type data. 1179 if exp_type in EXP_TYPE_LIST_R1RHO: 1180 # No offsets set. 1181 if not hasattr(cdp, 'spin_lock_offset_list'): 1182 yield 0.0, 0 1183 1184 # Loop over the offset data. 1185 else: 1186 for offset in cdp.spin_lock_offset_list: 1187 # Find a matching experiment ID. 1188 found = False 1189 for id in cdp.exp_type.keys(): 1190 # Skip non-matching experiments. 1191 if cdp.exp_type[id] != exp_type: 1192 continue 1193 1194 # Skip non-matching spectrometer frequencies. 1195 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 1196 continue 1197 1198 # Skip non-matching offsets. 1199 if cdp.spin_lock_offset[id] != offset: 1200 continue 1201 1202 # Found. 1203 found = True 1204 break 1205 1206 # No data. 1207 if not found: 1208 continue 1209 1210 # Increment the index. 1211 oi += 1 1212 1213 # Yield each unique field strength or frequency. 1214 if return_indices: 1215 yield offset, oi 1216 else: 1217 yield offset
1218 1219
1220 -def loop_offset_point(exp_type=None, frq=None, skip_ref=True, return_indices=False):
1221 """Generator method for looping over the offsets and dispersion points. 1222 1223 @keyword exp_type: The experiment type. 1224 @type exp_type: str 1225 @keyword frq: The spectrometer frequency. 1226 @type frq: float 1227 @keyword skip_ref: A flag which if True will cause the reference point to be skipped. 1228 @type skip_ref: bool 1229 @keyword return_indices: A flag which if True will cause the offset and dispersion point indices to be returned as well. 1230 @type return_indices: bool 1231 @return: The offsets in ppm and the dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), and the index if asked. 1232 @rtype: float, float, (int, int) 1233 """ 1234 1235 # First loop over the offsets. 1236 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 1237 # Then the dispersion points. 1238 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 1239 # Yield all data. 1240 if return_indices: 1241 yield offset, point, oi, di 1242 else: 1243 yield offset, point
1244 1245
1246 -def loop_point(exp_type=None, frq=None, offset=None, time=None, skip_ref=True, return_indices=False):
1247 """Generator method for looping over the dispersion points. 1248 1249 @keyword exp_type: The experiment type. 1250 @type exp_type: str 1251 @keyword frq: The spectrometer frequency. 1252 @type frq: float 1253 @keyword offset: The spin-lock or hard pulse offset value in ppm. 1254 @type offset: None or float 1255 @keyword time: The relaxation time period. 1256 @type time: float 1257 @keyword skip_ref: A flag which if True will cause the reference point to be skipped. 1258 @type skip_ref: bool 1259 @keyword return_indices: A flag which if True will cause the experiment type index to be returned as well. 1260 @type return_indices: bool 1261 @return: Dispersion point data for the given indices (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz), and the index if asked. 1262 @rtype: float, (int) 1263 """ 1264 1265 # Checks. 1266 if exp_type == None: 1267 raise RelaxError("The experiment type must be supplied.") 1268 if frq == None: 1269 raise RelaxError("The spectrometer frequency must be supplied.") 1270 if offset == None: 1271 raise RelaxError("The offset must be supplied.") 1272 1273 # Assemble the dispersion data. 1274 ref_flag = not skip_ref 1275 if exp_type in EXP_TYPE_LIST_CPMG: 1276 fields = return_cpmg_frqs_single(exp_type=exp_type, frq=frq, offset=offset, time=time, ref_flag=ref_flag) 1277 elif exp_type in EXP_TYPE_LIST_R1RHO: 1278 fields = return_spin_lock_nu1_single(exp_type=exp_type, frq=frq, offset=offset, ref_flag=ref_flag) 1279 else: 1280 raise RelaxError("The experiment type '%s' is unknown." % exp_type) 1281 1282 # Initialise the index. 1283 di = -1 1284 1285 # Loop over the field data. 1286 for field in fields: 1287 # Skip the reference (the None value will be converted to the numpy nan value). 1288 if skip_ref and isNaN(field): 1289 continue 1290 1291 # Increment the index. 1292 di += 1 1293 1294 # Yield each unique field strength or frequency. 1295 if return_indices: 1296 yield field, di 1297 else: 1298 yield field
1299 1300
1301 -def loop_spectrum_ids(exp_type=None, frq=None, offset=None, point=None, time=None):
1302 """Generator method for selectively looping over the spectrum IDs. 1303 1304 @keyword exp_type: The experiment type. 1305 @type exp_type: str 1306 @keyword frq: The spectrometer frequency. 1307 @type frq: float 1308 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm. 1309 @type offset: None or float 1310 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 1311 @type point: float 1312 @keyword time: The relaxation time period. 1313 @type time: float 1314 @return: The spectrum ID. 1315 @rtype: str 1316 """ 1317 1318 # Loop over all spectrum IDs. 1319 for id in cdp.spectrum_ids: 1320 # Experiment type filter. 1321 if exp_type != None: 1322 # No experiment type set. 1323 if not hasattr(cdp, 'exp_type') or id not in cdp.exp_type: 1324 continue 1325 1326 # No match. 1327 if cdp.exp_type[id] != exp_type: 1328 continue 1329 1330 # The frequency filter. 1331 if frq != None: 1332 # No frequency data set. 1333 if not hasattr(cdp, 'spectrometer_frq') or id not in cdp.spectrometer_frq: 1334 continue 1335 1336 # No match. 1337 if cdp.spectrometer_frq[id] != frq: 1338 continue 1339 1340 # The dispersion point filter. 1341 if point != None: 1342 # No experiment type set. 1343 if not hasattr(cdp, 'exp_type') or id not in cdp.exp_type: 1344 continue 1345 1346 # The experiment type. 1347 exp_type = cdp.exp_type[id] 1348 1349 # The CPMG dispersion data. 1350 if exp_type in EXP_TYPE_LIST_CPMG: 1351 # No dispersion point data set. 1352 if not hasattr(cdp, 'cpmg_frqs') or id not in cdp.cpmg_frqs: 1353 continue 1354 1355 # Alias the structure 1356 disp_data = cdp.cpmg_frqs 1357 1358 # The R1rho dispersion data. 1359 else: 1360 # No dispersion point data set. 1361 if not hasattr(cdp, 'spin_lock_nu1') or id not in cdp.spin_lock_nu1: 1362 continue 1363 1364 # Alias the structure 1365 disp_data = cdp.spin_lock_nu1 1366 1367 # No match. 1368 if disp_data[id] != point: 1369 continue 1370 1371 # The time filter. 1372 if time != None: 1373 # No time data set. 1374 if not hasattr(cdp, 'relax_times') or id not in cdp.relax_times: 1375 continue 1376 1377 # No match. 1378 if cdp.relax_times[id] != time: 1379 continue 1380 1381 # Yield the Id. 1382 yield id
1383 1384
1385 -def loop_time(exp_type=None, frq=None, offset=None, point=None, return_indices=False):
1386 """Generator method for looping over the relaxation times. 1387 1388 @keyword exp_type: The experiment type. 1389 @type exp_type: str 1390 @keyword frq: The spectrometer frequency in Hz. 1391 @type frq: float 1392 @keyword offset: The spin-lock or hard pulse offset value in ppm. 1393 @type offset: None or float 1394 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 1395 @type point: float 1396 @keyword return_indices: A flag which if True will cause the relaxation time index to be returned as well. 1397 @type return_indices: bool 1398 @return: The relaxation time. 1399 @rtype: float 1400 """ 1401 1402 # Initialise the index. 1403 ti = -1 1404 1405 # Loop over the time points. 1406 if hasattr(cdp, 'relax_time_list'): 1407 for time in cdp.relax_time_list: 1408 # Find a matching experiment ID. 1409 found = False 1410 for id in cdp.exp_type.keys(): 1411 # Skip non-matching experiments. 1412 if exp_type != None and cdp.exp_type[id] != exp_type: 1413 continue 1414 1415 # Skip non-matching spectrometer frequencies. 1416 if frq != None and hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 1417 continue 1418 1419 # Skip non-matching offsets. 1420 if offset != None and hasattr(cdp, 'spin_lock_offset') and cdp.spin_lock_offset[id] != offset: 1421 continue 1422 1423 # The dispersion point filter. 1424 if point != None: 1425 # No experiment type set. 1426 if not hasattr(cdp, 'exp_type') or id not in cdp.exp_type: 1427 continue 1428 1429 # The experiment type. 1430 exp_type = cdp.exp_type[id] 1431 1432 # The CPMG dispersion data. 1433 if exp_type in EXP_TYPE_LIST_CPMG: 1434 # No dispersion point data set. 1435 if hasattr(cdp, 'cpmg_frqs') and cdp.cpmg_frqs[id] != point: 1436 continue 1437 1438 # The R1rho data 1439 if exp_type in EXP_TYPE_R1RHO: 1440 if hasattr(cdp, 'spin_lock_nu1') and cdp.spin_lock_nu1[id] != point: 1441 continue 1442 1443 if time != cdp.relax_times[id]: 1444 continue 1445 1446 # Found. 1447 found = True 1448 break 1449 1450 # No data. 1451 if not found: 1452 continue 1453 1454 # Increment the index. 1455 ti += 1 1456 1457 # Yield each unique relaxation time. 1458 if return_indices: 1459 yield time, ti 1460 else: 1461 yield time 1462 1463 # No times set. 1464 else: 1465 if return_indices: 1466 yield None, None 1467 else: 1468 yield None
1469 1470
1471 -def num_exp_types():
1472 """Count the number of experiment types present. 1473 1474 @return: The number of experiment types. 1475 @rtype: int 1476 """ 1477 1478 # The count. 1479 count = len(cdp.exp_type_list) 1480 1481 # Return the count. 1482 return count
1483 1484
1485 -def pack_back_calc_r2eff(spin=None, spin_id=None, si=None, back_calc=None, proton_mmq_flag=False):
1486 """Store the back calculated R2eff data for the given spin. 1487 1488 @keyword spin: The spin data container to store the data in. 1489 @type spin: SpinContainer instance 1490 @keyword spin_id: The spin ID string. 1491 @type spin_id: str 1492 @keyword si: The index of the given spin in the cluster. 1493 @type si: int 1494 @keyword back_calc: The back calculated data. The first index corresponds to the experiment type, the second is the spin of the cluster, the third is the magnetic field strength, and the fourth is the dispersion point. 1495 @type back_calc: list of lists of lists of lists of float 1496 @keyword proton_mmq_flag: The flag specifying if proton SQ or MQ CPMG data exists for the spin. 1497 @type proton_mmq_flag: bool 1498 """ 1499 1500 # Get the attached proton. 1501 proton = None 1502 if proton_mmq_flag: 1503 proton = return_attached_protons(spin_id)[0] 1504 1505 # Loop over the R2eff data. 1506 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 1507 # The R2eff key. 1508 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 1509 1510 # Alias the correct spin. 1511 current_spin = spin 1512 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 1513 current_spin = proton 1514 1515 # Missing data. 1516 if not hasattr(current_spin, 'r2eff') or key not in current_spin.r2eff.keys(): 1517 continue 1518 1519 # Initialise. 1520 if not hasattr(current_spin, 'r2eff_bc'): 1521 current_spin.r2eff_bc = {} 1522 1523 # Store the back-calculated data. 1524 current_spin.r2eff_bc[key] = back_calc[ei][si][mi][oi][di]
1525 1526
1527 -def plot_disp_curves(dir=None, num_points=1000, extend=500.0, force=False):
1528 """Custom 2D Grace plotting function for the dispersion curves. 1529 1530 One file will be created per spin system. 1531 1532 A python "grace to PNG/EPS/SVG..." conversion script is created at the end 1533 1534 @keyword dir: The optional directory to place the file into. 1535 @type dir: str 1536 @keyword num_points: The number of points to generate the interpolated fitted curves with. 1537 @type num_points: int 1538 @keyword extend: How far to extend the interpolated fitted curves to (in Hz). 1539 @type extend: float 1540 @param force: Boolean argument which if True causes the files to be overwritten if it already exists. 1541 @type force: bool 1542 """ 1543 1544 # Checks. 1545 pipes.test() 1546 check_mol_res_spin_data() 1547 1548 # 1H MMQ flag. 1549 proton_mmq_flag = has_proton_mmq_cpmg() 1550 1551 # Default hardcoded colours (one colour for each magnetic field strength). 1552 colour_order = [4, 15, 2, 13, 11, 1, 3, 5, 6, 7, 8, 9, 10, 12, 14] * 1000 1553 1554 # Loop over each spin. 1555 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 1556 # Skip protons for MMQ data. 1557 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 1558 continue 1559 1560 # Initialise some data structures. 1561 data = [] 1562 set_labels = [] 1563 x_err_flag = False 1564 y_err_flag = False 1565 axis_labels = [] 1566 set_colours = [] 1567 x_axis_type_zero = [] 1568 symbols = [] 1569 symbol_sizes = [] 1570 linetype = [] 1571 linestyle = [] 1572 1573 # The unique file name. 1574 file_name = "disp%s.agr" % spin_id.replace('#', '_').replace(':', '_').replace('@', '_') 1575 1576 # Open the file for writing. 1577 file_path = get_file_path(file_name, dir) 1578 file = open_write_file(file_name, dir, force) 1579 1580 # Get the attached proton. 1581 proton = None 1582 if proton_mmq_flag: 1583 proton = return_attached_protons(spin_id)[0] 1584 1585 # Set up the interpolated curve data structures. 1586 interpolated_flag = False 1587 if not spin.model in [MODEL_R2EFF]: 1588 # Set the flag. 1589 interpolated_flag = True 1590 1591 # Initialise some structures. 1592 cpmg_frqs_new = None 1593 spin_lock_nu1_new = None 1594 1595 # Interpolate the CPMG frequencies (numeric models). 1596 if spin.model in MODEL_LIST_NUMERIC_CPMG: 1597 cpmg_frqs = return_cpmg_frqs(ref_flag=False) 1598 relax_times = return_relax_times() 1599 if cpmg_frqs != None and len(cpmg_frqs[0][0]): 1600 cpmg_frqs_new = [] 1601 for ei in range(len(cpmg_frqs)): 1602 # Add a new dimension. 1603 cpmg_frqs_new.append([]) 1604 1605 # Then loop over the spectrometer frequencies. 1606 for mi in range(len(cpmg_frqs[ei])): 1607 # Add a new dimension. 1608 cpmg_frqs_new[ei].append([]) 1609 1610 # Finally the offsets. 1611 for oi in range(len(cpmg_frqs[ei][mi])): 1612 # Add a new dimension. 1613 cpmg_frqs_new[ei][mi].append([]) 1614 1615 # No data. 1616 if not len(cpmg_frqs[ei][mi][oi]): 1617 continue 1618 1619 # The minimum frequency unit. 1620 min_frq = 1.0 / relax_times[ei][mi] 1621 max_frq = max(cpmg_frqs[ei][mi][oi]) + round(extend / min_frq) * min_frq 1622 num_points = int(round(max_frq / min_frq)) 1623 1624 # Interpolate (adding the extended amount to the end). 1625 for di in range(num_points): 1626 point = (di + 1) * min_frq 1627 cpmg_frqs_new[ei][mi][oi].append(point) 1628 1629 # Convert to a numpy array. 1630 cpmg_frqs_new[ei][mi][oi] = array(cpmg_frqs_new[ei][mi][oi], float64) 1631 1632 # Interpolate the CPMG frequencies (analytic models). 1633 else: 1634 cpmg_frqs = return_cpmg_frqs(ref_flag=False) 1635 if cpmg_frqs != None and len(cpmg_frqs[0][0]): 1636 cpmg_frqs_new = [] 1637 for ei in range(len(cpmg_frqs)): 1638 # Add a new dimension. 1639 cpmg_frqs_new.append([]) 1640 1641 # Then loop over the spectrometer frequencies. 1642 for mi in range(len(cpmg_frqs[ei])): 1643 # Add a new dimension. 1644 cpmg_frqs_new[ei].append([]) 1645 1646 # Finally the offsets. 1647 for oi in range(len(cpmg_frqs[ei][mi])): 1648 # Add a new dimension. 1649 cpmg_frqs_new[ei][mi].append([]) 1650 1651 # No data. 1652 if not len(cpmg_frqs[ei][mi][oi]): 1653 continue 1654 1655 # Interpolate (adding the extended amount to the end). 1656 for di in range(num_points): 1657 point = (di + 1) * (max(cpmg_frqs[ei][mi][oi])+extend) / num_points 1658 cpmg_frqs_new[ei][mi][oi].append(point) 1659 1660 # Convert to a numpy array. 1661 cpmg_frqs_new[ei][mi][oi] = array(cpmg_frqs_new[ei][mi][oi], float64) 1662 1663 # Interpolate the spin-lock field strengths. 1664 spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False) 1665 if spin_lock_nu1 != None and len(spin_lock_nu1[0][0][0]): 1666 spin_lock_nu1_new = [] 1667 for ei in range(len(spin_lock_nu1)): 1668 # Add a new dimension. 1669 spin_lock_nu1_new.append([]) 1670 1671 # Then loop over the spectrometer frequencies. 1672 for mi in range(len(spin_lock_nu1[ei])): 1673 # Add a new dimension. 1674 spin_lock_nu1_new[ei].append([]) 1675 1676 # Finally the offsets. 1677 for oi in range(len(spin_lock_nu1[ei][mi])): 1678 # Add a new dimension. 1679 spin_lock_nu1_new[ei][mi].append([]) 1680 1681 # No data. 1682 if not len(spin_lock_nu1[ei][mi][oi]): 1683 continue 1684 1685 # Interpolate (adding the extended amount to the end). 1686 for di in range(num_points): 1687 point = (di + 1) * (max(spin_lock_nu1[ei][mi][oi])+extend) / num_points 1688 spin_lock_nu1_new[ei][mi][oi].append(point) 1689 1690 # Convert to a numpy array. 1691 spin_lock_nu1_new[ei][mi][oi] = array(spin_lock_nu1_new[ei][mi][oi], float64) 1692 1693 # Back calculate R2eff data for the second sets of plots. 1694 back_calc = specific_analyses.relax_disp.optimisation.back_calc_r2eff(spin=spin, spin_id=spin_id, cpmg_frqs=cpmg_frqs_new, spin_lock_nu1=spin_lock_nu1_new) 1695 1696 # Loop over each experiment type. 1697 graph_index = 0 1698 for exp_type, ei in loop_exp(return_indices=True): 1699 # Update the structures. 1700 data.append([]) 1701 set_labels.append([]) 1702 set_colours.append([]) 1703 x_axis_type_zero.append([]) 1704 symbols.append([]) 1705 symbol_sizes.append([]) 1706 linetype.append([]) 1707 linestyle.append([]) 1708 1709 # Alias the correct spin. 1710 current_spin = spin 1711 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 1712 current_spin = proton 1713 1714 # Loop over the spectrometer frequencies and offsets. 1715 set_index = 0 1716 err = False 1717 colour_index = 0 1718 for frq, offset, mi, oi in loop_frq_offset(exp_type=exp_type, return_indices=True): 1719 # Add a new set for the data at each frequency and offset. 1720 data[graph_index].append([]) 1721 1722 # Add a new label. 1723 if exp_type in EXP_TYPE_LIST_CPMG: 1724 label = "R\\s2eff\\N" 1725 else: 1726 label = "R\\s1\\xr\\B\\N" 1727 if offset != None and frq != None: 1728 label += " (%.1f MHz, %.3f ppm)" % (frq / 1e6, offset) 1729 elif frq != None: 1730 label += " (%.1f MHz)" % (frq / 1e6) 1731 elif offset != None: 1732 label += " (%.3f ppm)" % (offset) 1733 set_labels[ei].append(label) 1734 1735 # The other settings. 1736 set_colours[graph_index].append(colour_order[colour_index]) 1737 x_axis_type_zero[graph_index].append(True) 1738 symbols[graph_index].append(1) 1739 symbol_sizes[graph_index].append(0.45) 1740 linetype[graph_index].append(0) 1741 linestyle[graph_index].append(0) 1742 1743 # Loop over the dispersion points. 1744 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 1745 # The data key. 1746 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 1747 1748 # No data present. 1749 if key not in current_spin.r2eff: 1750 continue 1751 1752 # Add the data. 1753 data[graph_index][set_index].append([point, current_spin.r2eff[key]]) 1754 1755 # Add the error. 1756 if hasattr(current_spin, 'r2eff_err') and key in current_spin.r2eff_err: 1757 err = True 1758 data[graph_index][set_index][-1].append(current_spin.r2eff_err[key]) 1759 1760 # Increment the graph set index. 1761 set_index += 1 1762 colour_index += 1 1763 1764 # Add the back calculated data. 1765 colour_index = 0 1766 for frq, offset, mi, oi in loop_frq_offset(exp_type=exp_type, return_indices=True): 1767 # Add a new set for the data at each frequency and offset. 1768 data[graph_index].append([]) 1769 1770 # Add a new label. 1771 if exp_type in EXP_TYPE_LIST_CPMG: 1772 label = "Back calculated R\\s2eff\\N" 1773 else: 1774 label = "Back calculated R\\s1\\xr\\B\\N" 1775 if offset != None and frq != None: 1776 label += " (%.1f MHz, %.3f ppm)" % (frq / 1e6, offset) 1777 elif frq != None: 1778 label += " (%.1f MHz)" % (frq / 1e6) 1779 elif offset != None: 1780 label += " (%.3f ppm)" % (offset) 1781 set_labels[ei].append(label) 1782 1783 # The other settings. 1784 set_colours[graph_index].append(colour_order[colour_index]) 1785 x_axis_type_zero[graph_index].append(True) 1786 symbols[graph_index].append(4) 1787 symbol_sizes[graph_index].append(0.45) 1788 linetype[graph_index].append(1) 1789 if interpolated_flag: 1790 linestyle[graph_index].append(2) 1791 else: 1792 linestyle[graph_index].append(1) 1793 1794 # Loop over the dispersion points. 1795 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 1796 # The data key. 1797 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 1798 1799 # No data present. 1800 if not hasattr(current_spin, 'r2eff_bc') or key not in current_spin.r2eff_bc: 1801 continue 1802 1803 # Add the data. 1804 data[graph_index][set_index].append([point, current_spin.r2eff_bc[key]]) 1805 1806 # Handle the errors. 1807 if err: 1808 data[graph_index][set_index][-1].append(None) 1809 1810 # Increment the graph set index. 1811 set_index += 1 1812 colour_index += 1 1813 1814 # Add the interpolated back calculated data. 1815 if interpolated_flag: 1816 colour_index = 0 1817 for frq, offset, mi, oi in loop_frq_offset(exp_type=exp_type, return_indices=True): 1818 # Add a new set for the data at each frequency and offset. 1819 data[graph_index].append([]) 1820 1821 # Add a new label. 1822 if exp_type in EXP_TYPE_LIST_CPMG: 1823 label = "R\\s2eff\\N interpolated curve" 1824 else: 1825 label = "R\\s1\\xr\\B\\N interpolated curve" 1826 if offset != None and frq != None: 1827 label += " (%.1f MHz, %.3f ppm)" % (frq / 1e6, offset) 1828 elif frq != None: 1829 label += " (%.1f MHz)" % (frq / 1e6) 1830 elif offset != None: 1831 label += " (%.3f ppm)" % (offset) 1832 set_labels[ei].append(label) 1833 1834 # The other settings. 1835 set_colours[graph_index].append(colour_order[colour_index]) 1836 x_axis_type_zero[graph_index].append(True) 1837 if spin.model in MODEL_LIST_NUMERIC_CPMG: 1838 symbols[graph_index].append(8) 1839 else: 1840 symbols[graph_index].append(0) 1841 symbol_sizes[graph_index].append(0.20) 1842 linetype[graph_index].append(1) 1843 linestyle[graph_index].append(1) 1844 1845 # Loop over the dispersion points. 1846 for di in range(len(back_calc[ei][0][mi][oi])): 1847 # Skip invalid points (values of 1e100). 1848 if back_calc[ei][0][mi][oi][di] > 1e50: 1849 continue 1850 1851 # The X point. 1852 if exp_type in EXP_TYPE_LIST_CPMG: 1853 point = cpmg_frqs_new[ei][mi][oi][di] 1854 else: 1855 point = spin_lock_nu1_new[ei][mi][oi][di] 1856 1857 # Add the data. 1858 data[graph_index][set_index].append([point, back_calc[ei][0][mi][oi][di]]) 1859 1860 # Handle the errors. 1861 if err: 1862 data[graph_index][set_index][-1].append(None) 1863 1864 # Increment the graph set index. 1865 set_index += 1 1866 colour_index += 1 1867 1868 # Add the residuals for statistical comparison. 1869 colour_index = 0 1870 for frq, offset, mi, oi in loop_frq_offset(exp_type=exp_type, return_indices=True): 1871 # Add a new set for the data at each frequency and offset. 1872 data[graph_index].append([]) 1873 1874 # Add a new label. 1875 label = "Residuals" 1876 if offset != None and frq != None: 1877 label += " (%.1f MHz, %.3f ppm)" % (frq / 1e6, offset) 1878 elif frq != None: 1879 label += " (%.1f MHz)" % (frq / 1e6) 1880 elif offset != None: 1881 label += " (%.3f ppm)" % (offset) 1882 set_labels[ei].append(label) 1883 1884 # The other settings. 1885 set_colours[graph_index].append(colour_order[colour_index]) 1886 x_axis_type_zero[graph_index].append(True) 1887 symbols[graph_index].append(9) 1888 symbol_sizes[graph_index].append(0.45) 1889 linetype[graph_index].append(1) 1890 linestyle[graph_index].append(3) 1891 1892 # Loop over the dispersion points. 1893 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 1894 # The data key. 1895 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 1896 1897 # No data present. 1898 if key not in current_spin.r2eff or not hasattr(current_spin, 'r2eff_bc') or key not in current_spin.r2eff_bc: 1899 continue 1900 1901 # Add the data. 1902 data[graph_index][set_index].append([point, current_spin.r2eff[key] - current_spin.r2eff_bc[key]]) 1903 1904 # Handle the errors. 1905 if err: 1906 err = True 1907 data[graph_index][set_index][-1].append(current_spin.r2eff_err[key]) 1908 1909 # Increment the graph set index. 1910 set_index += 1 1911 colour_index += 1 1912 1913 # Increment the graph index. 1914 graph_index += 1 1915 1916 # The axis labels. 1917 if exp_type in EXP_TYPE_LIST_CPMG: 1918 axis_labels.append(['\\qCPMG pulse train frequency (Hz)\\Q', '%s - \\qR\\s2,eff\\N\\Q (rad.s\\S-1\\N)'%exp_type]) 1919 else: 1920 axis_labels.append(['\\qSpin-lock field strength (Hz)\\Q', '\\qR\\s1\\xr\\B\\N\\Q (rad.s\\S-1\\N)']) 1921 1922 # Remove all NaN values. 1923 for i in range(len(data)): 1924 for j in range(len(data[i])): 1925 for k in range(len(data[i][j])): 1926 for l in range(len(data[i][j][k])): 1927 if isNaN(data[i][j][k][l]): 1928 data[i][j][k][l] = 0.0 1929 1930 # Write the header. 1931 title = "Relaxation dispersion plot" 1932 graph_num = len(data) 1933 sets = [] 1934 legend = [] 1935 for gi in range(len(data)): 1936 sets.append(len(data[gi])) 1937 legend.append(False) 1938 legend[0] = True 1939 write_xy_header(file=file, title=title, graph_num=graph_num, sets=sets, set_names=set_labels, set_colours=set_colours, x_axis_type_zero=x_axis_type_zero, symbols=symbols, symbol_sizes=symbol_sizes, linetype=linetype, linestyle=linestyle, axis_labels=axis_labels, legend=legend, legend_box_fill_pattern=[0]*graph_num, legend_char_size=[0.8]*graph_num) 1940 1941 # Write the data. 1942 graph_type = 'xy' 1943 if err: 1944 graph_type = 'xydy' 1945 write_xy_data(data, file=file, graph_type=graph_type) 1946 1947 # Close the file. 1948 file.close() 1949 1950 # Add the file to the results file list. 1951 add_result_file(type='grace', label='Grace', file=file_path) 1952 1953 # Write a python "grace to PNG/EPS/SVG..." conversion script. 1954 # Open the file for writing. 1955 file_name = "grace2images.py" 1956 file = open_write_file(file_name, dir, force) 1957 1958 # Write the file. 1959 script_grace2images(file=file) 1960 1961 # Close the batch script, then make it executable (expanding any ~ characters). 1962 file.close() 1963 if dir: 1964 dir = expanduser(dir) 1965 chmod(dir + sep + file_name, S_IRWXU|S_IRGRP|S_IROTH) 1966 else: 1967 file_name = expanduser(file_name) 1968 chmod(file_name, S_IRWXU|S_IRGRP|S_IROTH)
1969 1970
1971 -def plot_exp_curves(file=None, dir=None, force=None, norm=None):
1972 """Custom 2D Grace plotting function for the exponential curves. 1973 1974 @keyword file: The name of the Grace file to create. 1975 @type file: str 1976 @keyword dir: The optional directory to place the file into. 1977 @type dir: str 1978 @param force: Boolean argument which if True causes the file to be overwritten if it already exists. 1979 @type force: bool 1980 @keyword norm: The normalisation flag which if set to True will cause all graphs to be normalised to a starting value of 1. 1981 @type norm: bool 1982 """ 1983 1984 # Test if the current pipe exists. 1985 pipes.test() 1986 1987 # Test if the sequence data is loaded. 1988 if not exists_mol_res_spin_data(): 1989 raise RelaxNoSequenceError 1990 1991 # Open the file for writing. 1992 file_path = get_file_path(file, dir) 1993 file = open_write_file(file, dir, force) 1994 1995 # Initialise some data structures. 1996 data = [] 1997 set_labels = [] 1998 x_err_flag = False 1999 y_err_flag = False 2000 2001 # 1H MMQ flag. 2002 proton_mmq_flag = has_proton_mmq_cpmg() 2003 2004 # Loop over the spectrometer frequencies. 2005 graph_index = 0 2006 err = False 2007 for exp_type, frq, offset, ei, mi, oi in loop_exp_frq_offset(return_indices=True): 2008 # Loop over the dispersion points. 2009 for point, di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): 2010 # Create a new graph. 2011 data.append([]) 2012 2013 # Loop over each spin. 2014 for spin, id in spin_loop(return_id=True, skip_desel=True): 2015 # Skip protons for MMQ data. 2016 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 2017 continue 2018 2019 # No data present. 2020 if not hasattr(spin, 'intensities'): 2021 continue 2022 2023 # Get the attached proton. 2024 proton = None 2025 if proton_mmq_flag: 2026 proton = return_attached_protons(spin_id)[0] 2027 2028 # Alias the correct spin. 2029 current_spin = spin 2030 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 2031 current_spin = proton 2032 2033 # Append a new set structure and set the name to the spin ID. 2034 data[graph_index].append([]) 2035 if graph_index == 0: 2036 set_labels.append("Spin %s" % id) 2037 2038 # Loop over the relaxation time periods. 2039 for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, point=point): 2040 # The key. 2041 keys = find_intensity_keys(exp_type=exp_type, frq=frq, offset=offset, point=point, time=time) 2042 2043 # Loop over each key. 2044 for key in keys: 2045 # No key present. 2046 if key not in current_spin.intensities: 2047 continue 2048 2049 # Add the data. 2050 if hasattr(current_spin, 'intensity_err'): 2051 data[graph_index][-1].append([time, current_spin.intensities[key], spin.intensity_err[key]]) 2052 err = True 2053 else: 2054 data[graph_index][-1].append([time, current_spin.intensities[key]]) 2055 2056 # Increment the frq index. 2057 graph_index += 1 2058 2059 # The axis labels. 2060 axis_labels = ['Relaxation time period (s)', 'Peak intensities'] 2061 2062 # Write the header. 2063 graph_num = len(data) 2064 sets = [] 2065 for gi in range(graph_num): 2066 sets.append(len(data[gi])) 2067 write_xy_header(file=file, graph_num=graph_num, sets=sets, set_names=[set_labels]*graph_num, axis_labels=[axis_labels]*graph_num, norm=[norm]*graph_num) 2068 2069 # Write the data. 2070 graph_type = 'xy' 2071 if err: 2072 graph_type = 'xydy' 2073 write_xy_data(data, file=file, graph_type=graph_type, norm=[norm]*graph_num) 2074 2075 # Close the file. 2076 file.close() 2077 2078 # Add the file to the results file list. 2079 add_result_file(type='grace', label='Grace', file=file_path)
2080 2081
2082 -def r2eff_read(id=None, file=None, dir=None, disp_frq=None, offset=None, spin_id_col=None, mol_name_col=None, res_num_col=None, res_name_col=None, spin_num_col=None, spin_name_col=None, data_col=None, error_col=None, sep=None):
2083 """Read R2eff/R1rho values directly from a file whereby each row corresponds to a different spin. 2084 2085 @keyword id: The experiment ID string to associate the data with. 2086 @type id: str 2087 @keyword file: The name of the file to open. 2088 @type file: str 2089 @keyword dir: The directory containing the file (defaults to the current directory if None). 2090 @type dir: str or None 2091 @keyword disp_frq: For CPMG-type data, the frequency of the CPMG pulse train. For R1rho-type data, the spin-lock field strength (nu1). The units must be Hertz. 2092 @type disp_frq: float 2093 @keyword offset: For R1rho-type data, the spin-lock offset value in ppm. 2094 @type offset: None or float 2095 @keyword spin_id_col: The column containing the spin ID strings. If supplied, the mol_name_col, res_name_col, res_num_col, spin_name_col, and spin_num_col arguments must be none. 2096 @type spin_id_col: int or None 2097 @keyword mol_name_col: The column containing the molecule name information. If supplied, spin_id_col must be None. 2098 @type mol_name_col: int or None 2099 @keyword res_name_col: The column containing the residue name information. If supplied, spin_id_col must be None. 2100 @type res_name_col: int or None 2101 @keyword res_num_col: The column containing the residue number information. If supplied, spin_id_col must be None. 2102 @type res_num_col: int or None 2103 @keyword spin_name_col: The column containing the spin name information. If supplied, spin_id_col must be None. 2104 @type spin_name_col: int or None 2105 @keyword spin_num_col: The column containing the spin number information. If supplied, spin_id_col must be None. 2106 @type spin_num_col: int or None 2107 @keyword data_col: The column containing the R2eff/R1rho data in Hz. 2108 @type data_col: int or None 2109 @keyword error_col: The column containing the R2eff/R1rho errors. 2110 @type error_col: int or None 2111 @keyword sep: The column separator which, if None, defaults to whitespace. 2112 @type sep: str or None 2113 """ 2114 2115 # Data checks. 2116 pipes.test() 2117 check_mol_res_spin_data() 2118 check_frequency(id=id) 2119 check_exp_type(id=id) 2120 2121 # Store the spectrum ID. 2122 add_spectrum_id(id) 2123 2124 # Get the metadata. 2125 frq = get_frequency(id=id) 2126 exp_type = get_exp_type(id=id) 2127 2128 # Loop over the data. 2129 data_flag = False 2130 mol_names = [] 2131 res_nums = [] 2132 res_names = [] 2133 spin_nums = [] 2134 spin_names = [] 2135 values = [] 2136 errors = [] 2137 for data in read_spin_data(file=file, dir=dir, spin_id_col=spin_id_col, mol_name_col=mol_name_col, res_num_col=res_num_col, res_name_col=res_name_col, spin_num_col=spin_num_col, spin_name_col=spin_name_col, data_col=data_col, error_col=error_col, sep=sep): 2138 # Unpack. 2139 if data_col and error_col: 2140 mol_name, res_num, res_name, spin_num, spin_name, value, error = data 2141 elif data_col: 2142 mol_name, res_num, res_name, spin_num, spin_name, value = data 2143 error = None 2144 else: 2145 mol_name, res_num, res_name, spin_num, spin_name, error = data 2146 value = None 2147 2148 # Test the error value (cannot be 0.0). 2149 if error == 0.0: 2150 raise RelaxError("An invalid error value of zero has been encountered.") 2151 2152 # Get the corresponding spin container. 2153 spin_id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=spin_num, spin_name=spin_name) 2154 spin = return_spin(spin_id) 2155 if spin == None: 2156 warn(RelaxNoSpinWarning(spin_id)) 2157 continue 2158 2159 # The dispersion point key. 2160 point_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=disp_frq) 2161 2162 # Store the R2eff data. 2163 if data_col: 2164 # Initialise if necessary. 2165 if not hasattr(spin, 'r2eff'): 2166 spin.r2eff = {} 2167 2168 # Store. 2169 spin.r2eff[point_key] = value 2170 2171 # Store the R2eff error. 2172 if error_col: 2173 # Initialise if necessary. 2174 if not hasattr(spin, 'r2eff_err'): 2175 spin.r2eff_err = {} 2176 2177 # Store. 2178 spin.r2eff_err[point_key] = error 2179 2180 # Data added. 2181 data_flag = True 2182 2183 # Append the data for printout. 2184 mol_names.append(mol_name) 2185 res_nums.append(res_num) 2186 res_names.append(res_name) 2187 spin_nums.append(spin_num) 2188 spin_names.append(spin_name) 2189 values.append(value) 2190 errors.append(error) 2191 2192 # Print out. 2193 write_spin_data(file=sys.stdout, mol_names=mol_names, res_nums=res_nums, res_names=res_names, spin_nums=spin_nums, spin_names=spin_names, data=values, data_name='R2eff', error=errors, error_name='R2eff_error') 2194 2195 # Update the global structures. 2196 if data_flag: 2197 # Set the dispersion point frequency. 2198 if exp_type in EXP_TYPE_LIST_CPMG: 2199 cpmg_frq(spectrum_id=id, cpmg_frq=disp_frq) 2200 else: 2201 spin_lock_field(spectrum_id=id, field=disp_frq)
2202 2203
2204 -def r2eff_read_spin(id=None, spin_id=None, file=None, dir=None, disp_point_col=None, offset_col=None, data_col=None, error_col=None, sep=None):
2205 """Read R2eff/R1rho values from file whereby each row is a different dispersion point. 2206 2207 @keyword id: The experiment ID string to associate the data with. This will be modified to include the dispersion point data as "%s_%s" % (id, disp_point). 2208 @type id: str 2209 @keyword spin_id: The spin ID string. 2210 @type spin_id: str 2211 @keyword file: The name of the file to open. 2212 @type file: str 2213 @keyword dir: The directory containing the file (defaults to the current directory if None). 2214 @type dir: str or None 2215 @keyword disp_point_col: The column containing the dispersion point information. For CPMG-type data, this is the frequency of the CPMG pulse train. For R1rho-type data, this is the spin-lock field strength (nu1). The units must be Hertz. 2216 @type disp_point_col: int 2217 @keyword offset_col: This is for R1rho data - the dispersion point column can be substituted for the offset values in Hertz. 2218 @type offset_col: None or int 2219 @keyword data_col: The column containing the R2eff/R1rho data in Hz. 2220 @type data_col: int 2221 @keyword error_col: The column containing the R2eff/R1rho errors. 2222 @type error_col: int 2223 @keyword sep: The column separator which, if None, defaults to whitespace. 2224 @type sep: str or None 2225 """ 2226 2227 # Data checks. 2228 pipes.test() 2229 check_mol_res_spin_data() 2230 2231 # Get the spin. 2232 spin = return_spin(spin_id) 2233 if spin == None: 2234 raise RelaxNoSpinError(spin_id) 2235 2236 # Extract the data from the file, removing comments and blank lines. 2237 file_data = extract_data(file, dir, sep=sep) 2238 file_data = strip(file_data) 2239 2240 # Loop over the data. 2241 data = [] 2242 new_ids = [] 2243 for line in file_data: 2244 # Invalid columns. 2245 if disp_point_col != None and disp_point_col > len(line): 2246 warn(RelaxWarning("The data %s is invalid, no dispersion point column can be found." % line)) 2247 continue 2248 if offset_col != None and offset_col > len(line): 2249 warn(RelaxWarning("The data %s is invalid, no offset column can be found." % line)) 2250 continue 2251 if data_col > len(line): 2252 warn(RelaxWarning("The R2eff/R1rho data %s is invalid, no data column can be found." % line)) 2253 continue 2254 if error_col > len(line): 2255 warn(RelaxWarning("The R2eff/R1rho data %s is invalid, no error column can be found." % line)) 2256 continue 2257 2258 # Unpack. 2259 if disp_point_col != None: 2260 ref_data = line[disp_point_col-1] 2261 elif offset_col != None: 2262 ref_data = line[offset_col-1] 2263 value = line[data_col-1] 2264 error = line[error_col-1] 2265 2266 # Convert and check the dispersion point or offset. 2267 try: 2268 ref_data = float(ref_data) 2269 except ValueError: 2270 if disp_point_col != None: 2271 warn(RelaxWarning("The dispersion point data of the line %s is invalid." % line)) 2272 elif offset_col != None: 2273 warn(RelaxWarning("The offset data of the line %s is invalid." % line)) 2274 continue 2275 2276 # Convert and check the value. 2277 if value == 'None': 2278 value = None 2279 if value != None: 2280 try: 2281 value = float(value) 2282 except ValueError: 2283 warn(RelaxWarning("The R2eff/R1rho value of the line %s is invalid." % line)) 2284 continue 2285 2286 # Convert and check the error. 2287 if error == 'None': 2288 error = None 2289 if error != None: 2290 try: 2291 error = float(error) 2292 except ValueError: 2293 warn(RelaxWarning("The R2eff/R1rho error of the line %s is invalid." % line)) 2294 continue 2295 2296 # Test the error value (cannot be 0.0). 2297 if error == 0.0: 2298 raise RelaxError("An invalid error value of zero has been encountered.") 2299 2300 # Find the matching spectrum ID. 2301 new_id = None 2302 for spectrum_id in cdp.spectrum_ids: 2303 # Skip IDs which don't start with the base ID. 2304 if not search("^%s"%id, spectrum_id): 2305 continue 2306 2307 # Find a close enough dispersion point (to one decimal place to allow for user truncation). 2308 if disp_point_col != None: 2309 if hasattr(cdp, 'cpmg_frqs') and spectrum_id in cdp.cpmg_frqs: 2310 if abs(ref_data - cdp.cpmg_frqs[spectrum_id]) < 0.1: 2311 new_id = spectrum_id 2312 break 2313 if hasattr(cdp, 'spin_lock_nu1') and spectrum_id in cdp.spin_lock_nu1: 2314 if abs(ref_data - cdp.spin_lock_nu1[spectrum_id]) < 0.1: 2315 new_id = spectrum_id 2316 break 2317 2318 # Find a close enough offset (to one decimal place to allow for user truncation). 2319 elif offset_col != None: 2320 if hasattr(cdp, 'spin_lock_offset') and spectrum_id in cdp.spin_lock_offset: 2321 # The sign to multiply offsets by. 2322 sign = 1.0 2323 if spin.isotope == '15N': 2324 sign = -1.0 2325 2326 # Convert the data. 2327 data_new = sign * frequency_to_ppm(frq=ref_data, B0=cdp.spectrometer_frq[spectrum_id], isotope=spin.isotope) 2328 2329 # Store the ID. 2330 if abs(data_new - cdp.spin_lock_offset[spectrum_id]) < 0.1: 2331 new_id = spectrum_id 2332 break 2333 2334 # No match. 2335 if new_id == None: 2336 if disp_point_col != None: 2337 raise RelaxError("The experiment ID corresponding to the base ID '%s' and the dispersion point '%s' could not be found." % (id, ref_data)) 2338 if offset_col != None: 2339 raise RelaxError("The experiment ID corresponding to the base ID '%s' and the offset '%s' could not be found." % (id, data_new)) 2340 2341 # Add the ID to the list. 2342 new_ids.append(new_id) 2343 2344 # Data checks. 2345 check_frequency(id=new_id) 2346 check_exp_type(id=new_id) 2347 2348 # Store the spectrum ID. 2349 add_spectrum_id(new_id) 2350 2351 # Get the metadata. 2352 frq = get_frequency(id=new_id) 2353 exp_type = get_exp_type(id=new_id) 2354 2355 # The dispersion point key. 2356 if disp_point_col != None: 2357 disp_point = ref_data 2358 offset = 0.0 2359 elif offset_col != None: 2360 disp_point = cdp.spin_lock_nu1[new_id] 2361 offset = data_new 2362 point_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=disp_point) 2363 2364 # Store the R2eff data. 2365 if data_col: 2366 # Initialise if necessary. 2367 if not hasattr(spin, 'r2eff'): 2368 spin.r2eff = {} 2369 2370 # Store. 2371 spin.r2eff[point_key] = value 2372 2373 # Store the R2eff error. 2374 if error_col: 2375 # Initialise if necessary. 2376 if not hasattr(spin, 'r2eff_err'): 2377 spin.r2eff_err = {} 2378 2379 # Store. 2380 spin.r2eff_err[point_key] = error 2381 2382 # Append the data for printout. 2383 if disp_point_col != None: 2384 data.append(["%-40s" % point_key, "%20.15f" % disp_point, "%20.15f" % value, "%20.15f" % error]) 2385 else: 2386 data.append(["%-40s" % point_key, "%20.15f" % offset, "%20.15f" % value, "%20.15f" % error]) 2387 2388 # Data added. 2389 data_flag = True 2390 2391 # No data, so fail hard! 2392 if not len(data): 2393 raise RelaxError("No R2eff/R1rho data could be extracted.") 2394 2395 # Print out. 2396 print("The following R2eff/R1rho data has been loaded into the relax data store:\n") 2397 if disp_point_col != None: 2398 write_data(out=sys.stdout, headings=["R2eff_key", "Disp_point", "R2eff", "R2eff_error"], data=data) 2399 else: 2400 write_data(out=sys.stdout, headings=["R2eff_key", "Offset (ppm)", "R2eff", "R2eff_error"], data=data)
2401 2402
2403 -def randomise_R1(spin=None, ri_id=None, N=None):
2404 """Randomise the R1 data for the given spin for use in the Monte Carlo simulations. 2405 2406 @keyword spin: The spin container to randomise the data for. 2407 @type spin: SpinContainer instance 2408 @keyword ri_id: The relaxation data ID string. 2409 @type ri_id: str 2410 @keyword N: The number of randomisations to perform. 2411 @type N: int 2412 """ 2413 2414 # The data already exists. 2415 if hasattr(spin, 'ri_data_sim') and ri_id in spin.ri_data_sim: 2416 return 2417 2418 # Initialise the structure. 2419 if not hasattr(spin, 'ri_data_sim'): 2420 spin.ri_data_sim = {} 2421 spin.ri_data_sim[ri_id] = [] 2422 2423 # Randomise. 2424 for i in range(N): 2425 spin.ri_data_sim[ri_id].append(gauss(spin.ri_data[ri_id], spin.ri_data_err[ri_id]))
2426 2427
2428 -def relax_time(time=0.0, spectrum_id=None):
2429 """Set the relaxation time period associated with a given spectrum. 2430 2431 @keyword time: The time, in seconds, of the relaxation period. 2432 @type time: float 2433 @keyword spectrum_id: The spectrum identification string. 2434 @type spectrum_id: str 2435 """ 2436 2437 # Test if the spectrum id exists. 2438 if spectrum_id not in cdp.spectrum_ids: 2439 raise RelaxNoSpectraError(spectrum_id) 2440 2441 # Initialise the global relaxation time data structures if needed. 2442 if not hasattr(cdp, 'relax_times'): 2443 cdp.relax_times = {} 2444 if not hasattr(cdp, 'relax_time_list'): 2445 cdp.relax_time_list = [] 2446 2447 # Add the time, converting to a float if needed. 2448 cdp.relax_times[spectrum_id] = float(time) 2449 2450 # The unique time points. 2451 if cdp.relax_times[spectrum_id] not in cdp.relax_time_list: 2452 cdp.relax_time_list.append(cdp.relax_times[spectrum_id]) 2453 cdp.relax_time_list.sort() 2454 2455 # Update the exponential time point count. 2456 cdp.num_time_pts = len(cdp.relax_time_list) 2457 2458 # Printout. 2459 print("Setting the '%s' spectrum relaxation time period to %s s." % (spectrum_id, cdp.relax_times[spectrum_id]))
2460 2461
2462 -def return_cpmg_frqs(ref_flag=True):
2463 """Return the list of nu_CPMG frequencies. 2464 2465 @keyword ref_flag: A flag which if False will cause the reference spectrum frequency of None to be removed from the list. 2466 @type ref_flag: bool 2467 @return: The list of nu_CPMG frequencies in Hz. It has the dimensions {Ei, Mi, Oi}. 2468 @rtype: rank-2 list of numpy rank-1 float64 arrays 2469 """ 2470 2471 # No data. 2472 if not hasattr(cdp, 'cpmg_frqs_list'): 2473 return None 2474 2475 # Initialise. 2476 cpmg_frqs = [] 2477 2478 # First loop over the experiment types. 2479 for exp_type, ei in loop_exp(return_indices=True): 2480 # Add a new dimension. 2481 cpmg_frqs.append([]) 2482 2483 # Then loop over the spectrometer frequencies. 2484 for frq, mi in loop_frq(return_indices=True): 2485 # Add a new dimension. 2486 cpmg_frqs[ei].append([]) 2487 2488 # Loop over the offsets. 2489 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 2490 # Add a new dimension. 2491 cpmg_frqs[ei][mi].append([]) 2492 2493 # Loop over the fields. 2494 for point in cdp.cpmg_frqs_list: 2495 # Skip reference points. 2496 if (not ref_flag) and point == None: 2497 continue 2498 2499 # Find a matching experiment ID. 2500 found = False 2501 for id in cdp.exp_type.keys(): 2502 # Skip non-matching experiments. 2503 if cdp.exp_type[id] != exp_type: 2504 continue 2505 2506 # Skip non-matching spectrometer frequencies. 2507 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 2508 continue 2509 2510 # Skip non-matching points. 2511 if cdp.cpmg_frqs[id] != point: 2512 continue 2513 2514 # Found. 2515 found = True 2516 break 2517 2518 # No data. 2519 if not found: 2520 continue 2521 2522 # Add the data. 2523 cpmg_frqs[ei][mi][oi].append(point) 2524 2525 # Convert to a numpy array. 2526 cpmg_frqs[ei][mi][oi] = array(cpmg_frqs[ei][mi][oi], float64) 2527 2528 # Return the data. 2529 return cpmg_frqs
2530 2531
2532 -def return_cpmg_frqs_single(exp_type=None, frq=None, offset=None, time=None, ref_flag=True):
2533 """Return the list of nu_CPMG frequencies. 2534 2535 @keyword exp_type: The experiment type. 2536 @type exp_type: str 2537 @keyword frq: The spectrometer frequency in Hz. 2538 @type frq: float 2539 @keyword offset: The hard pulse offset, if desired. 2540 @type offset: None or float 2541 @keyword time: The relaxation time period. 2542 @type time: float 2543 @keyword ref_flag: A flag which if False will cause the reference spectrum frequency of None to be removed from the list. 2544 @type ref_flag: bool 2545 @return: The list of nu_CPMG frequencies in Hz. 2546 @rtype: numpy rank-1 float64 array 2547 """ 2548 2549 # No data. 2550 if not hasattr(cdp, 'cpmg_frqs_list'): 2551 return None 2552 2553 # Initialise. 2554 cpmg_frqs = [] 2555 2556 # Loop over the points. 2557 for point in cdp.cpmg_frqs_list: 2558 # Skip reference points. 2559 if (not ref_flag) and point == None: 2560 continue 2561 2562 # Find a matching experiment ID. 2563 found = False 2564 for id in cdp.exp_type.keys(): 2565 # Skip non-matching experiments. 2566 if cdp.exp_type[id] != exp_type: 2567 continue 2568 2569 # Skip non-matching spectrometer frequencies. 2570 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 2571 continue 2572 2573 # Skip non-matching offsets. 2574 if offset != None and hasattr(cdp, 'spin_lock_offset') and cdp.spin_lock_offset[id] != offset: 2575 continue 2576 2577 # Skip non-matching time points. 2578 if time != None and hasattr(cdp, 'relax_times') and cdp.relax_times[id] != time: 2579 continue 2580 2581 # Skip non-matching points. 2582 if cdp.cpmg_frqs[id] != point: 2583 continue 2584 2585 # Found. 2586 found = True 2587 break 2588 2589 # No data. 2590 if not found: 2591 continue 2592 2593 # Add the data. 2594 cpmg_frqs.append(point) 2595 2596 # Return the data as a numpy array. 2597 return array(cpmg_frqs, float64)
2598 2599
2600 -def return_index_from_disp_point(value, exp_type=None):
2601 """Convert the dispersion point data into the corresponding index. 2602 2603 @param value: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 2604 @type value: float 2605 @keyword exp_type: The experiment type. 2606 @type exp_type: str 2607 @return: The corresponding index. 2608 @rtype: int 2609 """ 2610 2611 # Check. 2612 if exp_type == None: 2613 raise RelaxError("The experiment type has not been supplied.") 2614 2615 # Initialise. 2616 index = 0 2617 ref_correction = False 2618 2619 # CPMG-type experiments. 2620 if exp_type in EXP_TYPE_LIST_CPMG: 2621 index = cdp.cpmg_frqs_list.index(value) 2622 if None in cdp.cpmg_frqs_list: 2623 ref_correction = True 2624 2625 # R1rho-type experiments. 2626 elif exp_type in EXP_TYPE_LIST_R1RHO: 2627 index = cdp.spin_lock_nu1_list.index(value) 2628 if None in cdp.spin_lock_nu1_list: 2629 ref_correction = True 2630 2631 # Remove the reference point (always at index 0). 2632 for id in loop_spectrum_ids(exp_type=exp_type): 2633 if ref_correction and get_curve_type(id) == 'fixed time': 2634 index -= 1 2635 break 2636 2637 # Return the index. 2638 return index
2639 2640
2641 -def return_index_from_exp_type(exp_type=None):
2642 """Convert the experiment type into the corresponding index. 2643 2644 @keyword exp_type: The experiment type. 2645 @type exp_type: str 2646 @return: The corresponding index. 2647 @rtype: int 2648 """ 2649 2650 # Check. 2651 if exp_type == None: 2652 raise RelaxError("The experiment type has not been supplied.") 2653 2654 # Return the index. 2655 if exp_type in cdp.exp_type_list: 2656 return cdp.exp_type_list.index(exp_type) 2657 2658 # The number of experiments. 2659 num = len(cdp.exp_type_list)
2660 2661
2662 -def return_index_from_frq(value):
2663 """Convert the dispersion point data into the corresponding index. 2664 2665 @param value: The spectrometer frequency in Hz. 2666 @type value: float 2667 @return: The corresponding index. 2668 @rtype: int 2669 """ 2670 2671 # No frequency present. 2672 if value == None: 2673 return 0 2674 2675 # Return the index. 2676 return cdp.spectrometer_frq_list.index(value)
2677 2678
2679 -def return_index_from_disp_point_key(key, exp_type=None):
2680 """Convert the dispersion point key into the corresponding index. 2681 2682 @keyword exp_type: The experiment type. 2683 @type exp_type: str 2684 @param key: The dispersion point or R2eff/R1rho key. 2685 @type key: str 2686 @return: The corresponding index. 2687 @rtype: int 2688 """ 2689 2690 # Check. 2691 if exp_type == None: 2692 raise RelaxError("The experiment type has not been supplied.") 2693 2694 # CPMG-type experiments. 2695 if exp_type in EXP_TYPE_LIST_CPMG: 2696 return return_index_from_disp_point(cdp.cpmg_frqs[key], exp_type=exp_type) 2697 2698 # R1rho-type experiments. 2699 elif exp_type in EXP_TYPE_LIST_R1RHO: 2700 return return_index_from_disp_point(cdp.spin_lock_nu1[key], exp_type=exp_type)
2701 2702
2703 -def return_key_from_di(mi=None, di=None):
2704 """Convert the dispersion point index into the corresponding key. 2705 2706 @keyword mi: The spectrometer frequency index. 2707 @type mi: int 2708 @keyword di: The dispersion point or R2eff/R1rho index. 2709 @type di: int 2710 @return: The corresponding key. 2711 @rtype: str 2712 """ 2713 2714 # Insert the reference point (always at index 0). 2715 if has_fixed_time_exp_type(): 2716 di += 1 2717 2718 # The frequency. 2719 frq = return_value_from_frq_index(mi) 2720 2721 # CPMG data. 2722 if exp_type in EXP_TYPE_LIST_CPMG: 2723 point = cdp.cpmg_frqs_list[di] 2724 points = cdp.cpmg_frqs 2725 2726 # R1rho data. 2727 else: 2728 point = cdp.spin_lock_nu1_list[di] 2729 points = cdp.spin_lock_nu1 2730 2731 # Find the keys matching the dispersion point. 2732 key_list = [] 2733 all_keys = points.keys() 2734 for key in all_keys: 2735 if points[key] == point: 2736 key_list.append(key) 2737 2738 # Return the key. 2739 return key
2740 2741
2742 -def return_offset_data(spins=None, spin_ids=None, field_count=None, fields=None):
2743 """Return numpy arrays of the chemical shifts, offsets and tilt angles. 2744 2745 Indices 2746 ======= 2747 2748 The data structures consist of many different index types. These are: 2749 2750 - Ei: The index for each experiment type. 2751 - Si: The index for each spin of the spin cluster. 2752 - Mi: The index for each magnetic field strength. 2753 - Oi: The index for each spin-lock offset. In the case of CPMG-type data, this index is always zero. 2754 - Di: The index for each dispersion point (either the spin-lock field strength or the nu_CPMG frequency). 2755 2756 2757 @keyword spins: The list of spin containers in the cluster. 2758 @type spins: list of SpinContainer instances 2759 @keyword spin_ids: The list of spin IDs for the cluster. 2760 @type spin_ids: list of str 2761 @keyword field_count: The number of spectrometer field strengths. This may not be equal to the length of the fields list as the user may not have set the field strength. 2762 @type field_count: int 2763 @keyword fields: The spin-lock field strengths to use instead of the user loaded values - to enable interpolation. The dimensions are {Ei, Mi}. 2764 @type fields: rank-2 list of floats 2765 @return: The numpy array structures of the chemical shifts in rad/s {Ei, Si, Mi}, spin-lock offsets in rad/s {Ei, Si, Mi, Oi}, rotating frame tilt angles {Ei, Si, Mi, Oi, Di}, the average resonance offset in the rotating frame in rad/s {Ei, Si, Mi, Oi, Di} and the effective field in rotating frame in rad/s {Ei, Si, Mi, Oi, Di}. 2766 @rtype: rank-3 list of floats, rank-4 list of floats, rank-5 list of floats 2767 """ 2768 2769 # The counts. 2770 exp_num = num_exp_types() 2771 spin_num = 0 2772 for spin in spins: 2773 if spin.select: 2774 spin_num += 1 2775 2776 # Initialise the data structures for the target function. 2777 fields_orig = fields 2778 shifts = [] 2779 offsets = [] 2780 theta = [] 2781 Domega = [] 2782 w_e = [] 2783 for exp_type, ei in loop_exp(return_indices=True): 2784 shifts.append([]) 2785 offsets.append([]) 2786 theta.append([]) 2787 Domega.append([]) 2788 w_e.append([]) 2789 for si in range(spin_num): 2790 shifts[ei].append([]) 2791 offsets[ei].append([]) 2792 theta[ei].append([]) 2793 Domega[ei].append([]) 2794 w_e[ei].append([]) 2795 for frq, mi in loop_frq(return_indices=True): 2796 shifts[ei][si].append(None) 2797 offsets[ei][si].append([]) 2798 theta[ei][si].append([]) 2799 Domega[ei][si].append([]) 2800 w_e[ei][si].append([]) 2801 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 2802 offsets[ei][si][mi].append(None) 2803 theta[ei][si][mi].append([]) 2804 Domega[ei][si][mi].append([]) 2805 w_e[ei][si][mi].append([]) 2806 2807 # Assemble the data. 2808 data_flag = False 2809 si = 0 2810 for spin_index in range(len(spins)): 2811 # Skip deselected spins. 2812 if not spins[spin_index].select: 2813 continue 2814 2815 # Alias the spin. 2816 spin = spins[spin_index] 2817 spin_id = spin_ids[spin_index] 2818 2819 # No data. 2820 shift = 0.0 2821 if hasattr(spin, 'chemical_shift'): 2822 shift = spin.chemical_shift 2823 elif has_r1rho_exp_type(): 2824 warn(RelaxWarning("The chemical shift for the spin '%s' cannot be found. Be careful, it is being set to 0.0 ppm so offset calculations will probably be wrong!" % spin_id)) 2825 2826 # Loop over the experiments and spectrometer frequencies. 2827 data_flag = True 2828 for exp_type, frq, offset, ei, mi, oi in loop_exp_frq_offset(return_indices=True): 2829 # The R1rho and off-resonance R1rho flag. 2830 r1rho_flag = False 2831 if exp_type in EXP_TYPE_LIST_R1RHO: 2832 r1rho_flag = True 2833 r1rho_off_flag = False 2834 if exp_type in [MODEL_DPL94, MODEL_TP02, MODEL_TAP03, MODEL_MP05, MODEL_NS_R1RHO_2SITE]: 2835 r1rho_off_flag = True 2836 2837 # Make sure offset data exists for off-resonance R1rho-type experiments. 2838 if r1rho_off_flag and not hasattr(cdp, 'spin_lock_offset'): 2839 raise RelaxError("The spin-lock offsets have not been set.") 2840 2841 # The spin-lock data. 2842 if fields_orig != None: 2843 fields = fields_orig[ei][mi][oi] 2844 else: 2845 if not r1rho_flag: 2846 fields = return_cpmg_frqs_single(exp_type=exp_type, frq=frq, offset=offset, ref_flag=False) 2847 else: 2848 fields = return_spin_lock_nu1_single(exp_type=exp_type, frq=frq, offset=offset, ref_flag=False) 2849 2850 # Convert the shift from ppm to rad/s and store it. 2851 shifts[ei][si][mi] = frequency_to_rad_per_s(frq=shift, B0=frq, isotope=spin.isotope) 2852 2853 # Find a matching experiment ID. 2854 found = False 2855 for id in cdp.exp_type.keys(): 2856 # Skip non-matching experiments. 2857 if cdp.exp_type[id] != exp_type: 2858 continue 2859 2860 # Skip non-matching spectrometer frequencies. 2861 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 2862 continue 2863 2864 # Skip non-matching offsets. 2865 if r1rho_flag and hasattr(cdp, 'spin_lock_offset') and cdp.spin_lock_offset[id] != offset: 2866 continue 2867 2868 # Found. 2869 found = True 2870 break 2871 2872 # No data. 2873 if not found: 2874 continue 2875 2876 # Store the offset in rad/s. Only once and using the first key. 2877 if offsets[ei][si][mi][oi] == None: 2878 if r1rho_flag and hasattr(cdp, 'spin_lock_offset'): 2879 offsets[ei][si][mi][oi] = frequency_to_rad_per_s(frq=cdp.spin_lock_offset[id], B0=frq, isotope=spin.isotope) 2880 else: 2881 offsets[ei][si][mi][oi] = 0.0 2882 2883 # Loop over the dispersion points. 2884 for di in range(len(fields)): 2885 # Alias the point. 2886 point = fields[di] 2887 2888 # Skip reference spectra. 2889 if point == None: 2890 continue 2891 2892 # Calculate the tilt angle. 2893 omega1 = point * 2.0 * pi 2894 Delta_omega = shifts[ei][si][mi] - offsets[ei][si][mi][oi] 2895 Domega[ei][si][mi][oi].append(Delta_omega) 2896 if Delta_omega == 0.0: 2897 theta[ei][si][mi][oi].append(pi / 2.0) 2898 # Calculate the theta angle describing the tilted rotating frame relative to the laboratory. 2899 # If Delta_omega is negative, there follow the symmetry of atan, that atan(-x) = - atan(x). 2900 # Then it should be: theta = pi + atan(-x) = pi - atan(x) = pi - abs(atan( +/- x)) 2901 elif omega1 / Delta_omega > 0 : 2902 theta[ei][si][mi][oi].append(atan(omega1 / Delta_omega)) 2903 else: 2904 theta[ei][si][mi][oi].append(pi + atan(omega1 / Delta_omega)) 2905 2906 # Calculate effective field in rotating frame 2907 w_eff = sqrt( Delta_omega*Delta_omega + omega1*omega1 ) 2908 w_e[ei][si][mi][oi].append(w_eff) 2909 2910 # Increment the spin index. 2911 si += 1 2912 2913 # No shift data for the spin cluster. 2914 if not data_flag: 2915 return None, None, None 2916 2917 # Convert to numpy arrays. 2918 #for ei in range(exp_num): 2919 # for si in range(spin_num): 2920 # for mi in range(field_count): 2921 # theta[ei][si][mi] = array(theta[ei][si][mi], float64) 2922 2923 # Return the structures. 2924 return shifts, offsets, theta, Domega, w_e
2925 2926
2927 -def return_param_key_from_data(exp_type=None, frq=0.0, offset=0.0, point=0.0):
2928 """Generate the unique key from the spectrometer frequency and dispersion point. 2929 2930 @keyword exp_type: The experiment type. 2931 @type exp_type: str 2932 @keyword frq: The spectrometer frequency in Hz. 2933 @type frq: float 2934 @keyword offset: The optional offset value for off-resonance R1rho-type data. 2935 @type offset: None or float 2936 @keyword point: The dispersion point data (either the spin-lock field strength in Hz or the nu_CPMG frequency in Hz). 2937 @type point: float 2938 @return: The unique key. 2939 @rtype: str 2940 """ 2941 2942 # Convert the experiment type. 2943 if exp_type == None: 2944 raise RelaxError("The experiment type must be supplied.") 2945 exp_type = exp_type.replace(' ', '_').lower() 2946 2947 # Convert None values. 2948 if frq == None: 2949 frq = 0.0 2950 if offset == None: 2951 offset = 0.0 2952 if point == None: 2953 point = 0.0 2954 2955 # Generate the unique key. 2956 key = "%s_%.8f_%.3f_%.3f" % (exp_type, frq/1e6, offset, point) 2957 2958 # Return the unique key. 2959 return key
2960 2961
2962 -def return_r1_data(spins=None, spin_ids=None, field_count=None, sim_index=None):
2963 """Return the R1 data structures for off-resonance R1rho experiments. 2964 2965 @keyword spins: The list of spin containers in the cluster. 2966 @type spins: list of SpinContainer instances 2967 @keyword spin_ids: The list of spin IDs for the cluster. 2968 @type spin_ids: list of str 2969 @keyword field_count: The number of spectrometer field strengths. This may not be equal to the length of the fields list as the user may not have set the field strength. 2970 @type field_count: int 2971 @keyword sim_index: The index of the simulation to return the R1 data of. This should be None if the normal data is required. 2972 @type sim_index: None or int 2973 @return: The R1 relaxation data. 2974 @rtype: numpy rank-2 float array 2975 """ 2976 2977 # The spin count. 2978 spin_num = count_spins(spins) 2979 2980 # Initialise the data structure. 2981 r1 = -ones((spin_num, field_count), float64) 2982 2983 # Check for the presence of data. 2984 if not hasattr(cdp, 'ri_ids'): 2985 if has_r1rho_exp_type(): 2986 warn(RelaxWarning("No R1 relaxation data has been loaded. This is essential for the proper handling of offsets in off-resonance R1rho experiments.")) 2987 return 0.0 * r1 2988 2989 # Loop over the Rx IDs. 2990 flags = [False]*field_count 2991 for ri_id in cdp.ri_ids: 2992 # Only use R1 data. 2993 if cdp.ri_type[ri_id] != 'R1': 2994 continue 2995 2996 # The frequency. 2997 frq = cdp.spectrometer_frq[ri_id] 2998 mi = return_index_from_frq(frq) 2999 3000 # Flip the flag. 3001 flags[mi] = True 3002 3003 # Spin loop. 3004 for si in range(spin_num): 3005 # FIXME: This is a kludge - the data randomisation needs to be incorporated into the dispersion base_data_loop() method and the standard Monte Carlo simulation pathway used. 3006 # Randomise the R1 data, when required. 3007 if sim_index != None and (not hasattr(spins[si], 'ri_data_sim') or ri_id not in spins[si].ri_data_sim): 3008 randomise_R1(spin=spins[si], ri_id=ri_id, N=cdp.sim_number) 3009 3010 # Store the data. 3011 if sim_index != None: 3012 r1[si, mi] = spins[si].ri_data_sim[ri_id][sim_index] 3013 else: 3014 r1[si, mi] = spins[si].ri_data[ri_id] 3015 3016 # Check the data to prevent user mistakes. 3017 for mi in range(field_count): 3018 # The frequency. 3019 frq = return_value_from_frq_index(mi=mi) 3020 3021 # Check for R1 data for this frequency. 3022 if not flags[mi]: 3023 raise RelaxError("R1 data for the %.1f MHz field strength cannot be found." % (frq/1e6)) 3024 3025 # Check the spin data. 3026 for si in range(spin_num): 3027 if r1[si, mi] == -1.0: 3028 raise RelaxError("R1 data for the '%s' spin at %.1f MHz field strength cannot be found." % (spin_ids[si], frq/1e6)) 3029 3030 # Return the data. 3031 return r1
3032 3033
3034 -def return_r2eff_arrays(spins=None, spin_ids=None, fields=None, field_count=None, sim_index=None):
3035 """Return numpy arrays of the R2eff/R1rho values and errors. 3036 3037 @keyword spins: The list of spin containers in the cluster. 3038 @type spins: list of SpinContainer instances 3039 @keyword spin_ids: The list of spin IDs for the cluster. 3040 @type spin_ids: list of str 3041 @keyword fields: The list of spectrometer field strengths. 3042 @type fields: list of float 3043 @keyword field_count: The number of spectrometer field strengths. This may not be equal to the length of the fields list as the user may not have set the field strength. 3044 @type field_count: int 3045 @keyword sim_index: The index of the simulation to return the data of. This should be None if the normal data is required. 3046 @type sim_index: None or int 3047 @return: The numpy array structures of the R2eff/R1rho values, errors, missing data, and corresponding Larmor frequencies. For each structure, the first dimension corresponds to the experiment types, the second the spins of a spin block, the third to the spectrometer field strength, and the fourth is the dispersion points. For the Larmor frequency structure, the fourth dimension is omitted. For R1rho-type data, an offset dimension is inserted between the spectrometer field strength and the dispersion points. 3048 @rtype: lists of numpy float arrays, lists of numpy float arrays, lists of numpy float arrays, numpy rank-2 int array 3049 """ 3050 3051 # The counts. 3052 exp_num = num_exp_types() 3053 spin_num = count_spins(spins) 3054 3055 # 1H MMQ flag. 3056 proton_mmq_flag = has_proton_mmq_cpmg() 3057 3058 # Initialise the data structures for the target function. 3059 exp_types = [] 3060 values = [] 3061 errors = [] 3062 missing = [] 3063 frqs = [] 3064 frqs_H = [] 3065 relax_times = [] 3066 for exp_type, ei in loop_exp(return_indices=True): 3067 values.append([]) 3068 errors.append([]) 3069 missing.append([]) 3070 frqs.append([]) 3071 frqs_H.append([]) 3072 relax_times.append([]) 3073 for si in range(spin_num): 3074 values[ei].append([]) 3075 errors[ei].append([]) 3076 missing[ei].append([]) 3077 frqs[ei].append([]) 3078 frqs_H[ei].append([]) 3079 for frq, mi in loop_frq(return_indices=True): 3080 values[ei][si].append([]) 3081 errors[ei][si].append([]) 3082 missing[ei][si].append([]) 3083 frqs[ei][si].append(0.0) 3084 frqs_H[ei][si].append(0.0) 3085 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 3086 values[ei][si][mi].append([]) 3087 errors[ei][si][mi].append([]) 3088 missing[ei][si][mi].append([]) 3089 for mi in range(field_count): 3090 relax_times[ei].append(None) 3091 3092 # Pack the R2eff/R1rho data. 3093 data_flag = False 3094 si = 0 3095 for spin_index in range(len(spins)): 3096 # Skip deselected spins. 3097 if not spins[spin_index].select: 3098 continue 3099 3100 # Alias the spin. 3101 spin = spins[spin_index] 3102 spin_id = spin_ids[spin_index] 3103 3104 # Get the attached proton. 3105 proton = None 3106 if proton_mmq_flag: 3107 # Get all protons. 3108 proton_spins = return_attached_protons(spin_id) 3109 3110 # Only one allowed. 3111 if len(proton_spins) > 1: 3112 raise RelaxError("Only one attached proton is supported for the MMQ-type models.") 3113 3114 # Missing proton. 3115 if not len(proton_spins): 3116 raise RelaxError("No proton attached to the spin '%s' could be found. This is required for the MMQ-type models." % spin_id) 3117 3118 # Alias the single proton. 3119 proton = proton_spins[0] 3120 3121 # No data. 3122 if not hasattr(spin, 'r2eff') and not hasattr(proton, 'r2eff'): 3123 continue 3124 data_flag = True 3125 3126 # No isotope information. 3127 if not hasattr(spin, 'isotope'): 3128 raise RelaxSpinTypeError(spin_id=spin_ids[si]) 3129 3130 # Loop over the R2eff data. 3131 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 3132 3133 # Alias the correct spin. 3134 current_spin = spin 3135 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 3136 current_spin = proton 3137 3138 # Add the experiment type. 3139 if exp_type not in exp_types: 3140 exp_types.append(exp_type) 3141 3142 # The key. 3143 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 3144 if mi == 0: 3145 fact = 60.83831274541046 3146 else: 3147 fact = 81.11775032721394 3148 3149 # The Larmor frequency for this spin (and that of an attached proton for the MMQ models) and field strength (in MHz*2pi to speed up the ppm to rad/s conversion). 3150 if frq != None: 3151 frqs[ei][si][mi] = 2.0 * pi * frq / g1H * return_gyromagnetic_ratio(spin.isotope) * 1e-6 3152 frqs_H[ei][si][mi] = 2.0 * pi * frq * 1e-6 3153 3154 # Missing data. 3155 if key not in current_spin.r2eff.keys(): 3156 values[ei][si][mi][oi].append(0.0) 3157 errors[ei][si][mi][oi].append(1.0) 3158 missing[ei][si][mi][oi].append(1) 3159 continue 3160 else: 3161 missing[ei][si][mi][oi].append(0) 3162 3163 # The values. 3164 if sim_index == None: 3165 values[ei][si][mi][oi].append(current_spin.r2eff[key]) 3166 else: 3167 values[ei][si][mi][oi].append(current_spin.r2eff_sim[sim_index][key]) 3168 3169 # The errors. 3170 errors[ei][si][mi][oi].append(current_spin.r2eff_err[key]) 3171 3172 # The relaxation times. 3173 for id in cdp.spectrum_ids: 3174 # Non-matching data. 3175 if cdp.spectrometer_frq[id] != frq: 3176 continue 3177 if cdp.exp_type[id] != exp_type: 3178 continue 3179 if exp_type in EXP_TYPE_LIST_CPMG: 3180 if id not in cdp.cpmg_frqs.keys() or cdp.cpmg_frqs[id] != point: 3181 continue 3182 else: 3183 if id not in cdp.spin_lock_nu1.keys() or cdp.spin_lock_nu1[id] != point: 3184 continue 3185 3186 # Found. 3187 relax_time = cdp.relax_times[id] 3188 break 3189 3190 # Check the value if already set. 3191 if relax_times[ei][mi] != None: 3192 if relax_times[ei][mi] != relax_time: 3193 raise RelaxError("The relaxation times do not match for all experiments.") 3194 continue 3195 3196 # Store the time. 3197 relax_times[ei][mi] = relax_time 3198 3199 # Increment the spin index. 3200 si += 1 3201 3202 # No R2eff/R1rho data for the spin cluster. 3203 if not data_flag: 3204 raise RelaxError("No R2eff/R1rho data could be found for the spin cluster %s." % spin_ids) 3205 3206 # Convert to numpy arrays. 3207 relax_times = array(relax_times, float64) 3208 for exp_type, ei in loop_exp(return_indices=True): 3209 for si in range(spin_num): 3210 for frq, mi in loop_frq(return_indices=True): 3211 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 3212 values[ei][si][mi][oi] = array(values[ei][si][mi][oi], float64) 3213 errors[ei][si][mi][oi] = array(errors[ei][si][mi][oi], float64) 3214 missing[ei][si][mi][oi] = array(missing[ei][si][mi][oi], int32) 3215 3216 # Return the structures. 3217 return values, errors, missing, frqs, frqs_H, exp_types, relax_times
3218 3219
3220 -def return_relax_times():
3221 """Return the list of relaxation times. 3222 3223 @return: The list of relaxation times in s. 3224 @rtype: numpy rank-2 float64 array 3225 """ 3226 3227 # No data. 3228 if not hasattr(cdp, 'relax_times'): 3229 return None 3230 3231 # Initialise. 3232 relax_times = zeros((count_exp(), count_frq()), float64) 3233 3234 # Loop over the experiment types. 3235 for exp_type, frq, point, time, ei, mi, di, ti in loop_exp_frq_point_time(return_indices=True): 3236 # Fetch all of the matching intensity keys. 3237 keys = find_intensity_keys(exp_type=exp_type, frq=frq, point=point, time=time, raise_error=False) 3238 3239 # No data. 3240 if not len(keys): 3241 continue 3242 3243 # Add the data. 3244 relax_times[ei][mi] = cdp.relax_times[keys[0]] 3245 3246 # Return the data. 3247 return relax_times
3248 3249
3250 -def return_spin_lock_nu1(ref_flag=True):
3251 """Return the list of spin-lock field strengths. 3252 3253 @keyword ref_flag: A flag which if False will cause the reference spectrum frequency of None to be removed from the list. 3254 @type ref_flag: bool 3255 @return: The list of spin-lock field strengths in Hz. It has the dimensions {Ei, Mi, Oi}. 3256 @rtype: rank-2 list of numpy rank-1 float64 arrays 3257 """ 3258 3259 # No data. 3260 if not hasattr(cdp, 'spin_lock_nu1_list'): 3261 return None 3262 3263 # Initialise. 3264 nu1 = [] 3265 3266 # First loop over the experiment types. 3267 for exp_type, ei in loop_exp(return_indices=True): 3268 # Add a new dimension. 3269 nu1.append([]) 3270 3271 # Then loop over the spectrometer frequencies. 3272 for frq, mi in loop_frq(return_indices=True): 3273 # Add a new dimension. 3274 nu1[ei].append([]) 3275 3276 # Loop over the offsets. 3277 for offset, oi in loop_offset(exp_type=exp_type, frq=frq, return_indices=True): 3278 # Add a new dimension. 3279 nu1[ei][mi].append([]) 3280 3281 # Loop over the fields. 3282 for point in cdp.spin_lock_nu1_list: 3283 # Skip reference points. 3284 if (not ref_flag) and point == None: 3285 continue 3286 3287 # Find a matching experiment ID. 3288 found = False 3289 for id in cdp.exp_type.keys(): 3290 # Skip non-matching experiments. 3291 if cdp.exp_type[id] != exp_type: 3292 continue 3293 3294 # Skip non-matching spectrometer frequencies. 3295 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 3296 continue 3297 3298 # Skip non-matching offsets. 3299 if offset != None and hasattr(cdp, 'spin_lock_offset') and cdp.spin_lock_offset[id] != offset: 3300 continue 3301 3302 # Skip non-matching points. 3303 if cdp.spin_lock_nu1[id] != point: 3304 continue 3305 3306 # Found. 3307 found = True 3308 break 3309 3310 # No data. 3311 if not found: 3312 continue 3313 3314 # Add the data. 3315 nu1[ei][mi][oi].append(point) 3316 3317 # Convert to a numpy array. 3318 nu1[ei][mi][oi] = array(nu1[ei][mi][oi], float64) 3319 3320 # Return the data. 3321 return nu1
3322 3323
3324 -def return_spin_lock_nu1_single(exp_type=None, frq=None, offset=None, ref_flag=True):
3325 """Return the list of spin-lock field strengths. 3326 3327 @keyword exp_type: The experiment type. 3328 @type exp_type: str 3329 @keyword frq: The spectrometer frequency in Hz. 3330 @type frq: float 3331 @keyword offset: The spin-lock offset. 3332 @type offset: None or float 3333 @keyword ref_flag: A flag which if False will cause the reference spectrum frequency of None to be removed from the list. 3334 @type ref_flag: bool 3335 @return: The list of spin-lock field strengths in Hz. 3336 @rtype: numpy rank-1 float64 array 3337 """ 3338 3339 # No data. 3340 if not hasattr(cdp, 'spin_lock_nu1_list'): 3341 return None 3342 3343 # Initialise. 3344 nu1 = [] 3345 3346 # Loop over the points. 3347 for point in cdp.spin_lock_nu1_list: 3348 # Skip reference points. 3349 if (not ref_flag) and point == None: 3350 continue 3351 3352 # Find a matching experiment ID. 3353 found = False 3354 for id in cdp.exp_type.keys(): 3355 # Skip non-matching experiments. 3356 if cdp.exp_type[id] != exp_type: 3357 continue 3358 3359 # Skip non-matching spectrometer frequencies. 3360 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 3361 continue 3362 3363 # Skip non-matching offsets. 3364 if offset != None and hasattr(cdp, 'spin_lock_offset') and cdp.spin_lock_offset[id] != offset: 3365 continue 3366 3367 # Skip non-matching points. 3368 if cdp.spin_lock_nu1[id] != point: 3369 continue 3370 3371 # Found. 3372 found = True 3373 break 3374 3375 # No data. 3376 if not found: 3377 continue 3378 3379 # Add the data. 3380 nu1.append(point) 3381 3382 # Return the data as a numpy array. 3383 return array(nu1, float64)
3384 3385
3386 -def return_value_from_frq_index(mi=None):
3387 """Return the spectrometer frequency corresponding to the frequency index. 3388 3389 @keyword mi: The spectrometer frequency index. 3390 @type mi: int 3391 @return: The spectrometer frequency in Hertz or None if no information is present. 3392 @rtype: float 3393 """ 3394 3395 # No data. 3396 if not hasattr(cdp, 'spectrometer_frq_list'): 3397 return None 3398 3399 # Return the field. 3400 return cdp.spectrometer_frq_list[mi]
3401 3402
3403 -def return_value_from_offset_index(ei=None, mi=None, oi=None):
3404 """Return the offset corresponding to the offset index. 3405 3406 @keyword ei: The experiment type index. 3407 @type ei: int 3408 @keyword mi: The spectrometer frequency index. 3409 @type mi: int 3410 @keyword oi: The offset index. 3411 @type oi: int 3412 @return: The offset in Hertz or None if no information is present. 3413 @rtype: float 3414 """ 3415 3416 # Checks. 3417 if ei == None: 3418 raise RelaxError("The experiment type index must be supplied.") 3419 if mi == None: 3420 raise RelaxError("The spectrometer frequency index must be supplied.") 3421 3422 # Initialise the index. 3423 new_oi = -1 3424 3425 # The experiment type and frequency. 3426 exp_type = cdp.exp_type_list[ei] 3427 frq = return_value_from_frq_index(mi) 3428 3429 # CPMG-type data. 3430 if exp_type in EXP_TYPE_LIST_CPMG: 3431 # Return zero until hard pulse offset handling is implemented. 3432 return 0.0 3433 3434 # R1rho-type data. 3435 if exp_type in EXP_TYPE_LIST_R1RHO: 3436 # No offsets. 3437 if not hasattr(cdp, 'spin_lock_offset'): 3438 return None 3439 3440 # Loop over the offset data. 3441 for offset in cdp.spin_lock_offset_list: 3442 # Increment the index. 3443 new_oi += 1 3444 3445 # Find a matching experiment ID. 3446 found = False 3447 for id in cdp.exp_type.keys(): 3448 # Skip non-matching experiments. 3449 if cdp.exp_type[id] != exp_type: 3450 continue 3451 3452 # Skip non-matching spectrometer frequencies. 3453 if hasattr(cdp, 'spectrometer_frq') and cdp.spectrometer_frq[id] != frq: 3454 continue 3455 3456 # Skip non-matching offsets. 3457 if new_oi != oi: 3458 continue 3459 3460 # Found. 3461 found = True 3462 break 3463 3464 # No data. 3465 if not found: 3466 continue 3467 3468 # Return the offset. 3469 return offset
3470 3471
3472 -def set_exp_type(spectrum_id=None, exp_type=None):
3473 """Select the relaxation dispersion experiment type performed. 3474 3475 @keyword spectrum_id: The spectrum ID string. 3476 @type spectrum_id: str 3477 @keyword exp: The relaxation dispersion experiment type. 3478 @type exp: str 3479 """ 3480 3481 # Data checks. 3482 pipes.test() 3483 3484 # Add the spectrum ID to the data store if needed. 3485 add_spectrum_id(spectrum_id) 3486 3487 # Check the experiment type. 3488 if exp_type not in EXP_TYPE_LIST: 3489 raise RelaxError("The relaxation dispersion experiment '%s' is invalid, it must be one of %s." % (exp_type, EXP_TYPE_LIST)) 3490 3491 # Initialise the experiment type data structures if needed. 3492 if not hasattr(cdp, 'exp_type'): 3493 cdp.exp_type = {} 3494 if not hasattr(cdp, 'exp_type_list'): 3495 cdp.exp_type_list = [] 3496 3497 # Store the value. 3498 cdp.exp_type[spectrum_id] = exp_type 3499 3500 # Unique experiments. 3501 if cdp.exp_type[spectrum_id] not in cdp.exp_type_list: 3502 cdp.exp_type_list.append(cdp.exp_type[spectrum_id]) 3503 3504 # Printout. 3505 text = "The spectrum ID '%s' is now set to " % spectrum_id 3506 if exp_type == EXP_TYPE_CPMG_SQ: 3507 text += EXP_TYPE_DESC_CPMG_SQ + "." 3508 elif exp_type == EXP_TYPE_CPMG_MQ: 3509 text += EXP_TYPE_DESC_CPMG_MQ + "." 3510 elif exp_type == EXP_TYPE_CPMG_DQ: 3511 text += EXP_TYPE_DESC_CPMG_DQ + "." 3512 elif exp_type == EXP_TYPE_CPMG_ZQ: 3513 text += EXP_TYPE_DESC_CPMG_ZQ + "." 3514 elif exp_type == EXP_TYPE_CPMG_PROTON_SQ: 3515 text += EXP_TYPE_DESC_CPMG_PROTON_SQ + "." 3516 elif exp_type == EXP_TYPE_CPMG_PROTON_MQ: 3517 text += EXP_TYPE_DESC_CPMG_PROTON_MQ + "." 3518 elif exp_type == EXP_TYPE_R1RHO: 3519 text += EXP_TYPE_DESC_R1RHO + "." 3520 print(text)
3521 3522
3523 -def spin_has_frq_data(spin=None, frq=None):
3524 """Determine if the spin has intensity data for the given spectrometer frequency. 3525 3526 @keyword spin: The specific spin data container. 3527 @type spin: SpinContainer instance 3528 @keyword frq: The spectrometer frequency. 3529 @type frq: float 3530 @return: True if data for that spectrometer frequency is present, False otherwise. 3531 @rtype: bool 3532 """ 3533 3534 # Loop over the intensity data. 3535 for key in spin.intensities.keys(): 3536 if key in cdp.spectrometer_frq and cdp.spectrometer_frq[key] == frq: 3537 return True 3538 3539 # No data. 3540 return False
3541 3542
3543 -def spin_ids_to_containers(spin_ids):
3544 """Take the list of spin IDs and return the corresponding spin containers. 3545 3546 This is useful for handling the data from the model_loop() method. 3547 3548 3549 @param spin_ids: The list of spin ID strings. 3550 @type spin_ids: list of str 3551 @return: The list of spin containers. 3552 @rtype: list of SpinContainer instances 3553 """ 3554 3555 # Loop over the IDs and fetch the container. 3556 spins = [] 3557 for id in spin_ids: 3558 spins.append(return_spin(id)) 3559 3560 # Return the containers. 3561 return spins
3562 3563
3564 -def spin_lock_field(spectrum_id=None, field=None):
3565 """Set the spin-lock field strength (nu1) for the given spectrum. 3566 3567 @keyword spectrum_id: The spectrum ID string. 3568 @type spectrum_id: str 3569 @keyword field: The spin-lock field strength (nu1) in Hz. 3570 @type field: int or float 3571 """ 3572 3573 # Test if the spectrum ID exists. 3574 if spectrum_id not in cdp.spectrum_ids: 3575 raise RelaxNoSpectraError(spectrum_id) 3576 3577 # Initialise the global nu1 data structures if needed. 3578 if not hasattr(cdp, 'spin_lock_nu1'): 3579 cdp.spin_lock_nu1 = {} 3580 if not hasattr(cdp, 'spin_lock_nu1_list'): 3581 cdp.spin_lock_nu1_list = [] 3582 3583 # Add the frequency, converting to a float if needed. 3584 if field == None: 3585 cdp.spin_lock_nu1[spectrum_id] = field 3586 else: 3587 cdp.spin_lock_nu1[spectrum_id] = float(field) 3588 3589 # The unique curves for the R2eff fitting (R1rho). 3590 if cdp.spin_lock_nu1[spectrum_id] not in cdp.spin_lock_nu1_list: 3591 cdp.spin_lock_nu1_list.append(cdp.spin_lock_nu1[spectrum_id]) 3592 3593 # Sort the list (handling None for Python 3). 3594 flag = False 3595 if None in cdp.spin_lock_nu1_list: 3596 cdp.spin_lock_nu1_list.pop(cdp.spin_lock_nu1_list.index(None)) 3597 flag = True 3598 cdp.spin_lock_nu1_list.sort() 3599 if flag: 3600 cdp.spin_lock_nu1_list.insert(0, None) 3601 3602 # Update the exponential curve count (skipping the reference if present). 3603 cdp.dispersion_points = len(cdp.spin_lock_nu1_list) 3604 if None in cdp.spin_lock_nu1_list: 3605 cdp.dispersion_points -= 1 3606 3607 # Printout. 3608 if field == None: 3609 print("The spectrum ID '%s' is set to the reference." % spectrum_id) 3610 else: 3611 print("The spectrum ID '%s' spin-lock field strength is set to %s kHz." % (spectrum_id, cdp.spin_lock_nu1[spectrum_id]/1000.0))
3612 3613
3614 -def spin_lock_offset(spectrum_id=None, offset=None):
3615 """Set the spin-lock offset (omega_rf) for the given spectrum. 3616 3617 @keyword spectrum_id: The spectrum ID string. 3618 @type spectrum_id: str 3619 @keyword offset: The spin-lock offset (omega_rf) in ppm. 3620 @type offset: int or float 3621 """ 3622 3623 # Test if the spectrum ID exists. 3624 if spectrum_id not in cdp.spectrum_ids: 3625 raise RelaxNoSpectraError(spectrum_id) 3626 3627 # Initialise the global offset data structures if needed. 3628 if not hasattr(cdp, 'spin_lock_offset'): 3629 cdp.spin_lock_offset = {} 3630 if not hasattr(cdp, 'spin_lock_offset_list'): 3631 cdp.spin_lock_offset_list = [] 3632 3633 # Add the offset, converting to a float if needed. 3634 if offset == None: 3635 raise RelaxError("The offset value must be provided.") 3636 cdp.spin_lock_offset[spectrum_id] = float(offset) 3637 3638 # The unique curves for the R2eff fitting (R1rho). 3639 if cdp.spin_lock_offset[spectrum_id] not in cdp.spin_lock_offset_list: 3640 cdp.spin_lock_offset_list.append(cdp.spin_lock_offset[spectrum_id]) 3641 3642 # Sort the list. 3643 cdp.spin_lock_offset_list.sort() 3644 3645 # Printout. 3646 print("Setting the '%s' spectrum spin-lock offset to %s ppm." % (spectrum_id, cdp.spin_lock_offset[spectrum_id]))
3647 3648
3649 -def write_disp_curves(dir=None, force=None):
3650 """Write out the dispersion curves to text files. 3651 3652 One file will be created per spin system. 3653 3654 3655 @keyword dir: The optional directory to place the file into. 3656 @type dir: str 3657 @param force: If True, the files will be overwritten if they already exists. 3658 @type force: bool 3659 """ 3660 3661 # Checks. 3662 pipes.test() 3663 check_mol_res_spin_data() 3664 3665 # The formatting strings. 3666 format_head = "# %-18s %-20s %-20s %-20s %-20s %-20s\n" 3667 format = "%-20s %20s %20s %20s %20s %20s\n" 3668 3669 # 1H MMQ flag. 3670 proton_mmq_flag = has_proton_mmq_cpmg() 3671 3672 # Loop over each spin. 3673 for spin, spin_id in spin_loop(return_id=True, skip_desel=True): 3674 # Skip protons for MMQ data. 3675 if spin.model in MODEL_LIST_MMQ and spin.isotope == '1H': 3676 continue 3677 3678 # Define writing variables. 3679 writing_vars = [['disp',("Experiment_name", "Field_strength_(MHz)", "Disp_point_(Hz)", "R2eff_(measured)", "R2eff_(back_calc)", "R2eff_errors")]] 3680 3681 # If the model is of R1rho type, then also write as R2eff as function of theta. 3682 if spin.model in MODEL_LIST_R1RHO_FULL and has_r1rho_exp_type() and hasattr(spin, 'isotope'): 3683 # Add additonal looping over writing parameters. 3684 writing_vars.append(['disp_theta',("Experiment_name", "Field_strength_(MHz)", "Tilt_angle_(rad)", "R2eff_(measured)", "R2eff_(back_calc)", "R2eff_errors")]) 3685 #writing_vars.append(['disp_w_eff',("Experiment_name", "Field_strength_(MHz)", "Effective_field_(rad_s-1))", "R2eff_(measured)", "R2eff_(back_calc)", "R2eff_errors")]) 3686 3687 # Loop over writing vars 3688 for wvar in writing_vars: 3689 # Get the attached proton. 3690 proton = None 3691 if proton_mmq_flag: 3692 proton = return_attached_protons(spin_id)[0] 3693 3694 # The unique file name. 3695 file_name = "%s%s.out" % (wvar[0], spin_id.replace('#', '_').replace(':', '_').replace('@', '_')) 3696 3697 # Open the file for writing. 3698 file_path = get_file_path(file_name, dir) 3699 file = open_write_file(file_name, dir, force) 3700 3701 # Write a header. 3702 file.write(format_head % wvar[1]) 3703 3704 # Loop over the dispersion points. 3705 for exp_type, frq, offset, point, ei, mi, oi, di in loop_exp_frq_offset_point(return_indices=True): 3706 # Alias the correct spin. 3707 current_spin = spin 3708 if exp_type in [EXP_TYPE_CPMG_PROTON_SQ, EXP_TYPE_CPMG_PROTON_MQ]: 3709 current_spin = proton 3710 3711 # The data key. 3712 key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) 3713 3714 # Format the R2eff data. 3715 r2eff = "-" 3716 if hasattr(current_spin, 'r2eff') and key in current_spin.r2eff: 3717 r2eff = "%.15f" % current_spin.r2eff[key] 3718 3719 # Format the R2eff back calc data. 3720 r2eff_bc = "-" 3721 if hasattr(current_spin, 'r2eff_bc') and key in current_spin.r2eff_bc: 3722 r2eff_bc = "%.15f" % current_spin.r2eff_bc[key] 3723 3724 # Format the R2eff errors. 3725 r2eff_err = "-" 3726 if hasattr(current_spin, 'r2eff_err') and key in current_spin.r2eff_err: 3727 r2eff_err = "%.15f" % current_spin.r2eff_err[key] 3728 3729 # Define value to be written. 3730 if wvar[0] == 'disp_theta': 3731 theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=spin) 3732 value = theta_spin_dic[key] 3733 elif wvar[0] == 'disp_w_eff': 3734 theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=spin) 3735 value = w_eff_spin_dic[key] 3736 # Else use the standard dispersion point data. 3737 else: 3738 value = point 3739 3740 # Write out the data. 3741 frq_text = "%.9f" % (frq/1e6) 3742 value_text = "%.6f" % value 3743 file.write(format % (repr(exp_type), frq_text, value_text, r2eff, r2eff_bc, r2eff_err)) 3744 3745 # Close the file. 3746 file.close() 3747 3748 # Add the file to the results file list. 3749 add_result_file(type='text', label='Text', file=file_path)
3750