Author: bugman Date: Wed Mar 26 17:28:19 2014 New Revision: 22535 URL: http://svn.gna.org/viewcvs/relax?rev=22535&view=rev Log: Huge speed up of the interatom.define user function. This is to fix bug #21862 (https://gna.org/bugs/?21862), the freezing up of relax when using the dipolar relaxation button in the model-free auto-analysis in the GUI. This involves a number of changes. The algorithm for the backend of the interatom.define user function has been broken into two separate parts. The first part is new and uses the internal structural object atom_loop() twice for each spin ID string. This then calls the new are_bonded_index() structural object method which uses atom indices to find if two atoms are bonded, as the atom indices are returned from the atom_loop(). The are_bonded_index() is orders of magnitude faster than are_bonded() as selection objects are not used and the bonded data structure can be directly accessed. The are_bonded() method has also been slightly speed up by improving its logic. The second part is to perform the original algorithm of two nested spin loops over each spin ID and using the are_bonded() structural method. This second part only happens if the first part finds nothing. The structural object atom_loop() method has been modified to be able to return the molecule index. These indices are needed for the new are_bonded_index() method. When running relax with the profile flag turned on, a simple script which loads the 'Ubiquitin2.bz2' saved state and then the "interatom.define(spin_id1='@N', spin_id2='@H', direct_bond=True)" user function decreases from a total time of 143 to 3.8 seconds. However there are no speed changes detectable in the relax test suite - on one computer the system, unit and GUI tests only only vary by a fraction of a second. Modified: trunk/lib/structure/internal/object.py trunk/pipe_control/interatomic.py Modified: trunk/lib/structure/internal/object.py URL: http://svn.gna.org/viewcvs/relax/trunk/lib/structure/internal/object.py?rev=22535&r1=22534&r2=22535&view=diff ============================================================================== --- trunk/lib/structure/internal/object.py (original) +++ trunk/lib/structure/internal/object.py Wed Mar 26 17:28:19 2014 @@ -1182,18 +1182,18 @@ # Find the first atom. index1 = None - for i in range(len(mol.atom_num)): - # Skip a non-matching first atom. - if sel_obj1.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): - index1 = i - break - - # Find the second atom. index2 = None for i in range(len(mol.atom_num)): - # Skip a non-matching first atom. - if sel_obj2.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): + # Matching first atom. + if index1 == None and sel_obj1.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): + index1 = i + + # Matching second atom. + if index2 == None and sel_obj2.contains_spin(mol.atom_num[i], mol.atom_name[i], mol.res_num[i], mol.res_name[i], mol.mol_name): index2 = i + + # Nothing left to do. + if index1 != None and index2 != None: break # Connectivities exist. @@ -1204,7 +1204,40 @@ return False - def atom_loop(self, atom_id=None, str_id=None, model_num=None, mol_name_flag=False, res_num_flag=False, res_name_flag=False, atom_num_flag=False, atom_name_flag=False, element_flag=False, pos_flag=False, index_flag=False, ave=False): + def are_bonded_index(self, mol_index1=None, atom_index1=None, mol_index2=None, atom_index2=None): + """Determine if two atoms, given as indices, are directly bonded to each other. + + @keyword mol_index1: The molecule index of the first atom. + @type mol_index1: int + @keyword atom_index1: The index of the first atom. + @type atom_index1: int + @keyword mol_index2: The molecule index of the second atom. + @type mol_index2: int + @keyword atom_index2: The index of the second atom. + @type atom_index2: int + @return: True if the atoms are directly bonded. + @rtype: bool + """ + + # Alias the molecule. + mol1 = self.structural_data[0].mol[mol_index1] + mol2 = self.structural_data[0].mol[mol_index2] + + # Build the connectivities if needed. + if not len(mol1.bonded[atom_index1]): + self._find_bonded_atoms(atom_index1, mol1, radius=2) + if not len(mol2.bonded[atom_index2]): + self._find_bonded_atoms(atom_index2, mol2, radius=2) + + # Is the second atom in the bonded list of the first? + if atom_index2 in mol1.bonded[atom_index1]: + return True + + # No! + return False + + + def atom_loop(self, atom_id=None, str_id=None, model_num=None, mol_name_flag=False, res_num_flag=False, res_name_flag=False, atom_num_flag=False, atom_name_flag=False, element_flag=False, pos_flag=False, mol_index_flag=False, index_flag=False, ave=False): """Generator function for looping over all atoms in the internal relax structural object. This method should be designed as a U{generator<http://www.python.org/dev/peps/pep-0255/>}. It should loop over all atoms of the system yielding the following atomic information, if the corresponding flag is True, in tuple form: @@ -1239,6 +1272,8 @@ @type element_flag: bool @keyword pos_flag: A flag which if True will cause the atomic position to be yielded. @type pos_flag: bool + @keyword mol_index_flag: A flag which if True will cause the molecule index to be yielded. + @type mol_index_flag: bool @keyword index_flag: A flag which if True will cause the atomic index to be yielded. @type index_flag: bool @keyword ave: A flag which if True will result in this method returning the average atom properties across all loaded structures. @@ -1337,6 +1372,8 @@ atomic_tuple = atomic_tuple + (element,) if pos_flag: atomic_tuple = atomic_tuple + (pos,) + if mol_index_flag: + atomic_tuple += (mol_index,) if index_flag: atomic_tuple += (i,) Modified: trunk/pipe_control/interatomic.py URL: http://svn.gna.org/viewcvs/relax/trunk/pipe_control/interatomic.py?rev=22535&r1=22534&r2=22535&view=diff ============================================================================== --- trunk/pipe_control/interatomic.py (original) +++ trunk/pipe_control/interatomic.py Wed Mar 26 17:28:19 2014 @@ -1,6 +1,6 @@ ############################################################################### # # -# Copyright (C) 2012-2013 Edward d'Auvergne # +# Copyright (C) 2012-2014 Edward d'Auvergne # # # # This file is part of the program relax (http://www.nmr-relax.com). # # # @@ -36,7 +36,7 @@ from lib.io import extract_data, strip, write_data from lib.warnings import RelaxWarning, RelaxZeroVectorWarning from pipe_control import pipes -from pipe_control.mol_res_spin import Selection, count_spins, exists_mol_res_spin_data, return_spin, spin_loop +from pipe_control.mol_res_spin import Selection, count_spins, exists_mol_res_spin_data, generate_spin_id_unique, return_spin, spin_loop def copy(pipe_from=None, pipe_to=None, spin_id1=None, spin_id2=None, verbose=True): @@ -216,23 +216,52 @@ # Get the data pipe. dp = pipes.get_pipe(pipe) - # Loop over both spin selections. + # Initialise the spin ID pairs list. ids = [] - for spin1, mol_name1, res_num1, res_name1, id1 in spin_loop(spin_id1, pipe=pipe, full_info=True, return_id=True): - for spin2, mol_name2, res_num2, res_name2, id2 in spin_loop(spin_id2, pipe=pipe, full_info=True, return_id=True): - # Directly bonded atoms. - if direct_bond: - # Different molecules. - if mol_name1 != mol_name2: + + # Use the structural data to find connected atoms. + if hasattr(dp, 'structure'): + # Loop over the atoms of the first spin selection. + for mol_name1, res_num1, res_name1, atom_num1, atom_name1, mol_index1, atom_index1 in dp.structure.atom_loop(atom_id=spin_id1, model_num=1, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, mol_index_flag=True, index_flag=True): + # Generate the first spin ID. + id1 = generate_spin_id_unique(pipe_cont=dp, mol_name=mol_name1, res_num=res_num1, res_name=res_name1, spin_num=atom_num1, spin_name=atom_name1) + + # Do the spin exist? + if not return_spin(id1): + continue + + # Loop over the atoms of the second spin selection. + for mol_name2, res_num2, res_name2, atom_num2, atom_name2, mol_index2, atom_index2 in dp.structure.atom_loop(atom_id=spin_id2, model_num=1, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, mol_index_flag=True, index_flag=True): + # Directly bonded atoms. + if direct_bond: + # Different molecules. + if mol_name1 != mol_name2: + continue + + # Skip non-bonded atom pairs. + if not dp.structure.are_bonded_index(mol_index1=mol_index1, atom_index1=atom_index1, mol_index2=mol_index2, atom_index2=atom_index2): + continue + + # Generate the second spin ID. + id2 = generate_spin_id_unique(pipe_cont=dp, mol_name=mol_name2, res_num=res_num2, res_name=res_name2, spin_num=atom_num2, spin_name=atom_name2) + + # Do the spin exist? + if not return_spin(id2): continue - # From structural info. - if hasattr(dp, 'structure') and dp.structure.get_molecule(mol_name1, model=1): - if not dp.structure.are_bonded(atom_id1=id1, atom_id2=id2): + # Store the IDs for the printout. + ids.append([id1, id2]) + + # No structural data present or the spin IDs are not in the structural data, so use spin loops and some basic rules. + if ids == []: + for spin1, mol_name1, res_num1, res_name1, id1 in spin_loop(spin_id1, pipe=pipe, full_info=True, return_id=True): + for spin2, mol_name2, res_num2, res_name2, id2 in spin_loop(spin_id2, pipe=pipe, full_info=True, return_id=True): + # Directly bonded atoms. + if direct_bond: + # Different molecules. + if mol_name1 != mol_name2: continue - # From the residue info. - else: # No element info. if not hasattr(spin1, 'element'): raise RelaxError("The spin '%s' does not have the element type set." % id1) @@ -252,22 +281,8 @@ elif pair and res_num1 == None and res_name1 != res_name2: continue - # Get the interatomic data object, if it exists. - interatom = return_interatom(id1, id2, pipe=pipe) - - # Create the container if needed. - if interatom == None: - interatom = create_interatom(spin_id1=id1, spin_id2=id2, pipe=pipe) - - # Check that this has not already been set up. - if interatom.dipole_pair: - raise RelaxError("The magnetic dipole-dipole interaction already exists between the spins '%s' and '%s'." % (id1, id2)) - - # Set a flag indicating that a dipole-dipole interaction is present. - interatom.dipole_pair = True - - # Store the IDs for the printout. - ids.append([repr(id1), repr(id2)]) + # Store the IDs for the printout. + ids.append([id1, id2]) # No matches, so fail! if not len(ids): @@ -289,8 +304,30 @@ else: raise RelaxError("Unknown error.") - # Print out. + # Define the interaction. + for id1, id2 in ids: + # Get the interatomic data object, if it exists. + interatom = return_interatom(id1, id2, pipe=pipe) + + # Create the container if needed. + if interatom == None: + interatom = create_interatom(spin_id1=id1, spin_id2=id2, pipe=pipe) + + # Check that this has not already been set up. + if interatom.dipole_pair: + raise RelaxError("The magnetic dipole-dipole interaction already exists between the spins '%s' and '%s'." % (id1, id2)) + + # Set a flag indicating that a dipole-dipole interaction is present. + interatom.dipole_pair = True + + # Printout. if verbose: + # Conversion. + for i in range(len(ids)): + ids[i][0] = repr(ids[i][0]) + ids[i][1] = repr(ids[i][1]) + + # The printout. print("Interatomic interactions are now defined for the following spins:\n") write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)