Search polymer neighborhood of all specified comps

For my polymer neighborhood search, I was previously looking at only the
molecule components that got changed by a transformation. This results
in searching of only polymer neighborhoods of the changed component, but
not other unchanged components which might nevertheless be constrained
by the reaction rule.
So I now look at all molecule components that are constrained as well.
This is done by looking at the template molecule corresponding to the
transformation reactant.

Alter search to only polymeric constrained comps and
Search neighborhood of only bonded polymeric comps

Trying to speed up the polymer search but not very useful

This is the commit message msneddon#2:

This is done to avoid slow searching of all specified components
of all molecules.

README.org

@@ -1,5 +1,9 @@
 #+TITLE: Notes for Rasi's modification to the program
 
+[2018-11-23 Fri]
+- For my polymer neighborhood search, I was previously looking at only the molecule components that got changed by a transformation. This results in searching of only polymer neighborhoods of the changed component, but not other unchanged components which might nevertheless be constrained by the reaction rule.
+- So I now look at all molecule components that are constrained as well. This is done by looking at the template molecule corresponding to the transformation reactant.
+
 [2018-08-23 Thu]
 - NFsim is working robustly now with or without polymer and network connectivity options.
 - As far as I can tell, both these modules are working as expected.

src/NFcore/NFcore.hh

@@ -960,10 +960,10 @@ namespace NFcore
 			void depthFirstSearch(vector <Molecule *> &members);
 
 			/* functions for searching a polymer molecule within the interaction distance
-			 * of the site where they  will undergo a change
+			 * of the site that is specified or changed
 			 * Arvind Rasi Subramaniam
 			 */
-			void traversePolymerNeighborhood(vector <Molecule *> &members, Mapping * mapping);
+			void traversePolymerNeighborhood(vector <Molecule *> &members, int cIndex);
 			/* This function is essentially same as breadthFirstSearch but adapted for
 			 * product retrieval of non-polymer molecules.
 			 * Arvind Rasi Subramaniam

src/NFcore/molecule.cpp

@@ -416,21 +416,22 @@ void Molecule::printBondDetails(ostream &o) {
  */
 void Molecule::printBondDetails(NFstream &o)
 {
-	int degree = 0;
 	o<< parentMoleculeType->getName() << "\t"<<ID_unique;
-	o<<"\t";
-	for(int c=0; c<numOfComponents; c++)
-	{
-		if(bond[c]==NOBOND) {continue;}
-		else {
-			o<<"||";
-			o << parentMoleculeType->getComponentName(c);
-			if (parentMoleculeType->getComponentStateName(c,component[c]) != "NO_STATE") {
-				o<< "-" << parentMoleculeType->getComponentStateName(c,component[c]);
+	if (this->getMoleculeType()->getSystem()->getTrackConnected()) {
+		o<<"\t";
+		for(int c=0; c<numOfComponents; c++)
+		{
+			if(bond[c]==NOBOND) {continue;}
+			else {
+				o<<"||";
+				o << parentMoleculeType->getComponentName(c);
+				if (parentMoleculeType->getComponentStateName(c,component[c]) != "NO_STATE") {
+					o<< "-" << parentMoleculeType->getComponentStateName(c,component[c]);
+				}
+				o<<":";
+				o<<bond[c]->getMoleculeType()->getComponentName(this->indexOfBond[c]);
+				o<<"-"<<bond[c]->getMoleculeTypeName()<<"_"<<bond[c]->getUniqueID();
 			}
-			o<<":";
-			o<<bond[c]->getMoleculeType()->getComponentName(this->indexOfBond[c]);
-			o<<"-"<<bond[c]->getMoleculeTypeName()<<"_"<<bond[c]->getUniqueID();
 		}
 	}
 	o.flush();
@@ -804,17 +805,15 @@ void Molecule::printMoleculeList(vector <Molecule *> &members)
  * @param mapping - contains the molecule component that changed.
  * @author Arvind Rasi Subramaniam
  */
-void Molecule::traversePolymerNeighborhood(vector <Molecule *> &members, Mapping * mapping) {
+void Molecule::traversePolymerNeighborhood(vector <Molecule *> &members, int cIndex) {
 
-	int cIndex;
 	int nearbycIndex;
 	int polymerType;
 	int polymerLocation;
 	int polymerInteractionDistance;
 
 	Molecule * mol = this;
 	MoleculeType * mt =  this->getMoleculeType();
-	cIndex = mapping->getIndex();
 	// The molecule has no bonds to check, for eg. dna in rasi's translation model
 	if (mol->numOfComponents == 0 | cIndex == -1) return;
 	// if molecule is not a polymer, check the binding partner

src/NFcore/templateMolecule.cpp

@@ -451,6 +451,10 @@ void TemplateMolecule::addBond(string thisBsiteName,
 	//If we called this, then we are adding a bond to a nonsymmetric site
 	int compIndex=moleculeType->getCompIndexFromName(thisBsiteName);
 
+	// Add only polymeric components for getting bonded neighbors during
+	// polymeric search
+	// Arvind Rasi Subramaniam Nov 24, 2018
+	if (moleculeType->getPolymerType(compIndex) > 0) specifiedComps.push_back(compIndex);
 	//Insert the new information
 	bondComp.push_back(compIndex);
 	bondCompName.push_back(thisBsiteName);

src/NFcore/templateMolecule.hh

@@ -140,6 +140,10 @@ namespace NFcore
 
 		bool isMoleculeTypeAndComponentPresent(MoleculeType * mt, int cIndex);
 
+		// To get all components that are specified
+		// Includes Bound, Empty, Constrained, Excluded and Bonded Comps
+		vector <int> getAllSpecifiedComps() {return specifiedComps;};
+
 	protected:
 
 		static int TotalTemplateMoleculeCount;
@@ -157,6 +161,12 @@ namespace NFcore
 		////  be handled separately...
 		////  1) unique components
 		////  2) symmetric components
+		///
+
+		// List of all components that are specified in the TemplateMolecule
+		// Includes Bound, Empty, Constrained, Excluded and Bonded Comps
+		// Arvind Rasi Subramaniam Nov 23, 2018
+		vector <int> specifiedComps;
 
 
 		// Which of the unique components must be empty (no bonds)

src/NFreactions/transformations/transformationSet.cpp

@@ -745,7 +745,26 @@ bool TransformationSet::getListOfProducts(MappingSet **mappingSets,
 				// If not, you can accidentally miss searching within the interaction distance
 				// of all site changes.
 				// Arvind Rasi Subramaniam
-				molecule->traversePolymerNeighborhood(products, mappingSets[r]->get(i));
+				molecule->traversePolymerNeighborhood(products, mappingSets[r]->get(i)->getIndex());
+
+				// Iterate over all constrained components of the corresponding TemplateMolecule
+				// even if these components are not changed by the transformation.
+				// This helps in the case of rxns such as endocleave_3_hit to update
+				// ribosomes that are bound but 3' to the cut site.
+				// In that case, the cut site and the ribosome binding site are on
+				// different polymers, and the ribosome binding site does not change,
+				// but it is constrained.
+				// Arvind Rasi Subramaniam Nov 23, 2018
+				TemplateMolecule * tm = transformations[r].at(i)->getTemplateMolecule();
+				// If there are no templatemolecules associated with the transformation
+				if (!tm) continue;
+				// Need to check neighborhood of only polymer molecules
+				if (!tm->getMoleculeType()->checkIfPolymer()) continue;
+				for (int cIndex : tm->getAllSpecifiedComps()) {
+					molecule->traversePolymerNeighborhood(products, cIndex);
+				}
+
+
 			}
 		}
 	}