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

Source Code for Module specific_analyses.model_free.api

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2007-2009,2011-2012,2014,2017 Edward d'Auvergne               # 
   4  # Copyright (C) 2007 Gary S Thompson                                          # 
   5  #                                                                             # 
   6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
   7  #                                                                             # 
   8  # This program is free software: you can redistribute it and/or modify        # 
   9  # it under the terms of the GNU General Public License as published by        # 
  10  # the Free Software Foundation, either version 3 of the License, or           # 
  11  # (at your option) any later version.                                         # 
  12  #                                                                             # 
  13  # This program is distributed in the hope that it will be useful,             # 
  14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
  15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
  16  # GNU General Public License for more details.                                # 
  17  #                                                                             # 
  18  # You should have received a copy of the GNU General Public License           # 
  19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
  20  #                                                                             # 
  21  ############################################################################### 
  22   
  23  # Module docstring. 
  24  """The Lipari-Szabo model-free analysis API object.""" 
  25   
  26   
  27  # Python module imports. 
  28  import bmrblib 
  29  from copy import deepcopy 
  30  from math import pi 
  31  from minfx.grid import grid_split 
  32  from numpy import array, dot, float64, int32, zeros 
  33  from numpy.linalg import inv 
  34  from re import match, search 
  35  import string 
  36  import sys 
  37  from types import MethodType 
  38  from warnings import warn 
  39   
  40  # relax module imports. 
  41  from lib.arg_check import is_num_list, is_str_list 
  42  from lib.errors import RelaxError, RelaxFault, RelaxNoModelError, RelaxNoSequenceError, RelaxNoTensorError 
  43  from lib.float import isInf 
  44  from lib.periodic_table import periodic_table 
  45  from lib.physical_constants import h_bar, mu0 
  46  from lib.text.sectioning import subsection 
  47  from lib.warnings import RelaxDeselectWarning, RelaxWarning 
  48  from multi import Processor_box 
  49  from pipe_control import diffusion_tensor, interatomic, mol_res_spin, pipes, relax_data, sequence 
  50  from pipe_control.bmrb import list_sample_conditions 
  51  from pipe_control.exp_info import bmrb_write_citations, bmrb_write_methods, bmrb_write_software 
  52  from pipe_control.interatomic import return_interatom_list 
  53  from pipe_control.mol_res_spin import count_spins, exists_mol_res_spin_data, find_index, get_molecule_names, return_spin, return_spin_from_index, return_spin_indices, spin_loop 
  54  from pipe_control.spectrometer import check_spectrometer_setup 
  55  from specific_analyses.api_base import API_base 
  56  from specific_analyses.api_common import API_common 
  57  from specific_analyses.model_free.bmrb import sf_csa_read, sf_model_free_read, to_bmrb_model 
  58  from specific_analyses.model_free.data import compare_objects 
  59  from specific_analyses.model_free.molmol import Molmol 
  60  from specific_analyses.model_free.model import determine_model_type 
  61  from specific_analyses.model_free.parameters import are_mf_params_set, assemble_param_names, assemble_param_vector, linear_constraints 
  62  from specific_analyses.model_free.optimisation import MF_grid_command, MF_memo, MF_minimise_command, minimise_data_setup, relax_data_opt_structs 
  63  from specific_analyses.model_free.parameter_object import Model_free_params 
  64  from specific_analyses.model_free.pymol import Pymol 
  65  from target_functions.mf import Mf 
  66   
  67   
68 -class Model_free(API_base, API_common):
69 """Parent class containing all the model-free specific functions.""" 70 71 # Class variable for storing the class instance (for the singleton design pattern). 72 instance = None 73
74 - def __init__(self):
75 """Initialise the class by placing API_common methods into the API.""" 76 77 # Place methods into the API. 78 self.base_data_loop = self._base_data_loop_spin 79 self.return_error = self._return_error_relax_data 80 self.return_value = self._return_value_general 81 self.sim_pack_data = self._sim_pack_relax_data 82 83 # Initialise the macro classes. 84 self._molmol_macros = Molmol() 85 self._pymol_macros = Pymol() 86 87 # Alias the macro creation methods. 88 self.pymol_macro = self._pymol_macros.create_macro 89 self.molmol_macro = self._molmol_macros.create_macro 90 91 # Place a copy of the parameter list object in the instance namespace. 92 self._PARAMS = Model_free_params()
93 94
95 - def back_calc_ri(self, spin_index=None, ri_id=None, ri_type=None, frq=None):
96 """Back-calculation of relaxation data from the model-free parameter values. 97 98 @keyword spin_index: The global spin index. 99 @type spin_index: int 100 @keyword ri_id: The relaxation data ID string. 101 @type ri_id: str 102 @keyword ri_type: The relaxation data type. 103 @type ri_type: str 104 @keyword frq: The field strength. 105 @type frq: float 106 @return: The back calculated relaxation data value corresponding to the index. 107 @rtype: float 108 """ 109 110 # Get the spin container. 111 spin, spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 112 113 # Missing interatomic vectors. 114 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'): 115 interatoms = interatomic.return_interatom_list(spin_hash=spin._hash) 116 for i in range(len(interatoms)): 117 # No dipolar relaxation mechanism. 118 if not interatoms[i].dipole_pair: 119 continue 120 121 # Check the unit vectors. 122 if not hasattr(interatoms[i], 'vector') or interatoms[i].vector is None: 123 warn(RelaxDeselectWarning(spin_id, 'missing structural data')) 124 return 125 126 # Execute the over-fit deselection. 127 self.overfit_deselect(data_check=False, verbose=False) 128 129 # Get the relaxation value from the minimise function. 130 value = self.minimise(min_algor='back_calc', min_options=(spin_index, ri_id, ri_type, frq)) 131 132 # Return the relaxation value. 133 return value
134 135
136 - def bmrb_read(self, file_path, version=None, sample_conditions=None):
137 """Read the model-free results from a BMRB NMR-STAR v3.1 formatted file. 138 139 @param file_path: The full file path. 140 @type file_path: str 141 @keyword version: The BMRB version to force the reading. 142 @type version: None or str 143 @keyword sample_conditions: The sample condition label to read. Only one sample condition can be read per data pipe. 144 @type sample_conditions: None or str 145 """ 146 147 # Initialise the NMR-STAR data object. 148 star = bmrblib.create_nmr_star('relax_model_free_results', file_path, version) 149 150 # Read the contents of the STAR formatted file. 151 star.read() 152 153 # The sample conditions. 154 sample_conds = list_sample_conditions(star) 155 if sample_conditions and sample_conditions not in sample_conds: 156 raise RelaxError("The sample conditions label '%s' does not correspond to any of the labels in the file: %s" % (sample_conditions, sample_conds)) 157 if not sample_conditions and len(sample_conds) > 1: 158 raise RelaxError("Only one of the sample conditions in %s can be loaded per relax data pipe." % sample_conds) 159 160 # The diffusion tensor. 161 diffusion_tensor.bmrb_read(star) 162 163 # Generate the molecule and residue containers from the entity records. 164 mol_res_spin.bmrb_read(star) 165 166 # Read the relaxation data saveframes. 167 relax_data.bmrb_read(star, sample_conditions=sample_conditions) 168 169 # Read the model-free data saveframes. 170 sf_model_free_read(star, sample_conditions=sample_conditions) 171 172 # Read the CSA data saveframes. 173 sf_csa_read(star)
174 175
176 - def bmrb_write(self, file_path, version=None):
177 """Write the model-free results to a BMRB NMR-STAR v3.1 formatted file. 178 179 @param file_path: The full file path. 180 @type file_path: str 181 @keyword version: The BMRB NMR-STAR dictionary format to output to. 182 @type version: str 183 """ 184 185 # Checks. 186 check_spectrometer_setup(escalate=2) 187 188 # Alias the current data pipe. 189 cdp = pipes.get_pipe() 190 191 # Initialise the NMR-STAR data object. 192 star = bmrblib.create_nmr_star('relax_model_free_results', file_path, version) 193 194 # Global minimisation stats. 195 global_chi2 = None 196 if hasattr(cdp, 'chi2'): 197 global_chi2 = cdp.chi2 198 199 # Rex frq. 200 rex_frq = cdp.spectrometer_frq[cdp.ri_ids[0]] 201 202 # Initialise the spin specific data lists. 203 mol_name_list = [] 204 res_num_list = [] 205 res_name_list = [] 206 atom_name_list = [] 207 208 csa_list = [] 209 r_list = [] 210 isotope_list = [] 211 element_list = [] 212 213 local_tm_list = [] 214 s2_list = [] 215 s2f_list = [] 216 s2s_list = [] 217 te_list = [] 218 tf_list = [] 219 ts_list = [] 220 rex_list = [] 221 222 local_tm_err_list = [] 223 s2_err_list = [] 224 s2f_err_list = [] 225 s2s_err_list = [] 226 te_err_list = [] 227 tf_err_list = [] 228 ts_err_list = [] 229 rex_err_list = [] 230 231 chi2_list = [] 232 model_list = [] 233 234 # Store the spin specific data in lists for later use. 235 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, skip_desel=True, return_id=True): 236 # Skip the protons. 237 if spin.name == 'H' or (hasattr(spin, 'element') and spin.element == 'H'): 238 warn(RelaxWarning("Skipping the proton spin '%s'." % spin_id)) 239 continue 240 241 # No model setup. 242 if not hasattr(spin, 'model'): 243 warn(RelaxWarning("Skipping the spin '%s' as no model-free model has been set up." % spin_id)) 244 continue 245 246 # Check the data for None (not allowed in BMRB!). 247 if res_num == None: 248 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id) 249 if res_name == None: 250 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id) 251 if spin.name == None: 252 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id) 253 if not hasattr(spin, 'isotope') or spin.isotope == None: 254 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id) 255 if not hasattr(spin, 'element') or spin.element == None: 256 raise RelaxError("For the BMRB, the spin element type of '%s' must be specified. Please use the spin user function for setting the element type." % spin_id) 257 258 # The molecule/residue/spin info. 259 mol_name_list.append(mol_name) 260 res_num_list.append(res_num) 261 res_name_list.append(res_name) 262 atom_name_list.append(spin.name) 263 264 # CSA values. 265 if hasattr(spin, 'csa'): 266 csa_list.append(spin.csa * 1e6) # In ppm. 267 else: 268 csa_list.append(None) 269 270 # Interatomic distances. 271 interatoms = return_interatom_list(spin_hash=spin._hash) 272 for i in range(len(interatoms)): 273 # No relaxation mechanism. 274 if not interatoms[i].dipole_pair: 275 continue 276 277 # Add the interatomic distance. 278 if hasattr(interatoms[i], 'r'): 279 r_list.append(interatoms[i].r) 280 else: 281 r_list.append(None) 282 283 # Stop adding. 284 break 285 286 # The nuclear isotope. 287 if hasattr(spin, 'isotope'): 288 isotope_list.append(int(spin.isotope.strip(string.ascii_letters))) 289 else: 290 isotope_list.append(None) 291 292 # The element. 293 if hasattr(spin, 'element'): 294 element_list.append(spin.element) 295 else: 296 element_list.append(None) 297 298 # Diffusion tensor. 299 if hasattr(spin, 'local_tm'): 300 local_tm_list.append(spin.local_tm) 301 else: 302 local_tm_list.append(None) 303 if hasattr(spin, 'local_tm_err'): 304 local_tm_err_list.append(spin.local_tm_err) 305 else: 306 local_tm_err_list.append(None) 307 308 # Model-free parameter values. 309 s2_list.append(spin.s2) 310 s2f_list.append(spin.s2f) 311 s2s_list.append(spin.s2s) 312 te_list.append(spin.te) 313 tf_list.append(spin.tf) 314 ts_list.append(spin.ts) 315 if spin.rex == None: 316 rex_list.append(None) 317 else: 318 rex_list.append(spin.rex / (2.0*pi*rex_frq**2)) 319 320 # Model-free parameter errors. 321 params = ['s2', 's2f', 's2s', 'te', 'tf', 'ts', 'rex'] 322 for param in params: 323 # The error list. 324 err_list = locals()[param+'_err_list'] 325 326 # Append the error. 327 if hasattr(spin, param+'_err'): 328 # The value. 329 val = getattr(spin, param+'_err') 330 331 # Scaling. 332 if param == 'rex' and val != None: 333 val = val / (2.0*pi*rex_frq**2) 334 335 # Append. 336 err_list.append(val) 337 338 # Or None. 339 else: 340 err_list.append(None) 341 342 343 # Opt stats. 344 chi2_list.append(spin.chi2) 345 346 # Model-free model. 347 model_list.append(to_bmrb_model(spin.model)) 348 349 # Check that spin data is present. 350 if not len(model_list): 351 raise RelaxError("No model-free data could be found for any of the currently selected spins.") 352 353 # Convert the molecule names into the entity IDs. 354 entity_ids = zeros(len(mol_name_list), int32) 355 mol_names = get_molecule_names() 356 for i in range(len(mol_name_list)): 357 for j in range(len(mol_names)): 358 if mol_name_list[i] == mol_names[j]: 359 entity_ids[i] = j+1 360 361 362 # Create Supergroup 2 : The citations. 363 ###################################### 364 365 # Generate the citations saveframe. 366 bmrb_write_citations(star) 367 368 369 # Create Supergroup 3 : The molecular assembly saveframes. 370 ########################################################## 371 372 # Generate the entity saveframe. 373 mol_res_spin.bmrb_write_entity(star) 374 375 376 # Create Supergroup 4: The experimental descriptions saveframes. 377 ################################################################# 378 379 # Generate the method saveframes. 380 bmrb_write_methods(star) 381 382 # Generate the software saveframe. 383 software_ids, software_labels = bmrb_write_software(star) 384 385 386 # Create Supergroup 5 : The NMR parameters saveframes. 387 ###################################################### 388 389 # Generate the CSA saveframe. 390 star.chem_shift_anisotropy.add(entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, csa=csa_list) 391 392 393 # Create Supergroup 6 : The kinetic data saveframes. 394 #################################################### 395 396 # Generate the relaxation data saveframes. 397 relax_data.bmrb_write(star) 398 399 400 # Create Supergroup 7 : The thermodynamics saveframes. 401 ###################################################### 402 403 # Get the relax software id. 404 for i in range(len(software_ids)): 405 if software_labels[i] == 'relax': 406 software_id = software_ids[i] 407 408 # Generate the model-free data saveframe. 409 star.model_free.add(global_chi2=global_chi2, software_ids=[software_id], software_labels=['relax'], entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, bond_length=r_list, local_tc=local_tm_list, s2=s2_list, s2f=s2f_list, s2s=s2s_list, te=te_list, tf=tf_list, ts=ts_list, rex=rex_list, local_tc_err=local_tm_err_list, s2_err=s2_err_list, s2f_err=s2f_err_list, s2s_err=s2s_err_list, te_err=te_err_list, tf_err=tf_err_list, ts_err=ts_err_list, rex_err=rex_err_list, rex_frq=rex_frq, chi2=chi2_list, model_fit=model_list) 410 411 412 # Create Supergroup 8 : The structure determination saveframes. 413 ############################################################### 414 415 # Generate the diffusion tensor saveframes. 416 diffusion_tensor.bmrb_write(star) 417 418 419 # Write the contents to the STAR formatted file. 420 star.write()
421 422
423 - def calculate(self, spin_id=None, scaling_matrix=None, verbosity=1, sim_index=None):
424 """Calculation of the model-free chi-squared value. 425 426 @keyword spin_id: The spin identification string. 427 @type spin_id: str 428 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 429 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 430 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 431 @type verbosity: int 432 @keyword sim_index: The optional MC simulation index. 433 @type sim_index: int 434 """ 435 436 # Test if sequence data is loaded. 437 if not exists_mol_res_spin_data(): 438 raise RelaxNoSequenceError 439 440 # Determine the model type. 441 model_type = determine_model_type() 442 443 # Test if diffusion tensor data exists. 444 if model_type != 'local_tm' and not diffusion_tensor.diff_data_exists(): 445 raise RelaxNoTensorError('diffusion') 446 447 # Test if the PDB file has been loaded. 448 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(cdp, 'structure'): 449 raise RelaxNoPdbError 450 451 # Loop over the spins. 452 for spin, id in spin_loop(spin_id, return_id=True): 453 # Skip deselected spins. 454 if not spin.select: 455 continue 456 457 # Test if the model-free model has been setup. 458 if not spin.model: 459 raise RelaxNoModelError 460 461 # Test if the nuclear isotope type has been set. 462 if not hasattr(spin, 'isotope'): 463 raise RelaxSpinTypeError 464 465 # Test if the model-free parameter values exist. 466 unset_param = are_mf_params_set(spin) 467 if unset_param != None: 468 raise RelaxNoValueError(unset_param) 469 470 # Test if the CSA value has been set. 471 if not hasattr(spin, 'csa') or spin.csa == None: 472 raise RelaxNoValueError("CSA") 473 474 # Test the interatomic data. 475 interatoms = return_interatom_list(spin_hash=spin._hash) 476 for interatom in interatoms: 477 # No relaxation mechanism. 478 if not interatom.dipole_pair: 479 continue 480 481 # Test if unit vectors exist. 482 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and not hasattr(interatom, 'vector'): 483 raise RelaxNoVectorsError 484 485 # Test if multiple unit vectors exist. 486 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere' and hasattr(interatom, 'vector') and is_num_list(interatom.vector[0], raise_error=False): 487 raise RelaxMultiVectorError 488 489 # The interacting spin. 490 if spin_id != interatom.spin_id1: 491 spin_id2 = interatom.spin_id1 492 else: 493 spin_id2 = interatom.spin_id2 494 spin2 = return_spin(spin_id=spin_id2) 495 496 # Test if the nuclear isotope type has been set. 497 if not hasattr(spin2, 'isotope'): 498 raise RelaxSpinTypeError 499 500 # Test if the interatomic distance has been set. 501 if not hasattr(interatom, 'r') or interatom.r == None: 502 raise RelaxNoValueError("interatomic distance", spin_id=spin_id, spin_id2=spin_id2) 503 504 # Skip spins where there is no data or errors. 505 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 506 continue 507 508 # Make sure that the errors are strictly positive numbers. 509 for ri_id in cdp.ri_ids: 510 # Alias. 511 err = spin.ri_data_err[ri_id] 512 513 # Checks. 514 if err != None and err == 0.0: 515 raise RelaxError("Zero error for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (id, ri_id)) 516 elif err != None and err < 0.0: 517 raise RelaxError("Negative error of %s for spin '%s' for the relaxation data ID '%s', minimisation not possible." % (err, id, ri_id)) 518 519 # Create the initial parameter vector. 520 param_vector = assemble_param_vector(spin=spin, sim_index=sim_index) 521 522 # The relaxation data optimisation structures. 523 data = relax_data_opt_structs(spin, sim_index=sim_index) 524 525 # The spin data. 526 ri_data = [array(data[0])] 527 ri_data_err = [array(data[1])] 528 num_frq = [data[2]] 529 num_ri = [data[3]] 530 ri_labels = [data[4]] 531 frq = [data[5]] 532 remap_table = [data[6]] 533 noe_r1_table = [data[7]] 534 gx = [periodic_table.gyromagnetic_ratio(spin.isotope)] 535 if sim_index == None: 536 csa = [spin.csa] 537 else: 538 csa = [spin.csa_sim[sim_index]] 539 540 # The interatomic data. 541 interatoms = return_interatom_list(spin_hash=spin._hash) 542 for i in range(len(interatoms)): 543 # No relaxation mechanism. 544 if not interatoms[i].dipole_pair: 545 continue 546 547 # The surrounding spins. 548 if spin_id != interatoms[i].spin_id1: 549 spin_id2 = interatoms[i].spin_id1 550 else: 551 spin_id2 = interatoms[i].spin_id2 552 spin2 = return_spin(spin_id=spin_id2) 553 554 # The data. 555 if sim_index == None: 556 r = [interatoms[i].r] 557 else: 558 r = [interatoms[i].r_sim[sim_index]] 559 560 # Vectors. 561 if model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 562 xh_unit_vectors = [interatoms[i].vector] 563 else: 564 xh_unit_vectors = [None] 565 566 # Gyromagnetic ratios. 567 gh = [periodic_table.gyromagnetic_ratio(spin2.isotope)] 568 569 # Count the number of model-free parameters for the residue index. 570 num_params = [len(spin.params)] 571 572 # Repackage the parameter values as a local model (ignore if the diffusion tensor is not fixed). 573 param_values = [assemble_param_vector(model_type='mf')] 574 575 # Package the diffusion tensor parameters. 576 if model_type == 'local_tm': 577 diff_params = [spin.local_tm] 578 diff_type = 'sphere' 579 else: 580 # Diff type. 581 diff_type = cdp.diff_tensor.type 582 583 # Spherical diffusion. 584 if diff_type == 'sphere': 585 diff_params = [cdp.diff_tensor.tm] 586 587 # Spheroidal diffusion. 588 elif diff_type == 'spheroid': 589 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.theta, cdp.diff_tensor.phi] 590 591 # Ellipsoidal diffusion. 592 elif diff_type == 'ellipsoid': 593 diff_params = [cdp.diff_tensor.tm, cdp.diff_tensor.Da, cdp.diff_tensor.Dr, cdp.diff_tensor.alpha, cdp.diff_tensor.beta, cdp.diff_tensor.gamma] 594 595 # Initialise the model-free function. 596 mf = Mf(init_params=param_vector, model_type='mf', diff_type=diff_type, diff_params=diff_params, num_spins=1, equations=[spin.equation], param_types=[spin.params], param_values=param_values, relax_data=ri_data, errors=ri_data_err, bond_length=r, csa=csa, num_frq=num_frq, frq=frq, num_ri=num_ri, remap_table=remap_table, noe_r1_table=noe_r1_table, ri_labels=ri_labels, gx=gx, gh=gh, h_bar=h_bar, mu0=mu0, num_params=num_params, vectors=xh_unit_vectors) 597 598 # Chi-squared calculation. 599 try: 600 chi2 = mf.func(param_vector) 601 except OverflowError: 602 chi2 = 1e200 603 604 # Global chi-squared value. 605 if model_type == 'all' or model_type == 'diff': 606 cdp.chi2 = cdp.chi2 + chi2 607 else: 608 spin.chi2 = chi2
609 610
611 - def create_mc_data(self, data_id=None):
612 """Create the Monte Carlo Ri data. 613 614 @keyword data_id: The spin identification string, as yielded by the base_data_loop() generator method. 615 @type data_id: str 616 @return: The Monte Carlo simulation data. 617 @rtype: list of floats 618 """ 619 620 # Initialise the MC data structure. 621 mc_data = [] 622 623 # Get the spin container and global spin index. 624 spin = return_spin(spin_id=data_id) 625 global_index = find_index(data_id) 626 627 # Skip deselected spins. 628 if not spin.select: 629 return 630 631 # Test if the model is set, returning None to signal that the spin should be skipped. 632 if not hasattr(spin, 'model') or not spin.model: 633 return None 634 635 # Loop over the relaxation data. 636 for ri_id in cdp.ri_ids: 637 # Back calculate the value. 638 value = self.back_calc_ri(spin_index=global_index, ri_id=ri_id, ri_type=cdp.ri_type[ri_id], frq=cdp.spectrometer_frq[ri_id]) 639 640 # Append the value. 641 mc_data.append(value) 642 643 # Return the data. 644 return mc_data
645 646
647 - def data_init(self, data, sim=False):
648 """Initialise the spin specific data structures. 649 650 @param data: The spin ID string from the _base_data_loop_spin() method. 651 @type data: str 652 @keyword sim: The Monte Carlo simulation flag, which if true will initialise the simulation data structure. 653 @type sim: bool 654 """ 655 656 # Get the spin container. 657 spin = return_spin(spin_id=data) 658 659 # Loop over the data structure names. 660 for name in self._PARAMS.loop(scope='spin'): 661 # Blacklisted data structures. 662 if name in ['ri_data', 'ri_data_bc', 'ri_data_err']: 663 continue 664 665 # Data structures which are initially empty arrays. 666 list_data = [ 'params' ] 667 if name in list_data: 668 init_data = [] 669 670 # Set everything else initially to None or False. 671 init_data = None 672 if self._PARAMS.type(name) == bool: 673 init_data = False 674 if name == 'select': 675 init_data = True 676 677 # If the name is not in the spin container, add it. 678 if not hasattr(spin, name): 679 setattr(spin, name, init_data)
680 681
682 - def deselect(self, sim_index=None, model_info=None):
683 """Deselect models or simulations. 684 685 @keyword sim_index: The optional Monte Carlo simulation index. If None, then models will be deselected, otherwise the given simulation will. 686 @type sim_index: None or int 687 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 688 @type model_info: int 689 """ 690 691 # Determine the model type. 692 model_type = determine_model_type() 693 694 # Local models. 695 if model_type == 'mf' or model_type == 'local_tm': 696 # Get the spin. 697 spin = return_spin_from_index(global_index=model_info) 698 699 # Spin deselection. 700 if sim_index == None: 701 spin.select = False 702 703 # Simulation deselection. 704 else: 705 spin.select_sim[sim_index] = False 706 707 # Global models. 708 else: 709 # Global model deselection. 710 if sim_index == None: 711 raise RelaxError("Cannot deselect the global model.") 712 713 # Simulation deselection. 714 else: 715 # Deselect. 716 cdp.select_sim[sim_index] = False
717 718
719 - def duplicate_data(self, pipe_from=None, pipe_to=None, model_info=None, global_stats=False, verbose=True):
720 """Duplicate the data specific to a single model-free model. 721 722 @keyword pipe_from: The data pipe to copy the data from. 723 @type pipe_from: str 724 @keyword pipe_to: The data pipe to copy the data to. 725 @type pipe_to: str 726 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 727 @type model_info: int 728 @keyword global_stats: The global statistics flag. 729 @type global_stats: bool 730 @keyword verbose: A flag which if True will cause info about each spin to be printed out as the sequence is generated. 731 @type verbose: bool 732 """ 733 734 # Arg tests. 735 if model_info == None: 736 raise RelaxError("The model_info argument cannot be None.") 737 738 # First create the pipe_to data pipe, if it doesn't exist, but don't switch to it. 739 if not pipes.has_pipe(pipe_to): 740 pipes.create(pipe_to, pipe_type='mf', switch=False) 741 742 # Get the data pipes. 743 dp_from = pipes.get_pipe(pipe_from) 744 dp_to = pipes.get_pipe(pipe_to) 745 746 # Duplicate all non-sequence specific data. 747 for data_name in dir(dp_from): 748 # Skip the container objects. 749 if data_name in ['diff_tensor', 'mol', 'interatomic', 'structure', 'exp_info']: 750 continue 751 752 # Skip special objects. 753 if search('^_', data_name) or data_name in dp_from.__class__.__dict__: 754 continue 755 756 # Get the original object. 757 data_from = getattr(dp_from, data_name) 758 759 # The data already exists. 760 if hasattr(dp_to, data_name): 761 # Get the object in the target pipe. 762 data_to = getattr(dp_to, data_name) 763 764 # The data must match! 765 if data_from != data_to: 766 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 767 768 # Skip the data. 769 continue 770 771 # Duplicate the data. 772 setattr(dp_to, data_name, deepcopy(data_from)) 773 774 # Diffusion tensor comparison. 775 if hasattr(dp_from, 'diff_tensor'): 776 # Duplicate the tensor if it doesn't exist. 777 if not hasattr(dp_to, 'diff_tensor'): 778 setattr(dp_to, 'diff_tensor', deepcopy(dp_from.diff_tensor)) 779 780 # Otherwise compare the objects inside the container. 781 else: 782 # Loop over the modifiable objects. 783 for data_name in dp_from.diff_tensor._mod_attr: 784 # Get the original object. 785 data_from = None 786 if hasattr(dp_from.diff_tensor, data_name): 787 data_from = getattr(dp_from.diff_tensor, data_name) 788 789 # Get the target object. 790 if data_from and not hasattr(dp_to.diff_tensor, data_name): 791 raise RelaxError("The diffusion tensor object " + repr(data_name) + " of the " + repr(pipe_from) + " data pipe is not located in the " + repr(pipe_to) + " data pipe.") 792 elif data_from: 793 data_to = getattr(dp_to.diff_tensor, data_name) 794 else: 795 continue 796 797 # The data must match! 798 if data_from != data_to: 799 raise RelaxError("The object " + repr(data_name) + " is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 800 801 # Structure comparison. 802 if hasattr(dp_from, 'structure'): 803 # Duplicate the tensor if it doesn't exist. 804 if not hasattr(dp_to, 'structure'): 805 setattr(dp_to, 'structure', deepcopy(dp_from.structure)) 806 807 # Otherwise compare the objects inside the container. 808 else: 809 # Modifiable object checks. 810 compare_objects(dp_from.structure, dp_to.structure, pipe_from, pipe_to) 811 812 # Tests for the model and molecule containers. 813 if len(dp_from.structure.structural_data) != len(dp_from.structure.structural_data): 814 raise RelaxError("The number of structural models is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 815 816 # Loop over the models. 817 for i in range(len(dp_from.structure.structural_data)): 818 # Alias. 819 model_from = dp_from.structure.structural_data[i] 820 model_to = dp_to.structure.structural_data[i] 821 822 # Model numbers. 823 if model_from.num != model_to.num: 824 raise RelaxError("The structure models are not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 825 826 # Molecule number. 827 if len(model_from.mol) != len(model_to.mol): 828 raise RelaxError("The number of molecules is not consistent between the pipes " + repr(pipe_from) + " and " + repr(pipe_to) + ".") 829 830 # Loop over the models. 831 for mol_index in range(len(model_from.mol)): 832 # Modifiable object checks. 833 compare_objects(model_from.mol[mol_index], model_to.mol[mol_index], pipe_from, pipe_to) 834 835 # No sequence data, so skip the rest. 836 if dp_from.mol.is_empty(): 837 return 838 839 # Duplicate the sequence data if it doesn't exist. 840 if dp_to.mol.is_empty(): 841 sequence.copy(pipe_from=pipe_from, pipe_to=pipe_to, preserve_select=True, verbose=verbose) 842 843 # Duplicate the interatomic data if it doesn't exist. 844 if dp_to.interatomic.is_empty(): 845 interatomic.copy(pipe_from=pipe_from, pipe_to=pipe_to, verbose=verbose) 846 847 # Determine the model type of the original data pipe. 848 pipes.switch(pipe_from) 849 model_type = determine_model_type() 850 851 # Sequence specific data. 852 spin, spin_id = return_spin_from_index(global_index=model_info, pipe=pipe_from, return_spin_id=True) 853 if model_type == 'mf' or (model_type == 'local_tm' and not global_stats): 854 # Duplicate the spin specific data. 855 for name in dir(spin): 856 # Skip special objects. 857 if search('^__', name): 858 continue 859 860 # Get the object. 861 obj = getattr(spin, name) 862 863 # Skip methods. 864 if isinstance(obj, MethodType): 865 continue 866 867 # Duplicate the object. 868 new_obj = deepcopy(getattr(spin, name)) 869 setattr(dp_to.mol[spin._mol_index].res[spin._res_index].spin[spin._spin_index], name, new_obj) 870 871 # Duplicate the relaxation active spins which have not been copied yet. 872 interatoms = interatomic.return_interatom_list(spin_hash=spin._hash) 873 for interatom in interatoms: 874 # No relaxation mechanism. 875 if not interatom.dipole_pair: 876 continue 877 878 # The interacting spin. 879 if spin_id != interatom.spin_id1: 880 spin_id2 = interatom.spin_id1 881 else: 882 spin_id2 = interatom.spin_id2 883 spin2 = return_spin(spin_id=spin_id2) 884 885 # Duplicate the spin specific data. 886 mol_index, res_index, spin_index = return_spin_indices(spin_id=spin_id2) 887 dp_to.mol[mol_index].res[res_index].spin[spin_index] = deepcopy(spin2) 888 889 # Other data types. 890 else: 891 # Duplicate all the spin specific data. 892 dp_to.mol = deepcopy(dp_from.mol)
893 894
895 - def eliminate(self, name, value, args, sim=None, model_info=None):
896 """Model-free model elimination, parameter by parameter. 897 898 @param name: The parameter name. 899 @type name: str 900 @param value: The parameter value. 901 @type value: float 902 @param args: The c1 and c2 elimination constant overrides. 903 @type args: None or tuple of float 904 @keyword sim: The Monte Carlo simulation index. 905 @type sim: int 906 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 907 @type model_info: int 908 @return: True if the model is to be eliminated, False otherwise. 909 @rtype: bool 910 """ 911 912 # Default values. 913 c1 = 50.0 * 1e-9 914 c2 = 1.5 915 916 # Depack the arguments. 917 if args != None: 918 c1, c2 = args 919 920 # Determine the model type. 921 model_type = determine_model_type() 922 923 # Can't handle this one yet! 924 if model_type != 'mf' and model_type != 'local_tm': 925 raise RelaxError("Elimination of the global model is not yet supported.") 926 927 # Get the spin and it's id string. 928 spin, spin_id = return_spin_from_index(global_index=model_info, return_spin_id=True) 929 930 # Get the tm value. 931 if model_type == 'local_tm': 932 tm = spin.local_tm 933 else: 934 tm = cdp.diff_tensor.tm 935 936 # No tm value set, so skip the tests (no elimination). 937 if tm == None: 938 return False 939 940 # Local tm. 941 if name == 'local_tm' and value >= c1: 942 if sim == None: 943 print("Data pipe '%s': The local tm parameter of %.5g is greater than %.5g, eliminating spin system '%s'." % (pipes.cdp_name(), value, c1, spin_id)) 944 else: 945 print("Data pipe '%s': The local tm parameter of %.5g is greater than %.5g, eliminating simulation %i of spin system '%s'." % (pipes.cdp_name(), value, c1, sim, spin_id)) 946 return True 947 948 # Internal correlation times. 949 if match('t[efs]', name) and value >= c2 * tm: 950 if sim == None: 951 print("Data pipe '%s': The %s value of %.5g is greater than %.5g, eliminating spin system '%s'." % (pipes.cdp_name(), name, value, c2*tm, spin_id)) 952 else: 953 print("Data pipe '%s': The %s value of %.5g is greater than %.5g, eliminating simulation %i of spin system '%s'." % (pipes.cdp_name(), name, value, c2*tm, sim, spin_id)) 954 return True 955 956 # Accept model. 957 return False
958 959
960 - def get_param_names(self, model_info=None):
961 """Return a vector of parameter names. 962 963 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 964 @type model_info: int 965 @return: The vector of parameter names. 966 @rtype: list of str 967 """ 968 969 # Determine the model type. 970 model_type = determine_model_type() 971 972 # Get the spin ids. 973 if model_type == 'mf' or model_type == 'local_tm': 974 # Get the spin and it's id string. 975 spin, spin_id = return_spin_from_index(global_index=model_info, return_spin_id=True) 976 else: 977 spin_id = None 978 979 # Assemble and return the parameter names. 980 return assemble_param_names(model_type, spin_id=spin_id)
981 982
983 - def get_param_values(self, model_info=None, sim_index=None):
984 """Return a vector of parameter values. 985 986 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 987 @type model_info: int 988 @keyword sim_index: The Monte Carlo simulation index. 989 @type sim_index: int 990 @return: The vector of parameter values. 991 @rtype: list of str 992 """ 993 994 # Test if the model-free models have been set up. 995 for spin in spin_loop(): 996 # Skip deselected spins. 997 if not spin.select: 998 continue 999 1000 # Not setup. 1001 if not spin.model: 1002 raise RelaxNoModelError 1003 1004 # Determine the model type. 1005 model_type = determine_model_type() 1006 1007 # Set the spin container (to None if the model is global). 1008 if model_type == 'mf' or model_type == 'local_tm': 1009 spin = return_spin_from_index(global_index=model_info) 1010 else: 1011 spin = None 1012 1013 # Assemble the parameter values and return them. 1014 return assemble_param_vector(spin=spin, sim_index=sim_index, model_type=model_type)
1015 1016
1017 - def grid_search(self, lower=None, upper=None, inc=None, scaling_matrix=None, constraints=True, verbosity=1, sim_index=None):
1018 """The model-free grid search function. 1019 1020 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model. 1021 @type lower: list of lists of numbers 1022 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model. 1023 @type upper: list of lists of numbers 1024 @keyword inc: The per-model increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. 1025 @type inc: list of lists of int 1026 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 1027 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 1028 @keyword constraints: If True, constraints are applied during the grid search (eliminating parts of the grid). If False, no constraints are used. 1029 @type constraints: bool 1030 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 1031 @type verbosity: int 1032 @keyword sim_index: The index of the simulation to apply the grid search to. If None, the normal model is optimised. 1033 @type sim_index: int 1034 """ 1035 1036 # Minimisation. 1037 self.minimise(min_algor='grid', lower=lower, upper=upper, inc=inc, scaling_matrix=scaling_matrix, constraints=constraints, verbosity=verbosity, sim_index=sim_index)
1038 1039
1040 - def map_bounds(self, param, spin_id=None):
1041 """Create bounds for the OpenDX mapping function. 1042 1043 @param param: The name of the parameter to return the lower and upper bounds of. 1044 @type param: str 1045 @param spin_id: The spin identification string. 1046 @type spin_id: str 1047 @return: The upper and lower bounds of the parameter. 1048 @rtype: list of float 1049 """ 1050 1051 # Diffusion tensor bounds. 1052 if self._PARAMS.scope(param) == 'global': 1053 return diffusion_tensor.map_bounds(param) 1054 1055 # Get the spin. 1056 spin = return_spin(spin_id=spin_id) 1057 1058 # {s2, s2f, s2s}. 1059 if search('^s2', param): 1060 return [0.0, 1.0] 1061 1062 # {local tm, te, tf, ts}. 1063 elif search('^t', param) or param == 'local_tm': 1064 return [0.0, 1e-8] 1065 1066 # Rex. 1067 elif param == 'rex': 1068 return [0.0, 30.0 / (2.0 * pi * cdp.spectrometer_frq[cdp.ri_ids[0]])**2] 1069 1070 # Interatomic distances. 1071 elif param == 'r': 1072 return [1.0 * 1e-10, 1.1 * 1e-10] 1073 1074 # CSA. 1075 elif param == 'csa': 1076 return [-100 * 1e-6, -300 * 1e-6]
1077 1078
1079 - def minimise(self, min_algor=None, min_options=None, func_tol=None, grad_tol=None, max_iterations=None, constraints=False, scaling_matrix=None, verbosity=0, sim_index=None, lower=None, upper=None, inc=None):
1080 """Model-free minimisation function. 1081 1082 Three categories of models exist for which the approach to minimisation is different. These 1083 are: 1084 1085 Single spin optimisations: The 'mf' and 'local_tm' model types which are the 1086 model-free parameters for single spins, optionally with a local tm parameter. These 1087 models have no global parameters. 1088 1089 Diffusion tensor optimisations: The 'diff' diffusion tensor model type. No spin 1090 specific parameters exist. 1091 1092 Optimisation of everything: The 'all' model type consisting of all model-free and all 1093 diffusion tensor parameters. 1094 1095 1096 @keyword min_algor: The minimisation algorithm to use. 1097 @type min_algor: str 1098 @keyword min_options: An array of options to be used by the minimisation algorithm. 1099 @type min_options: array of str 1100 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1101 @type func_tol: None or float 1102 @keyword grad_tol: The gradient tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 1103 @type grad_tol: None or float 1104 @keyword max_iterations: The maximum number of iterations for the algorithm. 1105 @type max_iterations: int 1106 @keyword constraints: If True, constraints are used during optimisation. 1107 @type constraints: bool 1108 @keyword scaling_matrix: The per-model list of diagonal and square scaling matrices. 1109 @type scaling_matrix: list of numpy rank-2, float64 array or list of None 1110 @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. 1111 @type verbosity: int 1112 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 1113 @type sim_index: None or int 1114 @keyword lower: The per-model lower bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 1115 @type lower: list of lists of numbers 1116 @keyword upper: The per-model upper bounds of the grid search which must be equal to the number of parameters in the model. This optional argument is only used when doing a grid search. 1117 @type upper: list of lists of numbers 1118 @keyword inc: The per-model increments for each dimension of the space for the grid search. The number of elements in the array must equal to the number of parameters in the model. This argument is only used when doing a grid search. 1119 @type inc: list of lists of int 1120 """ 1121 1122 # Test if sequence data is loaded. 1123 if not exists_mol_res_spin_data(): 1124 raise RelaxNoSequenceError 1125 1126 # Test if the model-free model has been setup, and that the nuclear isotope types have been set. 1127 for spin in spin_loop(): 1128 # Skip deselected spins. 1129 if not spin.select: 1130 continue 1131 1132 # Not setup. 1133 if not spin.model: 1134 raise RelaxNoModelError 1135 1136 # Test if the nuclear isotope type has been set. 1137 if not hasattr(spin, 'isotope'): 1138 raise RelaxSpinTypeError 1139 1140 # Containers for the model-free data and optimisation parameters. 1141 data_store = Data_container() 1142 opt_params = Data_container() 1143 1144 # Add the imported parameters. 1145 data_store.h_bar = h_bar 1146 data_store.mu0 = mu0 1147 opt_params.min_algor = min_algor 1148 opt_params.min_options = min_options 1149 opt_params.func_tol = func_tol 1150 opt_params.grad_tol = grad_tol 1151 opt_params.max_iterations = max_iterations 1152 1153 # Add the keyword args. 1154 opt_params.verbosity = verbosity 1155 1156 # Determine the model type. 1157 data_store.model_type = determine_model_type() 1158 if not data_store.model_type: 1159 return 1160 1161 # Model type for the back-calculate function. 1162 if min_algor == 'back_calc' and data_store.model_type != 'local_tm': 1163 data_store.model_type = 'mf' 1164 1165 # Test if diffusion tensor data exists. 1166 if data_store.model_type != 'local_tm' and not diffusion_tensor.diff_data_exists(): 1167 raise RelaxNoTensorError('diffusion') 1168 1169 # Tests for the PDB file and unit vectors. 1170 if data_store.model_type != 'local_tm' and cdp.diff_tensor.type != 'sphere': 1171 # Test if the structure file has been loaded. 1172 if not hasattr(cdp, 'structure'): 1173 raise RelaxNoPdbError 1174 1175 # Test if unit vectors exist. 1176 for spin, spin_id in spin_loop(return_id=True): 1177 # Skip deselected spins. 1178 if not spin.select: 1179 continue 1180 1181 # Get the interatomic data container. 1182 interatoms = return_interatom_list(spin_hash=spin._hash) 1183 1184 # Unit vectors. 1185 for i in range(len(interatoms)): 1186 # No relaxation mechanism. 1187 if not interatoms[i].dipole_pair: 1188 continue 1189 1190 # Check for the vectors. 1191 if not hasattr(interatoms[i], 'vector'): 1192 raise RelaxNoVectorsError 1193 1194 # Test if the model-free parameter values are set for minimising diffusion tensor parameters by themselves. 1195 if data_store.model_type == 'diff': 1196 # Loop over the sequence. 1197 for spin in spin_loop(): 1198 unset_param = are_mf_params_set(spin) 1199 if unset_param != None: 1200 raise RelaxNoValueError(unset_param) 1201 1202 # Print out. 1203 if verbosity >= 1: 1204 if data_store.model_type == 'mf': 1205 print("Only the model-free parameters for single spins will be used.") 1206 elif data_store.model_type == 'local_mf': 1207 print("Only a local tm value together with the model-free parameters for single spins will be used.") 1208 elif data_store.model_type == 'diff': 1209 print("Only diffusion tensor parameters will be used.") 1210 elif data_store.model_type == 'all': 1211 print("The diffusion tensor parameters together with the model-free parameters for all spins will be used.") 1212 1213 # Test if the CSA and interatomic distances have been set. 1214 for spin, spin_id in spin_loop(return_id=True): 1215 # Skip deselected spins. 1216 if not spin.select: 1217 continue 1218 1219 # CSA value. 1220 if not hasattr(spin, 'csa') or spin.csa == None: 1221 raise RelaxNoValueError("CSA") 1222 1223 # Get the interatomic data container. 1224 interatoms = return_interatom_list(spin_hash=spin._hash) 1225 1226 # Interatomic distances. 1227 count = 0 1228 for i in range(len(interatoms)): 1229 # No relaxation mechanism. 1230 if not interatoms[i].dipole_pair: 1231 continue 1232 1233 # Check for the distances. 1234 if not hasattr(interatoms[i], 'r') or interatoms[i].r == None: 1235 raise RelaxNoValueError("interatomic distance", spin_id=interatoms[i].spin_id1, spin_id2=interatoms[i].spin_id2) 1236 1237 # Count the number of interactions. 1238 count += 1 1239 1240 # Too many interactions. 1241 if count > 1: 1242 raise RelaxError("The spin '%s' has %s dipolar relaxation interactions defined, but only a maximum of one is currently supported." % (spin_id, count)) 1243 1244 # Number of spins, minimisation instances, and data sets for each model type. 1245 if data_store.model_type == 'mf' or data_store.model_type == 'local_tm': 1246 num_data_sets = 1 1247 data_store.num_spins = 1 1248 elif data_store.model_type == 'diff' or data_store.model_type == 'all': 1249 num_data_sets = count_spins(skip_desel=False) 1250 data_store.num_spins = count_spins() 1251 1252 # Number of spins, minimisation instances, and data sets for the back-calculate function. 1253 if min_algor == 'back_calc': 1254 num_data_sets = 0 1255 data_store.num_spins = 1 1256 1257 # Get the Processor box singleton (it contains the Processor instance) and alias the Processor. 1258 processor_box = Processor_box() 1259 processor = processor_box.processor 1260 1261 # Loop over the models. 1262 for index in self.model_loop(): 1263 # Get the spin container if required. 1264 if data_store.model_type == 'diff' or data_store.model_type == 'all': 1265 spin_index = None 1266 spin, data_store.spin_id = None, None 1267 elif min_algor == 'back_calc': 1268 spin_index = opt_params.min_options[0] 1269 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1270 else: 1271 spin_index = index 1272 spin, data_store.spin_id = return_spin_from_index(global_index=spin_index, return_spin_id=True) 1273 1274 # Individual spin stuff. 1275 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm') and not min_algor == 'back_calc': 1276 # Skip deselected spins. 1277 if not spin.select: 1278 continue 1279 1280 # Skip spins missing relaxation data or errors. 1281 if not hasattr(spin, 'ri_data') or not hasattr(spin, 'ri_data_err'): 1282 continue 1283 1284 # Skip spins missing the dipolar interaction. 1285 if spin and (data_store.model_type == 'mf' or data_store.model_type == 'local_tm'): 1286 interatoms = return_interatom_list(spin_hash=spin._hash) 1287 if not len(interatoms): 1288 continue 1289 1290 # Parameter vector and diagonal scaling. 1291 if min_algor == 'back_calc': 1292 # Create the initial parameter vector. 1293 opt_params.param_vector = assemble_param_vector(spin=spin, model_type=data_store.model_type) 1294 1295 # Diagonal scaling. 1296 data_store.scaling_matrix = None 1297 1298 else: 1299 # Create the initial parameter vector. 1300 opt_params.param_vector = assemble_param_vector(spin=spin, sim_index=sim_index) 1301 1302 # The number of parameters. 1303 num_params = len(opt_params.param_vector) 1304 1305 # Diagonal scaling. 1306 data_store.scaling_matrix = scaling_matrix[index] 1307 if data_store.scaling_matrix is not None: 1308 opt_params.param_vector = dot(inv(data_store.scaling_matrix), opt_params.param_vector) 1309 1310 # Store the grid search options. 1311 opt_params.lower, opt_params.upper, opt_params.inc = None, None, None 1312 if lower != None: 1313 opt_params.lower = lower[index] 1314 if upper != None: 1315 opt_params.upper = upper[index] 1316 if inc != None: 1317 opt_params.inc = inc[index] 1318 1319 # Scaling of values for the set function. 1320 if match('^[Ss]et', min_algor): 1321 opt_params.min_options = dot(inv(data_store.scaling_matrix), opt_params.min_options) 1322 1323 # Linear constraints. 1324 if constraints: 1325 opt_params.A, opt_params.b = linear_constraints(num_params, model_type=data_store.model_type, spin=spin, scaling_matrix=data_store.scaling_matrix) 1326 else: 1327 opt_params.A, opt_params.b = None, None 1328 1329 # Get the data for minimisation. 1330 minimise_data_setup(data_store, min_algor, num_data_sets, opt_params.min_options, spin=spin, sim_index=sim_index) 1331 1332 # Setup the minimisation algorithm when constraints are present. 1333 if constraints and not match('^[Gg]rid', min_algor): 1334 algor = opt_params.min_options[0] 1335 else: 1336 algor = min_algor 1337 1338 # Initialise the function to minimise (for back-calculation and LM minimisation). 1339 if min_algor == 'back_calc' or match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 1340 mf = Mf(init_params=opt_params.param_vector, model_type=data_store.model_type, diff_type=data_store.diff_type, diff_params=data_store.diff_params, scaling_matrix=data_store.scaling_matrix, num_spins=data_store.num_spins, equations=data_store.equations, param_types=data_store.param_types, param_values=data_store.param_values, relax_data=data_store.ri_data, errors=data_store.ri_data_err, bond_length=data_store.r, csa=data_store.csa, num_frq=data_store.num_frq, frq=data_store.frq, num_ri=data_store.num_ri, remap_table=data_store.remap_table, noe_r1_table=data_store.noe_r1_table, ri_labels=data_store.ri_types, gx=data_store.gx, gh=data_store.gh, h_bar=data_store.h_bar, mu0=data_store.mu0, num_params=data_store.num_params, vectors=data_store.xh_unit_vectors) 1341 1342 # Levenberg-Marquardt minimisation. 1343 if match('[Ll][Mm]$', algor) or match('[Ll]evenburg-[Mm]arquardt$', algor): 1344 # Total number of ri. 1345 number_ri = 0 1346 for k in range(len(ri_data_err)): 1347 number_ri = number_ri + len(ri_data_err[k]) 1348 1349 # Reconstruct the error data structure. 1350 lm_error = zeros(number_ri, float64) 1351 index = 0 1352 for k in range(len(ri_data_err)): 1353 lm_error[index:index+len(ri_data_err[k])] = ri_data_err[k] 1354 index = index + len(ri_data_err[k]) 1355 1356 opt_params.min_options = opt_params.min_options + (mf.lm_dri, lm_error) 1357 1358 # Back-calculation. 1359 if min_algor == 'back_calc': 1360 return mf.calc_ri() 1361 1362 # Parallelised grid search for the diffusion parameter space. 1363 if match('^[Gg]rid', min_algor) and data_store.model_type == 'diff': 1364 # Print out. 1365 print("Parallelised diffusion tensor grid search.") 1366 1367 # Loop over each grid subdivision. 1368 for subdivision in grid_split(divisions=processor.processor_size(), lower=opt_params.lower, upper=opt_params.upper, inc=opt_params.inc, verbosity=verbosity): 1369 # Set the points. 1370 opt_params.subdivision = subdivision 1371 1372 # Grid search initialisation. 1373 command = MF_grid_command() 1374 1375 # Pass in the data and optimisation parameters. 1376 command.store_data(deepcopy(data_store), deepcopy(opt_params)) 1377 1378 # Set up the model-free memo and add it to the processor queue. 1379 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling_matrix=data_store.scaling_matrix) 1380 processor.add_to_queue(command, memo) 1381 1382 # Execute the queued elements. 1383 processor.run_queue() 1384 1385 # Exit this method. 1386 return 1387 1388 # Normal grid search (command initialisation). 1389 if search('^[Gg]rid', min_algor): 1390 command = MF_grid_command() 1391 1392 # Minimisation of all other model types (command initialisation). 1393 else: 1394 command = MF_minimise_command() 1395 1396 # Pass in the data and optimisation parameters. 1397 command.store_data(deepcopy(data_store), deepcopy(opt_params)) 1398 1399 # Set up the model-free memo and add it to the processor queue. 1400 memo = MF_memo(model_free=self, model_type=data_store.model_type, spin=spin, sim_index=sim_index, scaling_matrix=data_store.scaling_matrix) 1401 processor.add_to_queue(command, memo) 1402 1403 # Execute the queued elements. 1404 processor.run_queue()
1405 1406
1407 - def model_desc(self, model_info=None):
1408 """Return a description of the model. 1409 1410 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1411 @type model_info: int 1412 @return: The model description. 1413 @rtype: str 1414 """ 1415 1416 # Determine the model type. 1417 model_type = determine_model_type() 1418 1419 # Global models. 1420 if model_type == 'all': 1421 return "Global model - all diffusion tensor parameters and spin specific model-free parameters." 1422 elif model_type == 'diff': 1423 return "Diffusion tensor model." 1424 1425 # Spin specific model. 1426 else: 1427 # Get the spin container. 1428 spin, spin_id = return_spin_from_index(model_info, return_spin_id=True) 1429 1430 # Return the description. 1431 return "Model-free model of spin '%s'." % spin_id
1432 1433
1434 - def model_loop(self):
1435 """Generator method for looping over the models (global or local). 1436 1437 If the model type is 'all' or 'diff', then this yields the single value of zero. Otherwise 1438 the global spin index is yielded. 1439 1440 1441 @return: The model index. This index is zero for the global models or equal to the global spin 1442 index (which covers the molecule, residue, and spin indices). 1443 @rtype: int 1444 """ 1445 1446 # Determine the model type. 1447 model_type = determine_model_type() 1448 1449 # Global model. 1450 if model_type == 'all' or model_type == 'diff': 1451 yield 0 1452 1453 # Spin specific models. 1454 else: 1455 # Loop over the spins. 1456 global_index = -1 1457 for spin in spin_loop(): 1458 # Increment the global spin index. 1459 global_index = global_index + 1 1460 1461 # Yield the spin index. 1462 yield global_index
1463 1464
1465 - def model_statistics(self, model_info=None, spin_id=None, global_stats=None):
1466 """Return the k, n, and chi2 model statistics. 1467 1468 k - number of parameters. 1469 n - number of data points. 1470 chi2 - the chi-squared value. 1471 1472 1473 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1474 @type model_info: int 1475 @keyword spin_id: The spin identification string. Either this or the instance keyword argument must be supplied. 1476 @type spin_id: None or str 1477 @keyword global_stats: A parameter which determines if global or local statistics are returned. If None, then the appropriateness of global or local statistics is automatically determined. 1478 @type global_stats: None or bool 1479 @return: The optimisation statistics, in tuple format, of the number of parameters (k), the number of data points (n), and the chi-squared value (chi2). 1480 @rtype: tuple of (int, int, float) 1481 """ 1482 1483 # Bad argument combination. 1484 if model_info == None and spin_id == None: 1485 raise RelaxError("Either the model_info or spin_id argument must be supplied.") 1486 elif model_info != None and spin_id != None: 1487 raise RelaxError("The model_info arg " + repr(model_info) + " and spin_id arg " + repr(spin_id) + " clash. Only one should be supplied.") 1488 1489 # Determine the model type. 1490 model_type = determine_model_type() 1491 1492 # Determine if local or global statistics will be returned. 1493 if global_stats == None: 1494 if model_type in ['mf', 'local_tm']: 1495 global_stats = False 1496 else: 1497 global_stats = True 1498 1499 # Statistics for a single residue. 1500 if not global_stats: 1501 # Get the SpinContainer. 1502 if spin_id: 1503 spin = return_spin(spin_id) 1504 else: 1505 spin = return_spin_from_index(model_info) 1506 1507 # Skip deselected residues. 1508 if not spin.select: 1509 return None, None, None 1510 1511 # Missing data sets. 1512 if not hasattr(spin, 'ri_data'): 1513 return None, None, None 1514 1515 # Count the number of parameters. 1516 param_vector = assemble_param_vector(spin=spin) 1517 k = len(param_vector) 1518 1519 # Count the number of data points. 1520 n = len(spin.ri_data) 1521 1522 # The chi2 value. 1523 chi2 = spin.chi2 1524 1525 # Global stats. 1526 elif global_stats: 1527 # Count the number of parameters. 1528 param_vector = assemble_param_vector() 1529 k = len(param_vector) 1530 1531 # Count the number of data points. 1532 n = 0 1533 chi2 = 0 1534 for spin in spin_loop(): 1535 # Skip deselected residues. 1536 if not spin.select: 1537 continue 1538 1539 # Skip residues with no relaxation data. 1540 if not hasattr(spin, 'ri_data') or not len(spin.ri_data): 1541 continue 1542 1543 n = n + len(spin.ri_data) 1544 1545 # Local tm models. 1546 if model_type == 'local_tm': 1547 chi2 = chi2 + spin.chi2 1548 1549 # The chi2 value. 1550 if model_type != 'local_tm': 1551 if not hasattr(cdp, 'chi2'): 1552 raise RelaxError("Global statistics are not available, most likely because the global model has not been optimised.") 1553 chi2 = cdp.chi2 1554 1555 # Return the data. 1556 return k, n, chi2
1557 1558
1559 - def model_type(self):
1560 """Return the type of the model, either being 'local' or 'global'. 1561 1562 @return: The model type, one of 'local' or 'global'. 1563 @rtype: str 1564 """ 1565 1566 # Determine the model type. 1567 model_type = determine_model_type() 1568 1569 # Global models. 1570 if model_type in ['all', 'diff']: 1571 return 'global' 1572 1573 # Local models. 1574 else: 1575 return 'local'
1576 1577
1578 - def num_instances(self):
1579 """Function for returning the number of instances. 1580 1581 @return: The number of instances used for optimisation. Either the number of spins if 1582 the local optimisations are setup ('mf' and 'local_tm'), or 1 for the global 1583 models. 1584 @rtype: int 1585 """ 1586 1587 # Test if sequence data exists. 1588 if not exists_mol_res_spin_data(): 1589 return 0 1590 1591 # Determine the model type. 1592 model_type = determine_model_type() 1593 1594 # Sequence specific data. 1595 if model_type == 'mf' or model_type == 'local_tm': 1596 return count_spins() 1597 1598 # Other data. 1599 elif model_type == 'diff' or model_type == 'all': 1600 return 1 1601 1602 # Should not be here. 1603 else: 1604 raise RelaxFault
1605 1606
1607 - def overfit_deselect(self, data_check=True, verbose=True):
1608 """Deselect spins which have insufficient data to support minimisation. 1609 1610 @keyword data_check: A flag to signal if the presence of base data is to be checked for. 1611 @type data_check: bool 1612 @keyword verbose: A flag which if True will allow printouts. 1613 @type verbose: bool 1614 """ 1615 1616 # Print out. 1617 if verbose: 1618 print("\nOver-fit spin deselection:") 1619 1620 # Test if sequence data exists. 1621 if not exists_mol_res_spin_data(): 1622 raise RelaxNoSequenceError 1623 1624 # Is structural data required? 1625 need_vect = False 1626 if hasattr(cdp, 'diff_tensor') and (cdp.diff_tensor.type == 'spheroid' or cdp.diff_tensor.type == 'ellipsoid'): 1627 need_vect = True 1628 1629 # Loop over the sequence. 1630 deselect_flag = False 1631 spin_count = 0 1632 for spin, spin_id in spin_loop(return_id=True): 1633 # Skip deselected spins. 1634 if not spin.select: 1635 continue 1636 1637 # The interatomic data. 1638 interatoms = interatomic.return_interatom_list(spin_hash=spin._hash) 1639 1640 # Loop over the interatomic data. 1641 dipole_relax = False 1642 for i in range(len(interatoms)): 1643 # No dipolar relaxation mechanism. 1644 if not interatoms[i].dipole_pair: 1645 continue 1646 1647 # The surrounding spins. 1648 if spin_id != interatoms[i].spin_id1: 1649 spin_id2 = interatoms[i].spin_id1 1650 else: 1651 spin_id2 = interatoms[i].spin_id2 1652 spin2 = return_spin(spin_id2) 1653 1654 # Dipolar relaxation flag. 1655 dipole_relax = True 1656 1657 # No relaxation mechanism. 1658 if not dipole_relax or not hasattr(spin, 'csa') or spin.csa == None: 1659 warn(RelaxDeselectWarning(spin_id, 'an absence of relaxation mechanisms')) 1660 spin.select = False 1661 deselect_flag = True 1662 continue 1663 1664 # Data checks. 1665 if data_check: 1666 # The number of relaxation data points (and for infinite data). 1667 data_points = 0 1668 inf_data = False 1669 if hasattr(cdp, 'ri_ids') and hasattr(spin, 'ri_data'): 1670 for id in cdp.ri_ids: 1671 if id in spin.ri_data and spin.ri_data[id] != None: 1672 data_points += 1 1673 1674 # Infinite data! 1675 if isInf(spin.ri_data[id]): 1676 inf_data = True 1677 1678 # Infinite data. 1679 if inf_data: 1680 warn(RelaxDeselectWarning(spin_id, 'infinite relaxation data')) 1681 spin.select = False 1682 deselect_flag = True 1683 continue 1684 1685 # Relaxation data must exist! 1686 if not hasattr(spin, 'ri_data'): 1687 warn(RelaxDeselectWarning(spin_id, 'missing relaxation data')) 1688 spin.select = False 1689 deselect_flag = True 1690 continue 1691 1692 # Require 3 or more relaxation data points. 1693 elif data_points < 3: 1694 warn(RelaxDeselectWarning(spin_id, 'insufficient relaxation data, 3 or more data points are required')) 1695 spin.select = False 1696 deselect_flag = True 1697 continue 1698 1699 # Require at least as many data points as params to prevent over-fitting. 1700 elif hasattr(spin, 'params') and spin.params and len(spin.params) > data_points: 1701 warn(RelaxDeselectWarning(spin_id, 'over-fitting - more parameters than data points')) 1702 spin.select = False 1703 deselect_flag = True 1704 continue 1705 1706 # Test for structural data if required. 1707 for i in range(len(interatoms)): 1708 # No dipolar relaxation mechanism. 1709 if not interatoms[i].dipole_pair: 1710 continue 1711 1712 # Check the unit vectors. 1713 if need_vect: 1714 if not hasattr(interatoms[i], 'vector') or interatoms[i].vector is None: 1715 warn(RelaxDeselectWarning(spin_id, 'missing structural data')) 1716 spin.select = False 1717 deselect_flag = True 1718 continue 1719 1720 # Increment the spin number. 1721 spin_count += 1 1722 1723 # No spins selected, so fail hard to prevent the user from going any further. 1724 if spin_count == 0: 1725 warn(RelaxWarning("No spins are selected therefore the optimisation or calculation cannot proceed.")) 1726 1727 # Final printout. 1728 if verbose and not deselect_flag: 1729 print("No spins have been deselected.")
1730 1731
1732 - def print_model_title(self, prefix=None, model_info=None):
1733 """Print out the model title. 1734 1735 @keyword prefix: The starting text of the title. This should be printed out first, followed by the model information text. 1736 @type prefix: str 1737 @keyword model_info: The model information from model_loop(). 1738 @type model_info: unknown 1739 """ 1740 1741 # Determine the model type. 1742 model_type = determine_model_type() 1743 1744 # Local models. 1745 if model_type == 'mf' or model_type == 'local_tm': 1746 spin, spin_id = return_spin_from_index(global_index=model_info, return_spin_id=True) 1747 text = "%sSpin '%s'" % (prefix, spin_id) 1748 1749 # Global models. 1750 else: 1751 text = prefix + "Global model" 1752 1753 # The printout. 1754 subsection(file=sys.stdout, text=text, prespace=2)
1755 1756
1757 - def set_error(self, index, error, model_info=None):
1758 """Set the parameter errors. 1759 1760 @param index: The index of the parameter to set the errors for. 1761 @type index: int 1762 @param error: The error value. 1763 @type error: float 1764 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1765 @type model_info: int 1766 """ 1767 1768 # Parameter increment counter. 1769 inc = 0 1770 1771 # Determine the model type. 1772 model_type = determine_model_type() 1773 1774 # Get the parameter object names. 1775 param_names = self.data_names(set='params', scope='spin') 1776 1777 1778 # Diffusion tensor parameter errors. 1779 #################################### 1780 1781 if model_type == 'diff' or model_type == 'all': 1782 # Spherical diffusion. 1783 if cdp.diff_tensor.type == 'sphere': 1784 # Return the parameter array. 1785 if index == 0: 1786 cdp.diff_tensor.set(param='tm', value=error, category='err') 1787 1788 # Increment. 1789 inc = inc + 1 1790 1791 # Spheroidal diffusion. 1792 elif cdp.diff_tensor.type == 'spheroid': 1793 # Return the parameter array. 1794 if index == 0: 1795 cdp.diff_tensor.set(param='tm', value=error, category='err') 1796 elif index == 1: 1797 cdp.diff_tensor.set(param='Da', value=error, category='err') 1798 elif index == 2: 1799 cdp.diff_tensor.set(param='theta', value=error, category='err') 1800 elif index == 3: 1801 cdp.diff_tensor.set(param='phi', value=error, category='err') 1802 1803 # Increment. 1804 inc = inc + 4 1805 1806 # Ellipsoidal diffusion. 1807 elif cdp.diff_tensor.type == 'ellipsoid': 1808 # Return the parameter array. 1809 if index == 0: 1810 cdp.diff_tensor.set(param='tm', value=error, category='err') 1811 elif index == 1: 1812 cdp.diff_tensor.set(param='Da', value=error, category='err') 1813 elif index == 2: 1814 cdp.diff_tensor.set(param='Dr', value=error, category='err') 1815 elif index == 3: 1816 cdp.diff_tensor.set(param='alpha', value=error, category='err') 1817 elif index == 4: 1818 cdp.diff_tensor.set(param='beta', value=error, category='err') 1819 elif index == 5: 1820 cdp.diff_tensor.set(param='gamma', value=error, category='err') 1821 1822 # Increment. 1823 inc = inc + 6 1824 1825 1826 # Model-free parameter errors for the model type 'all'. 1827 ####################################################### 1828 1829 if model_type == 'all': 1830 # Loop over the spins. 1831 for spin in spin_loop(): 1832 # Skip deselected spins. 1833 if not spin.select: 1834 continue 1835 1836 # Loop over the residue specific parameters. 1837 for param in param_names: 1838 # Return the parameter array. 1839 if index == inc: 1840 setattr(spin, param + "_err", error) 1841 1842 # Increment. 1843 inc = inc + 1 1844 1845 1846 # Model-free parameters for the model types 'mf' and 'local_tm'. 1847 ################################################################ 1848 1849 if model_type == 'mf' or model_type == 'local_tm': 1850 # Get the spin container. 1851 spin = return_spin_from_index(model_info) 1852 1853 # Skip deselected residues. 1854 if not spin.select: 1855 return 1856 1857 # Loop over the residue specific parameters. 1858 for param in param_names: 1859 # Return the parameter array. 1860 if index == inc: 1861 setattr(spin, param + "_err", error) 1862 1863 # Increment. 1864 inc = inc + 1
1865 1866
1867 - def set_param_values(self, param=None, value=None, index=None, spin_id=None, error=False, force=True):
1868 """Set the model-free parameter values. 1869 1870 @keyword param: The parameter name list. 1871 @type param: list of str 1872 @keyword value: The parameter value list. 1873 @type value: list 1874 @keyword index: The index for parameters which are of the list-type. This is unused. 1875 @type index: None or int 1876 @keyword spin_id: The spin identification string, only used for spin specific parameters. 1877 @type spin_id: None or str 1878 @keyword error: A flag which if True will allow the parameter errors to be set instead of the values. 1879 @type error: bool 1880 @keyword force: A flag which if True will cause current values to be overwritten. If False, a RelaxError will raised if the parameter value is already set. 1881 @type force: bool 1882 """ 1883 1884 # Checks. 1885 is_str_list(param, 'parameter name') 1886 1887 # Separate out the diffusion tensor parameters from the model-free parameters. 1888 diff_params = [] 1889 diff_vals = [] 1890 mf_params = [] 1891 mf_vals = [] 1892 for i in range(len(param)): 1893 # Diffusion tensor parameter. 1894 if self._PARAMS.scope(param[i]) == 'global': 1895 if error: 1896 diff_params.append(param[i] + '_err') 1897 else: 1898 diff_params.append(param[i]) 1899 diff_vals.append(value[i]) 1900 1901 # Model-free parameter. 1902 else: 1903 mf_params.append(param[i]) 1904 mf_vals.append(value[i]) 1905 1906 # Set the diffusion tensor parameters. 1907 if diff_params: 1908 diffusion_tensor.set(value=diff_vals, param=diff_params) 1909 1910 # Set the model-free parameters. 1911 for i in range(len(mf_params)): 1912 # Check if it is a model-free parameter. 1913 if mf_params[i] not in self.data_names(set='params', scope='spin') and mf_params[i] not in self.data_names(set='generic', scope='spin'): 1914 raise RelaxError("The parameter '%s' is unknown. It should be one of %s or %s" % (mf_params[i], self.data_names(set='params', scope='spin'), self.data_names(set='generic', scope='spin'))) 1915 1916 # The error object name. 1917 if error: 1918 mf_params[i] += '_err' 1919 1920 # Set the parameter. 1921 for spin in spin_loop(spin_id): 1922 setattr(spin, mf_params[i], mf_vals[i])
1923 1924
1925 - def set_selected_sim(self, select_sim, model_info=None):
1926 """Set all simulation selection flags. 1927 1928 @param select_sim: The selection flags. 1929 @type select_sim: bool 1930 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 1931 @type model_info: int 1932 """ 1933 1934 # Determine the model type. 1935 model_type = determine_model_type() 1936 1937 # Global model. 1938 if model_type == 'all' or model_type == 'diff': 1939 cdp.select_sim = select_sim 1940 1941 # Spin specific model. 1942 else: 1943 # Get the spin container. 1944 spin = return_spin_from_index(model_info) 1945 1946 # Skip if deselected. 1947 if not spin.select: 1948 return 1949 1950 # Set the simulation flags. 1951 spin.select_sim = deepcopy(select_sim)
1952 1953
1954 - def set_update(self, param, spin):
1955 """Function to update the other model-free parameters. 1956 1957 @param param: The name of the parameter which has been changed. 1958 @type param: str 1959 @param spin: The SpinContainer object. 1960 @type spin: SpinContainer 1961 """ 1962 1963 # S2f parameter. 1964 if param == 's2f': 1965 # Update S2 if S2s exists. 1966 if hasattr(spin, 's2s') and spin.s2s != None: 1967 spin.s2 = spin.s2f * spin.s2s 1968 1969 1970 # S2s parameter. 1971 if param == 's2s': 1972 # Update S2 if S2f exists. 1973 if hasattr(spin, 's2f') and spin.s2f != None: 1974 spin.s2 = spin.s2f * spin.s2s
1975 1976
1977 - def sim_init_values(self):
1978 """Initialise the Monte Carlo parameter values.""" 1979 1980 # Determine the model type. 1981 model_type = determine_model_type() 1982 1983 # Get the parameter object names. 1984 param_names = self.data_names(set='params', scope='spin') 1985 1986 # Get the minimisation statistic object names. 1987 min_names = self.data_names(set='min', scope='spin') 1988 1989 # List of diffusion tensor parameters. 1990 if model_type == 'diff' or model_type == 'all': 1991 # Spherical diffusion. 1992 if cdp.diff_tensor.type == 'sphere': 1993 diff_params = ['tm'] 1994 1995 # Spheroidal diffusion. 1996 elif cdp.diff_tensor.type == 'spheroid': 1997 diff_params = ['tm', 'Da', 'theta', 'phi'] 1998 1999 # Ellipsoidal diffusion. 2000 elif cdp.diff_tensor.type == 'ellipsoid': 2001 diff_params = ['tm', 'Da', 'Dr', 'alpha', 'beta', 'gamma'] 2002 2003 2004 # Test if Monte Carlo parameter values have already been set. 2005 ############################################################# 2006 2007 # Diffusion tensor parameters and non spin specific minimisation statistics. 2008 if model_type == 'diff' or model_type == 'all': 2009 # Loop over the parameters. 2010 for object_name in diff_params: 2011 # Name for the simulation object. 2012 sim_object_name = object_name + '_sim' 2013 2014 # Test if the simulation object already exists. 2015 if hasattr(cdp.diff_tensor, sim_object_name): 2016 raise RelaxError("Monte Carlo parameter values have already been set.") 2017 2018 # Loop over the minimisation stats objects. 2019 for object_name in min_names: 2020 # Name for the simulation object. 2021 sim_object_name = object_name + '_sim' 2022 2023 # Test if the simulation object already exists. 2024 if hasattr(cdp, sim_object_name): 2025 raise RelaxError("Monte Carlo parameter values have already been set.") 2026 2027 # Spin specific parameters. 2028 if model_type != 'diff': 2029 for spin in spin_loop(): 2030 # Skip deselected spins. 2031 if not spin.select: 2032 continue 2033 2034 # Loop over all the parameter names. 2035 for object_name in param_names: 2036 # Name for the simulation object. 2037 sim_object_name = object_name + '_sim' 2038 2039 # Test if the simulation object already exists. 2040 if hasattr(spin, sim_object_name): 2041 raise RelaxError("Monte Carlo parameter values have already been set.") 2042 2043 2044 # Set the Monte Carlo parameter values. 2045 ####################################### 2046 2047 # Loop over the global minimisation stats objects. 2048 for object_name in min_names: 2049 # Skip non-existent objects. 2050 if not hasattr(cdp, object_name): 2051 continue 2052 2053 # Name for the simulation object. 2054 sim_object_name = object_name + '_sim' 2055 2056 # Create the simulation object. 2057 setattr(cdp, sim_object_name, []) 2058 2059 # Get the simulation object. 2060 sim_object = getattr(cdp, sim_object_name) 2061 2062 # Loop over the simulations. 2063 for j in range(cdp.sim_number): 2064 # Get the object. 2065 object = getattr(cdp, object_name) 2066 2067 # Copy and append the data. 2068 sim_object.append(deepcopy(object)) 2069 2070 # Diffusion tensor parameters and non spin specific minimisation statistics. 2071 if model_type == 'diff' or model_type == 'all': 2072 # Set up the number of simulations. 2073 cdp.diff_tensor.set_sim_num(cdp.sim_number) 2074 2075 # Loop over the parameters, setting the initial simulation values to those of the parameter value. 2076 for object_name in diff_params: 2077 for j in range(cdp.sim_number): 2078 cdp.diff_tensor.set(param=object_name, value=deepcopy(getattr(cdp.diff_tensor, object_name)), category='sim', sim_index=j) 2079 2080 # Spin specific parameters. 2081 if model_type != 'diff': 2082 for spin in spin_loop(): 2083 # Skip deselected spins. 2084 if not spin.select: 2085 continue 2086 2087 # Loop over all the data names. 2088 for object_name in param_names: 2089 # Name for the simulation object. 2090 sim_object_name = object_name + '_sim' 2091 2092 # Create the simulation object. 2093 setattr(spin, sim_object_name, []) 2094 2095 # Get the simulation object. 2096 sim_object = getattr(spin, sim_object_name) 2097 2098 # Loop over the simulations. 2099 for j in range(cdp.sim_number): 2100 # Copy and append the data. 2101 sim_object.append(deepcopy(getattr(spin, object_name))) 2102 2103 # Loop over all the minimisation object names. 2104 for object_name in min_names: 2105 # Name for the simulation object. 2106 sim_object_name = object_name + '_sim' 2107 2108 # Create the simulation object. 2109 setattr(spin, sim_object_name, []) 2110 2111 # Get the simulation object. 2112 sim_object = getattr(spin, sim_object_name) 2113 2114 # Loop over the simulations. 2115 for j in range(cdp.sim_number): 2116 # Copy and append the data. 2117 sim_object.append(deepcopy(getattr(spin, object_name)))
2118 2119
2120 - def sim_return_chi2(self, index=None, model_info=None):
2121 """Return the simulation chi-squared values. 2122 2123 @keyword index: The optional simulation index. 2124 @type index: int 2125 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 2126 @type model_info: int 2127 @return: The list of simulation chi-squared values. If the index is supplied, only a single value will be returned. 2128 @rtype: list of float or float 2129 """ 2130 2131 # Determine the model type. 2132 model_type = determine_model_type() 2133 2134 # Single instance. 2135 if model_type == 'all' or model_type == 'diff': 2136 return cdp.chi2_sim 2137 2138 # Multiple instances. 2139 else: 2140 # Get the spin container. 2141 spin = return_spin_from_index(model_info) 2142 2143 # Return the list. 2144 return spin.chi2_sim
2145 2146
2147 - def sim_return_param(self, index, model_info=None):
2148 """Return the array of simulation parameter values. 2149 2150 @param index: The index of the parameter to return the array of values for. 2151 @type index: int 2152 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 2153 @type model_info: int 2154 @return: The array of simulation parameter values. 2155 @rtype: list of float 2156 """ 2157 2158 # Parameter increment counter. 2159 inc = 0 2160 2161 # Determine the model type. 2162 model_type = determine_model_type() 2163 2164 # Get the parameter object names. 2165 param_names = self.data_names(set='params', scope='spin') 2166 2167 2168 # Diffusion tensor parameters. 2169 ############################## 2170 2171 if model_type == 'diff' or model_type == 'all': 2172 # Spherical diffusion. 2173 if cdp.diff_tensor.type == 'sphere': 2174 # Return the parameter array. 2175 if index == 0: 2176 return cdp.diff_tensor.tm_sim 2177 2178 # Increment. 2179 inc = inc + 1 2180 2181 # Spheroidal diffusion. 2182 elif cdp.diff_tensor.type == 'spheroid': 2183 # Return the parameter array. 2184 if index == 0: 2185 return cdp.diff_tensor.tm_sim 2186 elif index == 1: 2187 return cdp.diff_tensor.Da_sim 2188 elif index == 2: 2189 return cdp.diff_tensor.theta_sim 2190 elif index == 3: 2191 return cdp.diff_tensor.phi_sim 2192 2193 # Increment. 2194 inc = inc + 4 2195 2196 # Ellipsoidal diffusion. 2197 elif cdp.diff_tensor.type == 'ellipsoid': 2198 # Return the parameter array. 2199 if index == 0: 2200 return cdp.diff_tensor.tm_sim 2201 elif index == 1: 2202 return cdp.diff_tensor.Da_sim 2203 elif index == 2: 2204 return cdp.diff_tensor.Dr_sim 2205 elif index == 3: 2206 return cdp.diff_tensor.alpha_sim 2207 elif index == 4: 2208 return cdp.diff_tensor.beta_sim 2209 elif index == 5: 2210 return cdp.diff_tensor.gamma_sim 2211 2212 # Increment. 2213 inc = inc + 6 2214 2215 2216 # Model-free parameters for the model type 'all'. 2217 ################################################# 2218 2219 if model_type == 'all': 2220 # Loop over the spins. 2221 for spin in spin_loop(): 2222 # Skip deselected spins. 2223 if not spin.select: 2224 continue 2225 2226 # Loop over the spin specific parameters. 2227 for param in param_names: 2228 # Return the parameter array. 2229 if index == inc: 2230 return getattr(spin, param + "_sim") 2231 2232 # Increment. 2233 inc = inc + 1 2234 2235 2236 # Model-free parameters for the model types 'mf' and 'local_tm'. 2237 ################################################################ 2238 2239 if model_type == 'mf' or model_type == 'local_tm': 2240 # Get the spin container. 2241 spin = return_spin_from_index(model_info) 2242 2243 # Skip deselected spins. 2244 if not spin.select: 2245 return 2246 2247 # Loop over the spin specific parameters. 2248 for param in param_names: 2249 # Return the parameter array. 2250 if index == inc: 2251 return getattr(spin, param + "_sim") 2252 2253 # Increment. 2254 inc = inc + 1
2255 2256
2257 - def sim_return_selected(self, model_info=None):
2258 """Return the array of selected simulation flags for the spin. 2259 2260 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 2261 @type model_info: int 2262 @return: The array of selected simulation flags. 2263 @rtype: list of int 2264 """ 2265 2266 # Determine the model type. 2267 model_type = determine_model_type() 2268 2269 # Single instance. 2270 if model_type == 'all' or model_type == 'diff': 2271 return cdp.select_sim 2272 2273 # Multiple instances. 2274 else: 2275 # Get the spin container. 2276 spin = return_spin_from_index(model_info) 2277 2278 # Skip if deselected. 2279 if not spin.select: 2280 return 2281 2282 # Return the list. 2283 return spin.select_sim
2284 2285
2286 - def skip_function(self, model_info=None):
2287 """Skip certain data. 2288 2289 @keyword model_info: The model information from model_loop(). This index is zero for the global models or equal to the global spin index (which covers the molecule, residue, and spin indices). 2290 @type model_info: int 2291 @return: True if the data should be skipped, False otherwise. 2292 @rtype: bool 2293 """ 2294 2295 # Determine the model type. 2296 model_type = determine_model_type() 2297 2298 # Sequence specific data. 2299 if (model_type == 'mf' or model_type == 'local_tm') and not return_spin_from_index(model_info).select: 2300 return True 2301 2302 # Don't skip. 2303 return False
2304 2305 2306
2307 -class Data_container:
2308 """Empty class to be used for data storage."""
2309