Checking if generated complexes match any product pattern of the rule.

NFsim_v1.11/src/NFcore/NFcore.hh

@@ -288,6 +288,9 @@ namespace NFcore
 			bool hasVCellAnchors() { return this->bHasVCellAnchors; };
 			bool setVCellAnchors( bool val ) { this->bHasVCellAnchors = val; };
 
+			bool isCheckingProductMaching() { return this->bCheckingProductMatching; };
+			bool setCheckingProductMaching( bool val ) { this->bCheckingProductMatching = val; };
+
 			void setEvaluateComplexScopedLocalFunctions( bool val ) { evaluateComplexScopedLocalFunctions = val; };
 			bool getEvaluateComplexScopedLocalFunctions( ) const { return evaluateComplexScopedLocalFunctions; };
 
@@ -417,6 +420,7 @@ namespace NFcore
 			bool useComplex;     /*!< sets whether or not to dynamically track complexes */
 			bool useVCellCompartments;     /*!< sets whether or not to apply vcell compartment propagation */
 			bool bHasVCellAnchors;     /*!< sets whether this models has moleculeTpyes with anchors */
+			bool bCheckingProductMatching;   /*!< sets whether or not we verify that generated complexes match product patterns */
 			bool useBinaryOutput; /*!< set to true to turn on binary output of data */
 			bool evaluateComplexScopedLocalFunctions; /*!< set to true to turn on enable complex-scoped local functions */
 			int universalTraversalLimit; /*!< sets depth to traverse molecules when updating reactant lists */

NFsim_v1.11/src/NFcore/reactionClass.cpp

@@ -455,44 +455,46 @@ void ReactionClass::fire(double random_A_number) {
 		}
 	}
 
-	for ( complexIter = productComplexes.begin(); complexIter != productComplexes.end(); ++complexIter ) {
-		// update all species observables for this complex
-		Complex* c = *complexIter;
-		set <string> complexMoleculeNames;
-		list <Molecule *>::iterator itm;
-		for( itm = c->complexMembers.begin(); itm != c->complexMembers.end(); itm++ )
-		{
-			string moleculeName = (*itm)->getMoleculeTypeName();
-			complexMoleculeNames.insert(moleculeName);
-		}
-
-		// see if any product pattern is satisfied by this complex
-		list<set<string>> listOfProductSets = this->transformationSet->getProductSets();
-		list<set<string>>::iterator it;
-		bool matchingSetFound = false;
-		for (it=listOfProductSets.begin(); it != listOfProductSets.end(); it++) {
-			set<string> productMoleculeNames = *it;
-			set<string>::iterator it;
-			bool setNamesFound = true;
-			for (it=productMoleculeNames.begin(); it != productMoleculeNames.end(); it++) {
-				string moleculeName = *it;
-				if ( complexMoleculeNames.find(moleculeName) == complexMoleculeNames.end() ) {
-					setNamesFound = false;
-					break;		// a molecule in this product pattern is not in the complex
-				}
-			}
-			if(setNamesFound == true) {
-				matchingSetFound = true; // all molecules in the product pattern found in the complex
-				break;					 // look no more
-			}
-		}
-		complexMoleculeNames.clear();
+    if (this->system->isCheckingProductMaching()) {
+        for (complexIter = productComplexes.begin(); complexIter != productComplexes.end(); ++complexIter) {
+            // update all species observables for this complex
+            Complex *c = *complexIter;
+            set<string> complexMoleculeNames;
+            list<Molecule *>::iterator itm;
+            for (itm = c->complexMembers.begin(); itm != c->complexMembers.end(); itm++) {
+                string moleculeName = (*itm)->getMoleculeTypeName();
+                complexMoleculeNames.insert(moleculeName);
+            }
+
+            // see if any product pattern is satisfied by this complex
+            list<set<string>> listOfProductSets = this->transformationSet->getProductSets();
+            list<set<string>>::iterator it;
+            bool matchingSetFound = false;
+            for (it = listOfProductSets.begin(); it != listOfProductSets.end(); it++) {
+                set<string> productMoleculeNames = *it;
+                set<string>::iterator it;
+                bool setNamesFound = true;
+                for (it = productMoleculeNames.begin(); it != productMoleculeNames.end(); it++) {
+                    string moleculeName = *it;
+                    if (complexMoleculeNames.find(moleculeName) == complexMoleculeNames.end()) {
+                        setNamesFound = false;
+                        break;        // a molecule in this product pattern is not in the complex
+                    }
+                }
+                if (setNamesFound == true) {
+                    matchingSetFound = true; // all molecules in the product pattern found in the complex
+                    break;                     // look no more
+                }
+            }
+            complexMoleculeNames.clear();
 //		cout<<endl;
-		if(matchingSetFound == false) {	 // this complex matches no product pattern, it shouldn't exist
-			cerr << "\nReaction '"<< name <<"' produced a complex that does not match any product pattern.\n"<< endl;
-			exit(1);
-		}
-	}
+            if (matchingSetFound == false) {     // this complex matches no product pattern, it shouldn't exist
+                cerr << "\n\nReaction '" << name << "' produced a complex that does not match any product pattern.\n"
+                     << endl;
+                exit(1);
+            }
+        }
+    }
 
 		// If we're handling observables on the fly, tell each molecule to add itself to observables.
 	if (onTheFlyObservables) {

NFsim_v1.11/src/NFcore/system.cpp

@@ -26,6 +26,7 @@ System::System(string name)
 	nextReaction = 0;
 	this->useVCellCompartments = false;
 	this->bHasVCellAnchors = false;
+	this->bCheckingProductMatching = false;
 	this->useComplex = false;     // NETGEN -- is this needed?
 	// NETGEN
 	allComplexes.setSystem( this );
@@ -54,6 +55,7 @@ System::System(string name, bool useComplex)
 
 	this->useVCellCompartments = false;
 	this->bHasVCellAnchors = false;
+	this->bCheckingProductMatching = false;
 	this->useComplex = useComplex;    // NETGEN -- is this needed?
 	// NETGEN
 	allComplexes.setSystem( this );
@@ -81,6 +83,7 @@ System::System(string name, bool useComplex, int globalMoleculeLimit)
 	nextReaction = 0;
 	this->useVCellCompartments = false;
 	this->bHasVCellAnchors = false;
+	this->bCheckingProductMatching = false;
 	this->useComplex = useComplex;  // NETGEN -- is this needed?
 	// NETGEN
 	allComplexes.setSystem( this );

NFsim_v1.11/src/NFsim.cpp

@@ -89,6 +89,8 @@
  *           cleared for each product complex (after assigning a consistent location for
  *           althrow std::runtime_error("error in NFutil::convertToDouble(\"" + s + "\")");l molecules in a complex.
  *
+ *	-pcmatch = tries to match generated complexes to the product patterns of the rule
+ *
  *  \section devel_sec Developers
  * To begin developing and extending NFsim, the best place to start looking is in
  * the src/NFtest/simple_system directory. Here you'll find two files, simple_system.hh
@@ -401,6 +403,14 @@ System *initSystemFromFlags(map<string,string> argMap, bool verbose)
 					if(verbose) cout<<"\tVCell Compartment trick (-vcell) flag detected."<<endl<<endl;
 				}
 
+				// Check if the complexes generated by applying a rule do match the product patterns
+				if (argMap.find("pcmatch")!=argMap.end()) {
+					s->setCheckingProductMaching(true);
+					if(verbose) cout<<"\tEnable (-pcmatch) the matching of complex with product pattern."<<endl<<endl;
+				} else {
+					s->setCheckingProductMaching(false);
+				}
+
 				//If requested, be sure to output the values of global functions
 				if (argMap.find("ogf")!=argMap.end()) {
 					s->turnOnGlobalFuncOut();