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

.gitignore

@@ -15,3 +15,4 @@
 /cmake-build*
 /build*
 /.idea/
+/nfsim/

NFsim_v1.11/src/NFcore/reactionClass.cpp

@@ -455,7 +455,46 @@ void ReactionClass::fire(double random_A_number) {
 		}
 	}
 
-	// If we're handling observables on the fly, tell each molecule to add itself to observables.
+	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 we're handling observables on the fly, tell each molecule to add itself to observables.
 	if (onTheFlyObservables) {
 
 		// molecule observables..

NFsim_v1.11/src/NFinput/NFinput.cpp

@@ -1113,7 +1113,6 @@ bool NFinput::initReactionRules(
 				vector <TemplateMolecule *> addmol_templates;
 
 
-
 				//  Read in the Reactant Patterns for this rule
 				TiXmlElement *pListOfReactantPatterns = pRxnRule->FirstChildElement("ListOfReactantPatterns");
 				if(!pListOfReactantPatterns) {
@@ -1149,7 +1148,51 @@ bool NFinput::initReactionRules(
 				//		cout << (*it).first << " => " << (*it).second->getMoleculeType()->getName() << endl;
 
 
-				///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+                //  Read in the Product Patterns for this rule
+                list <set<string>> listOfProductSets;
+                TiXmlElement *pListOfProductPatterns = pRxnRule->FirstChildElement("ListOfProductPatterns");
+                if(pListOfProductPatterns != nullptr) {
+                    TiXmlElement *pProduct;
+                    for ( pProduct = pListOfProductPatterns->FirstChildElement("ProductPattern"); pProduct != 0; pProduct = pProduct->NextSiblingElement("ProductPattern"))
+                    {
+                        const char *productName = pProduct->Attribute("id");
+                        if(!productName) {
+                            continue;
+                        }
+                        TiXmlElement *pListOfMols = pProduct->FirstChildElement("ListOfMolecules");
+                        if(!pListOfMols) {
+                            continue;
+                        }
+                        TiXmlElement *pMol;
+                        set <string> moleculeNames;
+                        for ( pMol = pListOfMols->FirstChildElement("Molecule"); pMol != 0; pMol = pMol->NextSiblingElement("Molecule"))
+                        {
+                            string moleculeName;
+                            if(!pMol->Attribute("name")) {
+                                continue;
+                            }
+                            moleculeName = pMol->Attribute("name");
+                            if(moleculeName=="Trash" || moleculeName=="TRASH" || moleculeName=="trash") {
+                                continue;   // we ignore this, may be dealt with elsewhere
+                            }
+                            moleculeNames.insert(moleculeName);
+                        }
+                        listOfProductSets.push_back(moleculeNames);
+                    }
+                }
+
+//				list <set<string>>::iterator it;
+//                for (it=listOfProductSets.begin(); it != listOfProductSets.end(); it++) {
+//                    set<string> moleculeNames = *it;
+//                    set<string>::iterator it;
+//                    for (it=moleculeNames.begin(); it != moleculeNames.end(); it++) {
+//                        string moleculeName = *it;
+//                        cout<<moleculeName<<" ";
+//                    }
+//                    cout<<endl;
+//                }
+
+                ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 				//Read in the list of operations we need to perform in this rule
 				TiXmlElement *pListOfOperations = pRxnRule->FirstChildElement("ListOfOperations");
 				if ( !pListOfOperations )
@@ -1364,6 +1407,8 @@ bool NFinput::initReactionRules(
 					ts->setSymmetryFactor(symmetryFactor);
 				}
 
+				ts->addProductSets(listOfProductSets);
+
 
 				//Next extract out the state changes
 				TiXmlElement *pStateChange;

NFsim_v1.11/src/NFreactions/transformations/transformationSet.cpp

@@ -174,6 +174,11 @@ bool TransformationSet::addLocalFunctionReference(TemplateMolecule *t, string Po
 	return true;
 }
 
+bool TransformationSet::addProductSets(list<set<string>> &productSets) {
+	this->productSets = productSets;
+	return true;
+}
+
 
 bool TransformationSet::addIncrementStateTransform(TemplateMolecule *t, string cName)
 {

NFsim_v1.11/src/NFreactions/transformations/transformationSet.hh

@@ -161,7 +161,13 @@ namespace NFcore
 			*/
 			bool addLocalFunctionReference(TemplateMolecule *t, string PointerName, int scope);
 
-
+			/*!
+				Adds lists of molecule names for each product pattern
+				Used to crudely verify that a generated complex matches a product pattern
+				@author Dan Vasilescu
+			*/
+			bool addProductSets(list<set<string>> &productSets);
+			list<set<string>> getProductSets() const { return this->productSets; };
 			/*!
 				Call this (in a ReactionClass) to transform the array of MappingSets (one
 				MappingSet per reactant in the correct position in the array, please!).
@@ -323,6 +329,9 @@ namespace NFcore
 			/*!	The array of TemplateMolecules that represent the reactants */
 			TemplateMolecule ** reactants;
 
+			/*! names of the molecules in every product pattern of the reaction	*/
+			list<set<string>> productSets;
+
 			/*!	The array of TemplateMolecules that represent the added molecules.
 			 *   Not sure if this will be used.  --Justin */
 			TemplateMolecule ** addmol;