Fix issues #63 and #65: implement product-side include/exclude filtering and add species observable dirty-flag cache

src/NFcore/NFcore.hh

@@ -1531,12 +1531,18 @@ namespace NFcore
 			// unset canonical flag
 			void unsetCanonical ( ) { is_canonical = false; };
 
+			// Species observable cache (Issue #65)
+			void setSpeciesObsDirty() { _speciesObsDirty = true; }
+			bool isSpeciesObsDirty() const { return _speciesObsDirty; }
+			int* getSpeciesObsCache() { return _speciesObsCache; }
+			int getSpeciesObsCacheSize() const { return _speciesObsCacheSize; }
+			void ensureSpeciesObsCache(int requiredSize);
+			void clearSpeciesObsDirty() { _speciesObsDirty = false; }
+
 			//This is public so that anybody can access the molecules quickly
 			list <Molecule *> complexMembers;
 			list <Molecule *>::iterator molIter;
 
-
-
 		protected:
 			// generate a canonical label using Nauty
 			void   generateCanonicalLabel ( );
@@ -1547,6 +1553,10 @@ namespace NFcore
 			bool    is_canonical;
 			string  canonical_label;
 
+			int*   _speciesObsCache;
+			int    _speciesObsCacheSize;
+			bool   _speciesObsDirty;
+
 		private:
 
 	};

src/NFcore/complex.cpp

@@ -13,7 +13,8 @@ using namespace NFcore;
 const int Node::IS_MOLECULE = -1;
 
 Complex::Complex(System * s, int ID_complex, Molecule * m)
-	: is_canonical( false ), canonical_label("")
+	: is_canonical( false ), canonical_label(""),
+	  _speciesObsCache(0), _speciesObsCacheSize(0), _speciesObsDirty(true)
 {
 	this->system = s;
 	this->ID_complex = ID_complex;
@@ -22,6 +23,15 @@ Complex::Complex(System * s, int ID_complex, Molecule * m)
 
 Complex::~Complex()
 {
+	delete[] _speciesObsCache;
+}
+
+void Complex::ensureSpeciesObsCache(int requiredSize) {
+	if (_speciesObsCacheSize >= requiredSize) return;
+	delete[] _speciesObsCache;
+	_speciesObsCacheSize = requiredSize;
+	_speciesObsCache = new int[_speciesObsCacheSize];
+	_speciesObsDirty = true;
 }
 
 bool Complex::isAlive() {
@@ -110,6 +120,9 @@ void Complex::mergeWithList(Complex * c)
 	this->unsetCanonical();
 	c->unsetCanonical();
 
+	// invalidate species observable cache
+	this->setSpeciesObsDirty();
+
 	// move molecules in c to this complex
 	c->refactorToNewComplex(this->ID_complex);
     this->complexMembers.splice(complexMembers.end(),c->complexMembers);
@@ -140,6 +153,7 @@ void Complex::updateComplexMembership(Molecule * m)
 	if(m->getComplexID()!=this->ID_complex) { cerr<< "ERROR IN COMPLEX!!! "<<endl; return; }
 
 	unsetCanonical();
+	setSpeciesObsDirty();
 
 	//Get list of things this molecule is still connected to
     list <Molecule *> members;
@@ -163,6 +177,7 @@ void Complex::updateComplexMembership(Molecule * m)
 	//Get the next available complex
 	// NETGEN -- redirected call to ComplexList object at system->allComplexes
 	Complex *newComplex = (system->getAllComplexes()).getNextAvailableComplex();
+	newComplex->setSpeciesObsDirty();
 	//cout<<" forming new complex:  next available: " <<newComplex->getComplexID()<<endl;
 
 	//renumber our complex elements

src/NFcore/molecule.cpp

@@ -267,18 +267,18 @@ bool Molecule::decrementPopulation()
 void Molecule::setComponentState(int cIndex, int newValue)
 {
 	this->component[cIndex]=newValue;
-	if (useComplex)
-		// Need to manually unset canonical flag since we're not calling a Complex method
+	if (useComplex) {
 		getComplex()->unsetCanonical();
-
+		getComplex()->setSpeciesObsDirty();
+	}
 }
 void Molecule::setComponentState(string cName, int newValue) {
 	this->component[this->parentMoleculeType->getCompIndexFromName(cName)]=newValue;
 
-	if (useComplex)
-		// Need to manually unset canonical flag since we're not calling a Complex method
+	if (useComplex) {
 		getComplex()->unsetCanonical();
-
+		getComplex()->setSpeciesObsDirty();
+	}
 }
 
 
@@ -475,8 +475,11 @@ void Molecule::bind(Molecule *m1, int cIndex1, Molecule *m2, int cIndex2)
 			m1->getComplex()->mergeWithList(m2->getComplex());
 		}
 		else
+		{
 			// Need to manually unset canonical flag since we're not calling a Complex method
 			m1->getComplex()->unsetCanonical();
+			m1->getComplex()->setSpeciesObsDirty();
+		}
 	}
 }
 

src/NFcore/reactionClass.cpp

@@ -423,6 +423,14 @@ string ReactionClass::fire(double random_A_number, bool track) {
 	// (excluding new molecules, we'll get those later --Justin)
 	this->transformationSet->getListOfProducts(mappingSet,products,traversalLimit);
 
+	// Check product-side filters (include_products / exclude_products).
+	// If the resulting products don't pass the filter, treat this as a null event.
+	if (!transformationSet->checkProductFilters(products)) {
+		products.clear();
+		++(System::NULL_EVENT_COUNTER);
+		return string("");
+	}
+
 	// Loop through the products (excluding added molecules) and remove from observables
 	if (this->onTheFlyObservables) {
 		std::unordered_set<int> updatedComplexIds;

src/NFcore/system.cpp

@@ -1396,14 +1396,17 @@ void System::recalculateAllObservables() {
 	}
 
 	int match = 0;
+	int nSpeciesObs = (int)speciesObservables.size();
 	Complex * complex;
 	allComplexes.resetComplexIter();
 	while ((complex = allComplexes.nextComplex())) {
 		if (complex->isAlive()) {
-			for (auto obsIter = speciesObservables.begin(); obsIter != speciesObservables.end(); ++obsIter) {
-				match = (*obsIter)->isObservable(complex);
-				for (int k = 0; k < match; k++) (*obsIter)->straightAdd();
+			complex->ensureSpeciesObsCache(nSpeciesObs);
+			for (int i = 0; i < nSpeciesObs; i++) {
+				match = speciesObservables[i]->isObservable(complex);
+				complex->getSpeciesObsCache()[i] = match;
 			}
+			complex->clearSpeciesObsDirty();
 		}
 	}
 }
@@ -1543,31 +1546,30 @@ void System::outputAllObservableCounts(double cSampleTime, int eventCounter)
 		{	(*molTypeIter)->addAllToObservables(); 	}
 
 		int match = 0;
+		int nSpeciesObs = (int)speciesObservables.size();
 
-	  	// NETGEN -- this bit replaces the commented block below
 	  	Complex * complex;
 	  	allComplexes.resetComplexIter();
 	  	while(  (complex = allComplexes.nextComplex()) )
 	  	{
 	  		if( complex->isAlive() )
 	  		{
-	  			for(obsIter = speciesObservables.begin(); obsIter != speciesObservables.end(); obsIter++)
-	  			{
-	  				match = (*obsIter)->isObservable( complex );
-	  				for (int k=0; k<match; k++) (*obsIter)->straightAdd();
-	  			}
+				complex->ensureSpeciesObsCache(nSpeciesObs);
+
+				if (complex->isSpeciesObsDirty()) {
+					for (int i=0; i<nSpeciesObs; i++) {
+						match = speciesObservables[i]->isObservable(complex);
+						complex->getSpeciesObsCache()[i] = match;
+					}
+					complex->clearSpeciesObsDirty();
+				}
+
+				for (int i=0; i<nSpeciesObs; i++) {
+					match = complex->getSpeciesObsCache()[i];
+					for (int k=0; k<match; k++) speciesObservables[i]->straightAdd();
+				}
 	  		}
 	  	}
-		/*
-		for(complexIter = allComplexes.begin(); complexIter != allComplexes.end(); complexIter++) {
-			if((*complexIter)->isAlive()) {
-				for(obsIter = speciesObservables.begin(); obsIter != speciesObservables.end(); obsIter++) {
-					match = (*obsIter)->isObservable((*complexIter));
-					for(int k=0; k<match; k++) (*obsIter)->straightAdd();
-				}
-			}
-		}
-		*/
 	}
 
 

src/NFinput/NFinput.cpp

@@ -1746,14 +1746,83 @@ bool NFinput::initReactionRulePermutation(
 					}
 				}
 
-				if (pRxnRule->FirstChildElement("ListOfExcludeProducts") ||
-					pRxnRule->FirstChildElement("ListOfIncludeProducts"))
+			// Parse ListOfExcludeProducts
+			// Format: <ListOfExcludeProducts> contains <Pattern> children directly,
+			// with 'id' matching the product pattern number.
+			{
+				TiXmlElement *pExcludeProducts;
+				for (pExcludeProducts = pRxnRule->FirstChildElement("ListOfExcludeProducts");
+					 pExcludeProducts != 0; pExcludeProducts = pExcludeProducts->NextSiblingElement("ListOfExcludeProducts"))
 				{
-					cerr << "Error:: ReactionRule " << rxnName
-					     << " uses include_products()/exclude_products(), which are not yet enforced in NFsim." << endl;
-					cerr << "Error:: Aborting to avoid silently incorrect results." << endl;
-					return false;
+					if (!pExcludeProducts->Attribute("id")) {
+						cerr << "Error:: ListOfExcludeProducts in " << rxnName << " has no id attribute!" << endl;
+						return false;
+					}
+					string productId = pExcludeProducts->Attribute("id");
+					int productIndex;
+					try {
+						productIndex = stoi(productId) - 1;
+					} catch (...) {
+						cerr << "Error:: ListOfExcludeProducts id '" << productId << "' is not a valid product index in reaction " << rxnName << endl;
+						return false;
+					}
+
+					for (TiXmlElement *pPat = pExcludeProducts->FirstChildElement("Pattern"); pPat != 0; pPat = pPat->NextSiblingElement("Pattern")) {
+						string patternId = pPat->Attribute("id");
+						TiXmlElement *pListOfMols = pPat->FirstChildElement("ListOfMolecules");
+						if (pListOfMols) {
+							map<string, component> dummyComps, dummySymMap;
+							map<string, TemplateMolecule*> dummyTemplates;
+							TemplateMolecule *tm = readPattern(pListOfMols, s, parameter, allowedStates, patternId, dummyTemplates, dummyComps, dummySymMap, verbose, suggestedTraversalLimit);
+							if (tm != NULL) {
+								ts->addExcludeProduct(productIndex, tm, dummyTemplates);
+							} else {
+								cerr << "Error reading pattern for exclude products in reaction " << rxnName << endl;
+								return false;
+							}
+						}
+					}
+				}
+			}
+
+			// Parse ListOfIncludeProducts
+			// Format: <ListOfIncludeProducts> contains <Pattern> children directly,
+			// with 'id' matching the product pattern number.
+			{
+				TiXmlElement *pIncludeProducts;
+				for (pIncludeProducts = pRxnRule->FirstChildElement("ListOfIncludeProducts");
+					 pIncludeProducts != 0; pIncludeProducts = pIncludeProducts->NextSiblingElement("ListOfIncludeProducts"))
+				{
+					if (!pIncludeProducts->Attribute("id")) {
+						cerr << "Error:: ListOfIncludeProducts in " << rxnName << " has no id attribute!" << endl;
+						return false;
+					}
+					string productId = pIncludeProducts->Attribute("id");
+					int productIndex;
+					try {
+						productIndex = stoi(productId) - 1;
+					} catch (...) {
+						cerr << "Error:: ListOfIncludeProducts id '" << productId << "' is not a valid product index in reaction " << rxnName << endl;
+						return false;
+					}
+
+					for (TiXmlElement *pPat = pIncludeProducts->FirstChildElement("Pattern"); pPat != 0; pPat = pPat->NextSiblingElement("Pattern")) {
+						string patternId = pPat->Attribute("id");
+						TiXmlElement *pListOfMols = pPat->FirstChildElement("ListOfMolecules");
+						if (pListOfMols) {
+							map<string, component> dummyComps, dummySymMap;
+							map<string, TemplateMolecule*> dummyTemplates;
+							TemplateMolecule *tm = readPattern(pListOfMols, s, parameter, allowedStates, patternId, dummyTemplates, dummyComps, dummySymMap, verbose, suggestedTraversalLimit);
+							if (tm != NULL) {
+								ts->addIncludeProduct(productIndex, tm, dummyTemplates);
+							} else {
+								cerr << "Error reading pattern for include products in reaction " << rxnName << endl;
+								return false;
+							}
+						}
+					}
 				}
+			}
 
 
 				//Next extract out the state changes

src/NFreactions/transformations/transformationSet.cpp

@@ -109,6 +109,7 @@ TransformationSet::~TransformationSet()
 	// (SIGABRT/SIGSEGV). rf.pattern is one of these templates, so it must not be
 	// deleted here either.
 	reactantFilters.clear();
+	productFilters.clear();
 
 	delete [] transformations;
 	delete [] reactants;
@@ -1070,3 +1071,57 @@ bool TransformationSet::checkReactantFilters(int reactantIndex, Molecule *mol) c
 	}
 	return true;
 }
+
+void TransformationSet::addExcludeProduct(int productIndex, TemplateMolecule *pattern, const map<string, TemplateMolecule*>& parsedTemplates) {
+	ProductFilter pf;
+	pf.productIndex = productIndex;
+	pf.pattern = pattern;
+	pf.isExclude = true;
+	pf.parsedTemplates = parsedTemplates;
+	productFilters.push_back(pf);
+}
+
+void TransformationSet::addIncludeProduct(int productIndex, TemplateMolecule *pattern, const map<string, TemplateMolecule*>& parsedTemplates) {
+	ProductFilter pf;
+	pf.productIndex = productIndex;
+	pf.pattern = pattern;
+	pf.isExclude = false;
+	pf.parsedTemplates = parsedTemplates;
+	productFilters.push_back(pf);
+}
+
+bool TransformationSet::checkProductFilters(const list<Molecule *> &products) const {
+	if (productFilters.empty()) return true;
+
+	// Collect unique complexes from the product molecule list
+	unordered_set<Complex*> productComplexes;
+	for (Molecule *mol : products) {
+		if (mol == 0 || !mol->isAlive()) continue;
+		productComplexes.insert(mol->getComplex());
+	}
+
+	for (const auto &pf : productFilters) {
+		bool anyComplexMatches = false;
+		for (Complex *c : productComplexes) {
+			bool patternMatches = false;
+			for (Molecule *cm : c->complexMembers) {
+				if (pf.pattern->compare(cm)) {
+					patternMatches = true;
+					break;
+				}
+			}
+			if (patternMatches) {
+				anyComplexMatches = true;
+				break;
+			}
+		}
+
+		if (pf.isExclude && anyComplexMatches) {
+			return false;
+		}
+		if (!pf.isExclude && !anyComplexMatches) {
+			return false;
+		}
+	}
+	return true;
+}

src/NFreactions/transformations/transformationSet.hh

@@ -50,6 +50,13 @@ namespace NFcore
 			map<string, TemplateMolecule*> parsedTemplates;
 		};
 
+		struct ProductFilter {
+			int productIndex;
+			TemplateMolecule *pattern;
+			bool isExclude;
+			map<string, TemplateMolecule*> parsedTemplates;
+		};
+
 		public:
 
 			/*!
@@ -315,6 +322,10 @@ namespace NFcore
 			void addIncludeReactant(int reactantIndex, TemplateMolecule *pattern, const map<string, TemplateMolecule*>& parsedTemplates);
 			bool checkReactantFilters(int reactantIndex, Molecule *mol) const;
 
+			void addExcludeProduct(int productIndex, TemplateMolecule *pattern, const map<string, TemplateMolecule*>& parsedTemplates);
+			void addIncludeProduct(int productIndex, TemplateMolecule *pattern, const map<string, TemplateMolecule*>& parsedTemplates);
+			bool checkProductFilters(const list<Molecule *> &products) const;
+
 		protected:
 			bool addBindingTransformImpl(TemplateMolecule *t1, string bSiteName1, TemplateMolecule *t2, string bSiteName2, bool isNewMolecule);
 
@@ -395,6 +406,7 @@ namespace NFcore
 			vector < pair<int,int> >  collision_pairs;
 
 			vector <ReactantFilter> reactantFilters;
+			vector <ProductFilter> productFilters;
 
 		private:
 			void initCommon();